用最小二乘法拟合二元多次曲线

引用
http://blog.sina.com.cn/s/blog_6e51df7f0100thie.html
对代码稍作修改和注释，防止链接失效。
  1         ///<summary>
  2         ///用最小二乘法拟合二元多次曲线
  3         ///例如y=ax+b
  4         ///其中MultiLine将返回a，b两个参数。
  5         ///a对应MultiLine[1]
  6         ///b对应MultiLine[0]
  7         ///</summary>
  8         ///<param name="arrX">已知点的x坐标集合</param>
  9         ///<param name="arrY">已知点的y坐标集合</param>
 10         ///<param name="length">已知点的个数</param>
 11         ///<param name="dimension">方程的最高次数</param>
 12         public  double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
 13         {
 14             int n = dimension + 1;                  //dimension次方程需要求 dimension+1个 系数
 15             double[,] Guass = new double[n, n + 1];      //高斯矩阵 例如：y=a0+a1*x+a2*x*x
 16             for (int i = 0; i < n; i++)
 17             {
 18                 int j;
 19                 for (j = 0; j < n; j++)
 20                 {
 21                     Guass[i, j] = SumArr(arrX, j + i, length);
 22                 }
 23                 Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
 24             }
 25 
 26             return ComputGauss(Guass, n);
 27 
 28         }
 29         private  double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和
 30         {
 31             double s = 0;
 32             for (int i = 0; i < length; i++)
 33             {
 34                 if (arr[i] != 0 || n != 0)
 35                     s = s + Math.Pow(arr[i], n);
 36                 else
 37                     s = s + 1;
 38             }
 39             return s;
 40         }
 41         private  double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
 42         {
 43             double s = 0;
 44             for (int i = 0; i < length; i++)
 45             {
 46                 if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
 47                     s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
 48                 else
 49                     s = s + 1;
 50             }
 51             return s;
 52 
 53         }
 54         private  double[] ComputGauss(double[,] Guass, int n)
 55         {
 56             int i, j;
 57             int k, m;
 58             double temp;
 59             double max;
 60             double s;
 61             double[] x = new double[n];
 62 
 63             for (i = 0; i < n; i++) x[i] = 0.0;//初始化
 64 
 65 
 66             for (j = 0; j < n; j++)
 67             {
 68                 max = 0;
 69 
 70                 k = j;
 71                 for (i = j; i < n; i++)
 72                 {
 73                     if (Math.Abs(Guass[i, j]) > max)
 74                     {
 75                         max = Guass[i, j];
 76                         k = i;
 77                     }
 78                 }
 79 
 80 
 81 
 82                 if (k != j)
 83                 {
 84                     for (m = j; m < n + 1; m++)
 85                     {
 86                         temp = Guass[j, m];
 87                         Guass[j, m] = Guass[k, m];
 88                         Guass[k, m] = temp;
 89 
 90                     }
 91                 }
 92 
 93                 if (0 == max)
 94                 {
 95                     // "此线性方程为奇异线性方程" 
 96 
 97                     return x;
 98                 }
 99 
100 
101                 for (i = j + 1; i < n; i++)
102                 {
103 
104                     s = Guass[i, j];
105                     for (m = j; m < n + 1; m++)
106                     {
107                         Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);
108 
109                     }
110                 }
111 
112 
113             }//结束for (j=0;j<n;j++)
114 
115 
116             for (i = n - 1; i >= 0; i--)
117             {
118                 s = 0;
119                 for (j = i + 1; j < n; j++)
120                 {
121                     s = s + Guass[i, j] * x[j];
122                 }
123 
124                 x[i] = (Guass[i, n] - s) / Guass[i, i];
125 
126             }
127 
128             return x;
129         }//返回值是函数的系数
例如：y=a0+a1*x 返回值则为a0 a1
例如：y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2