常用算法之：1、最小二乘法（2）

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

c++实现代码如下：

[cpp]  view plain  copy

1. #include <iostream>  
2. #include <vector>  
3. #include <cmath>  
4. using namespace std;  
5.   
6. //最小二乘拟合相关函数定义  
7. double sum(vector<double> Vnum, int n);  
8. double MutilSum(vector<double> Vx, vector<double> Vy, int n);  
9. double RelatePow(vector<double> Vx, int n, int ex);  
10. double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex);  
11. void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[]);  
12. void CalEquation(int exp, double coefficient[]);  
13. double F(double c[],int l,int m);  
14. double Em[6][4];  
15.   
16. //主函数，这里将数据拟合成二次曲线  
17. int main(int argc, char* argv[])  
18. {  
19.     double arry1[5]={0,0.25,0,5,0.75};  
20.     double arry2[5]={1,1.283,1.649,2.212,2.178};  
21.     double coefficient[5];  
22.     memset(coefficient,0,sizeof(double)*5);  
23.     vector<double> vx,vy;  
24.     for (int i=0; i<5; i++)  
25.     {  
26.         vx.push_back(arry1[i]);  
27.         vy.push_back(arry2[i]);  
28.     }  
29.     EMatrix(vx,vy,5,3,coefficient);  
30.     printf("拟合方程为：y = %lf + %lfx + %lfx^2 \n",coefficient[1],coefficient[2],coefficient[3]);  
31.     return 0;  
32. }  
33. //累加  
34. double sum(vector<double> Vnum, int n)  
35. {  
36.     double dsum=0;  
37.     for (int i=0; i<n; i++)  
38.     {  
39.         dsum+=Vnum[i];  
40.     }  
41.     return dsum;  
42. }  
43. //乘积和  
44. double MutilSum(vector<double> Vx, vector<double> Vy, int n)  
45. {  
46.     double dMultiSum=0;  
47.     for (int i=0; i<n; i++)  
48.     {  
49.         dMultiSum+=Vx[i]*Vy[i];  
50.     }  
51.     return dMultiSum;  
52. }  
53. //ex次方和  
54. double RelatePow(vector<double> Vx, int n, int ex)  
55. {  
56.     double ReSum=0;  
57.     for (int i=0; i<n; i++)  
58.     {  
59.         ReSum+=pow(Vx[i],ex);  
60.     }  
61.     return ReSum;  
62. }  
63. //x的ex次方与y的乘积的累加  
64. double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex)  
65. {  
66.     double dReMultiSum=0;  
67.     for (int i=0; i<n; i++)  
68.     {  
69.         dReMultiSum+=pow(Vx[i],ex)*Vy[i];  
70.     }  
71.     return dReMultiSum;  
72. }  
73. //计算方程组的增广矩阵  
74. void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[])  
75. {  
76.     for (int i=1; i<=ex; i++)  
77.     {  
78.         for (int j=1; j<=ex; j++)  
79.         {  
80.             Em[i][j]=RelatePow(Vx,n,i+j-2);  
81.         }  
82.         Em[i][ex+1]=RelateMutiXY(Vx,Vy,n,i-1);  
83.     }  
84.     Em[1][1]=n;  
85.     CalEquation(ex,coefficient);  
86. }  
87. //求解方程  
88. void CalEquation(int exp, double coefficient[])  
89. {  
90.     for(int k=1;k<exp;k++) //消元过程  
91.     {  
92.         for(int i=k+1;i<exp+1;i++)  
93.         {  
94.             double p1=0;  
95.   
96.             if(Em[k][k]!=0)  
97.                 p1=Em[i][k]/Em[k][k];  
98.   
99.             for(int j=k;j<exp+2;j++)   
100.                 Em[i][j]=Em[i][j]-Em[k][j]*p1;  
101.         }  
102.     }  
103.     coefficient[exp]=Em[exp][exp+1]/Em[exp][exp];  
104.     for(int l=exp-1;l>=1;l--)   //回代求解  
105.         coefficient[l]=(Em[l][exp+1]-F(coefficient,l+1,exp))/Em[l][l];  
106. }  
107. //供CalEquation函数调用  
108. double F(double c[],int l,int m)  
109. {  
110.     double sum=0;  
111.     for(int i=l;i<=m;i++)  
112.         sum+=Em[l-1][i]*c[i];  
113.     return sum;   
114. }