最小二乘法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. publicstaticdouble[]MultiLine(double[] arrX,double[] arrY,int length,int dimension)//二元多次线性方程拟合曲线
14. {
15. int n = dimension+1;//dimension次方程需要求 dimension+1个 系数
16. double[,]Guass=newdouble[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. returnComputGauss(Guass, n);
27. }
29. privatestaticdoubleSumArr(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. privatestaticdoubleSumArr(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. privatestaticdouble[]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 =newdouble[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