最小二乘法C#源码转

1. #region 最小二乘法拟合
2. ///<summary>
3. ///用最小二乘法拟合二元多次曲线
4. ///例如y=ax+b
5. ///其中MultiLine将返回a，b两个参数。
6. ///a对应MultiLine[1]
7. ///b对应MultiLine[0]
8. ///</summary>
9. ///<param name="arrX">已知点的x坐标集合</param>
10. ///<param name="arrY">已知点的y坐标集合</param>
11. ///<param name="length">已知点的个数</param>
12. ///<param name="dimension">方程的最高次数</param>
13. public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
14. {
15. int n = dimension + 1; //dimension次方程需要求 dimension+1个 系数
16. double[,] Guass = new double[n, n + 1]; //高斯矩阵 例如：y=a0+a1*x+a2*x*x
17. for (int i = 0; i < n; i++)
18. {
19. int j;
20. for (j = 0; j < n; j++)
21. {
22. Guass[i, j] = SumArr(arrX, j + i, length);
23. }
24. Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
25. }
26. return ComputGauss(Guass, n);
27. }
29. private static 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. }
42. private static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
43. {
44. double s = 0;
45. for (int i = 0; i < length; i++)
46. {
47. if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
48. s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
49. else
50. s = s + 1;
51. }
52. return s;
53. }
55. private static double[] ComputGauss(double[,] Guass, int n)
56. {
57. int i, j;
58. int k, m;
59. double temp;
60. double max;
61. double s;
62. double[] x = new double[n];
64. for (i = 0; i < n; i++) x[i] = 0.0;//初始化
65. for (j = 0; j < n; j++)
66. {
67. max = 0;
68. k = j;
69. for (i = j; i < n; i++)
70. {
71. if (Math.Abs(Guass[i, j]) > max)
72. {
73. max = Guass[i, j];
74. k = i;
75. }
76. }
77. if (k != j)
78. {
79. for (m = j; m < n + 1; m++)
80. {
81. temp = Guass[j, m];
82. Guass[j, m] = Guass[k, m];
83. Guass[k, m] = temp;
84. }
85. }
86. if (0 == max)
87. {
88. // "此线性方程为奇异线性方程"
89. return x;
90. }
91. for (i = j + 1; i < n; i++)
92. {
93. s = Guass[i, j];
94. for (m = j; m < n + 1; m++)
95. {
96. Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);
98. }
99. }
100. }
101. //结束for (j=0;j<n;j++)
103. for (i = n - 1; i >= 0; i--)
104. {
105. s = 0;
106. for (j = i + 1; j < n; j++)
107. {
108. s = s + Guass[i, j] * x[j];
109. }
111. x[i] = (Guass[i, n] - s) / Guass[i, i];
113. }
115. return x;
116. }//返回值是函数的系数
117. #endregion