最小二乘法拟合
/*** 最小二乘法拟合
*
*
* @author Ken转发
*
*/
public class Linest {
/**
* <p>
* 函数功能：最小二乘法曲线拟合
* </p>
*
* @param x
* 实型一维数组，长度为 n 。存放给定 n 个数据点的 X 坐标
* @param y
* 实型一维数组，长度为 n 。存放给定 n 个数据点的 Y 坐标
* @param n
* 变量。给定数据点的个数
* @param a
* 实型一维数组，长度为 m 。返回 m-1 次拟合多项式的 m 个系数
* @param m
* 拟合多项式的项数，即拟合多项式的最高次数为 m-1. 要求 m<=n 且m<=20。若 m>n 或 m>20
* ，则本函数自动按 m=min{n,20} 处理.
* <p>
* Date:2007-12-25 16:21 PM
* </p>
* @author qingbao-gao
* @return
*/
public static double[] PolyFit(double x[], double y[], int n, double a[],
int m) {
int i, j, k;
double z, p, c, g, q = 0, d1, d2;
double[] s = new double;
double[] t = new double;
double[] b = new double;
double[] dt = new double;
for (i = 0; i <= m - 1; i++) {
a = 0.0;
}
if (m > n) {
m = n;
}
if (m > 20) {
m = 20;
}
z = 0.0;
for (i = 0; i <= n - 1; i++) {
z = z + x / (1.0 * n);
}
b = 1.0;
d1 = 1.0 * n;
p = 0.0;
c = 0.0;
for (i = 0; i <= n - 1; i++) {
p = p + (x - z);
c = c + y;
}
c = c / d1;
p = p / d1;
a = c * b;
if (m > 1) {
t = 1.0;
t = -p;
d2 = 0.0;
c = 0.0;
g = 0.0;
for (i = 0; i <= n - 1; i++) {
q = x - z - p;
d2 = d2 + q * q;
c = c + y * q;
g = g + (x - z) * q * q;
}
c = c / d2;
p = g / d2;
q = d2 / d1;
d1 = d2;
a = c * t;
a = c * t + a;
}
for (j = 2; j <= m - 1; j++) {
s = t;
s = -p * t + t;
if (j >= 3)
for (k = j - 2; k >= 1; k--) {
s = -p * t + t - q * b;
}
s = -p * t - q * b;
d2 = 0.0;
c = 0.0;
g = 0.0;
for (i = 0; i <= n - 1; i++) {
q = s;
for (k = j - 1; k >= 0; k--) {
q = q * (x - z) + s;
}
d2 = d2 + q * q;
c = c + y * q;
g = g + (x - z) * q * q;
}
c = c / d2;
p = g / d2;
q = d2 / d1;
d1 = d2;
a = c * s;
t = s;
for (k = j - 1; k >= 0; k--) {
a = c * s + a;
b = t;
t = s;
}
}
dt = 0.0;
dt = 0.0;
dt = 0.0;
for (i = 0; i <= n - 1; i++) {
q = a;
for (k = m - 2; k >= 0; k--) {
q = a + q * (x - z);
}
p = q - y;
if (Math.abs(p) > dt) {
dt = Math.abs(p);
}
dt = dt + p * p;
dt = dt + Math.abs(p);
}
return a;
}
/**
* <p>
* 对X轴数据节点球平均值
* </p>
*
* @param x
* 存储X轴节点的数组
* <p>
* Date:2007-12-25 20:21 PM
* </p>
* @author qingbao-gao
* @return 平均值
*/
public static double ave(double[] x) {
double ave = 0;
double sum = 0;
if (x != null) {
for (int i = 0; i < x.length; i++) {
sum += x;
}
System.out.println("sum-->" + sum);
ave = sum / x.length;
System.out.println("ave" + ave + "x.length" + x.length);
}
return ave;
}
/**
* <p>
* 由X值获得Y值
* </p>
*
* @param x
* 当前X轴输入值,即为预测的月份
* @param xx
* 当前X轴输入值的前X数据点
* @param a
* 存储多项式系数的数组
* @param m
* 存储多项式的最高次数的数组
* <p>
* Date:2007-12-25 PM 20:07
* </p>
* <P>
* Author:qingbao-gao
* </P>
* @return 对应X轴节点值的Y轴值
*/
public static double getY(double x, double[] xx, double[] a, int m) {
double y = 0;
double ave = ave(xx);
double l = 0;
for (int i = 0; i < m; i++) {
l = a;
if (i > 0) {
y += a * Math.pow((x - ave), i);
System.out.println(i + "--|-->" + y + "--a--" + a);
}
System.out.println("a|" + a);
}
System.out.println("l--|" + (l));
return (y + l);
}
// --------------------------------------------测试代码
public static void main(String[] args) {
double[] x = { 200401, 200402, 200403, 200404, 200405, 200406, 200407,
200408, 200409, 2004010, 2004011, 2004012, 200501, 200502,
200503, 200504 };
double[] y = { 51, 51, 53, 53, 54, 55, 57, 60, 63, 64, 66, 66, 69, 71,
72, 75 };
double[] a = new double;
double[] aa = PolyFit(x, y, 16, a, 3);
double yy = 0;
System.out.println("拟合-->" + getY(200505, x, aa, 3));
}
}
页:
[1]