数据曲线拟合，数据趋势判断

最近在弄一些数据趋势判断该方面的东西，下面代码是网友使用最小二乘法做的数据拟合算法的java实现

public class PolyFitForcast {
	public PolyFitForcast() {
	}

	/**
	 * <p>
	 * 函数功能：最小二乘法曲线拟合
	 * </p>
	 * <p>
	 * 方程:Y = a(0) + a(1) * (X - X1)+ a(2) * (X - X1)^2 + ..... .+ a(m) * (X -
	 * X1)^m X1为X的平均数
	 * </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[20];
		double[] t = new double[20];
		double[] b = new double[20];
		double[] dt = new double[3];
		for (i = 0; i <= m - 1; i++) {
			a[i] = 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[i] / (1.0 * n);
		}
		b[0] = 1.0;
		d1 = 1.0 * n;
		p = 0.0;
		c = 0.0;
		for (i = 0; i <= n - 1; i++) {
			p = p + (x[i] - z);
			c = c + y[i];
		}
		c = c / d1;
		p = p / d1;
		a[0] = c * b[0];
		if (m > 1) {
			t[1] = 1.0;
			t[0] = -p;
			d2 = 0.0;
			c = 0.0;
			g = 0.0;
			for (i = 0; i <= n - 1; i++) {
				q = x[i] - z - p;
				d2 = d2 + q * q;
				c = c + y[i] * q;
				g = g + (x[i] - z) * q * q;
			}
			c = c / d2;
			p = g / d2;
			q = d2 / d1;
			d1 = d2;
			a[1] = c * t[1];
			a[0] = c * t[0] + a[0];
		}
		for (j = 2; j <= m - 1; j++) {
			s[j] = t[j - 1];
			s[j - 1] = -p * t[j - 1] + t[j - 2];
			if (j >= 3)
				for (k = j - 2; k >= 1; k--) {
					s[k] = -p * t[k] + t[k - 1] - q * b[k];
				}
			s[0] = -p * t[0] - q * b[0];
			d2 = 0.0;
			c = 0.0;
			g = 0.0;
			for (i = 0; i <= n - 1; i++) {
				q = s[j];
				for (k = j - 1; k >= 0; k--) {
					q = q * (x[i] - z) + s[k];
				}
				d2 = d2 + q * q;
				c = c + y[i] * q;
				g = g + (x[i] - z) * q * q;
			}
			c = c / d2;
			p = g / d2;
			q = d2 / d1;
			d1 = d2;
			a[j] = c * s[j];
			t[j] = s[j];
			for (k = j - 1; k >= 0; k--) {
				a[k] = c * s[k] + a[k];
				b[k] = t[k];
				t[k] = s[k];
			}
		}
		dt[0] = 0.0;
		dt[1] = 0.0;
		dt[2] = 0.0;
		for (i = 0; i <= n - 1; i++) {
			q = a[m - 1];
			for (k = m - 2; k >= 0; k--) {
				q = a[k] + q * (x[i] - z);
			}
			p = q - y[i];
			if (Math.abs(p) > dt[2]) {
				dt[2] = Math.abs(p);
			}
			dt[0] = dt[0] + p * p;
			dt[1] = dt[1] + Math.abs(p);
		}
		return a;
	}// end

	/**
	 * <p>
	 * 对X轴数据节点球平均值
	 * </p>
	 * 
	 * @param x
	 *            存储X轴节点的数组
	 *            <p>
	 *            Date:2007-12-25 20:21 PM
	 *            </p>
	 * @author qingbao-gao
	 * @return 平均值
	 */
	public static double average(double[] x) {
		double ave = 0;
		double sum = 0;
		if (x != null) {
			for (int i = 0; i < x.length; i++) {
				sum += x[i];
			}
			ave = sum / 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>
	 * @return 对应X轴节点值的Y轴值
	 */
	public static double getY(double x, double[] xx, double[] a, int m) {
		double y = 0;
		double ave = average(xx);

		double l = 0;
		for (int i = 0; i < m; i++) {
			l = a[0];
			if (i > 0) {
				y += a[i] * Math.pow((x - ave), i);
			}
		}
		return (y + l);
	}

	public static void main(String[] args) throws Exception{
	    
		
		double[] test_xx = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
				             11, 12, 13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,
				             31,32,33,34,35,36,37,38,39,40,
				             41,42,43,44,45,46,47,48,49,50,
				             51,52,53,54,55,56,57,58,59};// 15
		double[] ss  = {
				9880,10000,9880,10088,9958,9869,9712,9670,9719,9625,9416,
				9238,8962,9034,9375,9555,8162,8589,8119,8112,8092,8105,
				7821,7943,8086,7825,7914,7804,7789,7903,7630,7666,
				7545,7632,7619,7582,7427,7499,7463,7384,7315,7092,
				7123,6911,7024,7016,7009,6929,7093,6972,7008,6689,
				6793,6805,6835,6553,6659,6713,6757,6581};
		
		int mmm = ss.length;
		int m = 3;
		double[] a = new double[test_xx.length];
		double[] aa = PolyFit(test_xx, ss, mmm, a, m);
		for(int i=0;i<mmm;i++){
		   System.out.println("拟合-->" + getY(i, test_xx, aa, m));
		}
	}

}