引用
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