最小二乘曲线拟合

最小二乘曲线拟合

   在上一篇博客中我们介绍了最小二乘法的原理，以及代码实现的例子。

      http://blog.csdn.net/beijingmake209/article/details/27565125

   本次我们再给出一个程序实现的例子。编译环境VC6.0

  先给出一组需要拟合的数据：

  xx[]=  { 0.995119, 2.001185, 2.999068, 4.001035, 4.999859, 6.004461, 6.999335, 7.999433,

             9.002257, 10.003888, 11.004076, 12.001602, 13.003390, 14.001623, 15.003034,  
             16.002561, 17.003010, 18.003897, 19.002563, 20.003530};
 yy[] = { -7.6200, -2.460, 108.7600, 625.020, 2170.500, 5814.5800,13191.8400,26622.060,

            49230.2200, 85066.5000, 139226.2800, 217970.1400, 328843.8600,480798.4200,

            684310.00, 951499.9800, 1296254.9400, 1734346.6600, 2283552.1200, 2963773.50};

实现代码如下：

[html] view plain copy

print?

1. #include <stdio.h>  
2. #include <stdlib.h>  
3.   
4. #include "math.h"  
5.   
6. void polyfit(int n, double *x, double *y, int poly_n, double p[]);  
7. void gauss_solve(int n,double A[],double x[],double b[]);  
8.   
9. /*==================polyfit(n,x,y,poly_n,a)===================*/  
10. /*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/  
11. /*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/  
12. /*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/  
13. void polyfit(int n,double x[],double y[],int poly_n,double p[])  
14. {  
15.     int i,j;  
16.     double *tempx,*tempy,*sumxx,*sumxy,*ata;  
17.       
18.     tempx = (double *)calloc(n , sizeof(double));  
19.     sumxx = (double *)calloc((poly_n*2+1) , sizeof(double));  
20.     tempy = (double *)calloc(n , sizeof(double));  
21.     sumxy = (double *)calloc((poly_n+1) , sizeof(double));  
22.     ata = (double *)calloc( (poly_n+1)*(poly_n+1) , sizeof(double) );  
23.     for (i=0;i<n;i++)  
24.     {  
25.         tempx[i]=1;  
26.         tempy[i]=y[i];  
27.     }  
28.     for (i=0;i<2*poly_n+1;i++)  
29.     {  
30.         for (sumxx[i]=0,j=0;j<n;j++)  
31.         {  
32.             sumxx[i]+=tempx[j];  
33.             tempx[j]*=x[j];  
34.         }  
35.     }  
36.   
37.     for (i=0;i<poly_n+1;i++)  
38.     {  
39.         for (sumxy[i]=0,j=0;j<n;j++)  
40.         {  
41.             sumxy[i]+=tempy[j];  
42.             tempy[j]*=x[j];  
43.         }  
44.     }  
45.   
46.     for (i=0;i<poly_n+1;i++)  
47.     {  
48.         for (j=0;j<poly_n+1;j++)  
49.         {  
50.             ata[i*(poly_n+1)+j]=sumxx[i+j];  
51.         }  
52.     }  
53.     gauss_solve(poly_n+1,ata,p,sumxy);  
54.       
55.     free(tempx);  
56.     free(sumxx);  
57.     free(tempy);  
58.     free(sumxy);  
59.     free(ata);  
60. }  
61.   
62. /*============================================================*/  
63. ////    高斯消元法计算得到   n 次多项式的系数  
64. ////    n: 系数的个数  
65. ////    ata: 线性矩阵  
66. ////    sumxy: 线性方程组的Y值  
67. ////    p: 返回拟合的结果  
68. /*============================================================*/  
69. void gauss_solve(int n,double A[],double x[],double b[])  
70. {  
71.     int i,j,k,r;  
72.     double max;  
73.     for (k=0;k<n-1;k++)  
74.     {  
75.         max=fabs(A[k*n+k]);                 // find maxmum   
76.         r=k;  
77.         for (i=k+1;i<n-1;i++)  
78.         {  
79.             if (max<fabs(A[i*n+i]))  
80.             {  
81.                 max=fabs(A[i*n+i]);  
82.                 r=i;  
83.             }  
84.         }  
85.         if (r!=k)  
86.         {  
87.             for (i=0;i<n;i++)        //change array:A[k]&A[r]  
88.             {  
89.                 max=A[k*n+i];  
90.                 A[k*n+i]=A[r*n+i];  
91.                 A[r*n+i]=max;  
92.             }  
93.   
94.             max=b[k];                    //change array:b[k]&b[r]  
95.             b[k]=b[r];  
96.             b[r]=max;  
97.         }  
98.           
99.         for (i=k+1;i<n;i++)  
100.         {  
101.             for (j=k+1;j<n;j++)  
102.                 A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];  
103.             b[i]-=A[i*n+k]*b[k]/A[k*n+k];  
104.         }  
105.     }   
106.       
107.     for (i=n-1;i>=0;x[i]/=A[i*n+i],i--)  
108.     {  
109.         for (j=i+1,x[i]=b[i];j<n;j++)  
110.             x[i]-=A[i*n+j]*x[j];  
111.     }  
112.   
113. }  
114.   
115. void main()  
116. {  
117.     int i, sizenum;  
118.     double P[6];  
119.     int dimension = 5;  //5次多项式拟合  
120.     //  要拟合的数据  
121.     double xx[]=  {0.995119, 2.001185, 2.999068, 4.001035, 4.999859, 6.004461, 6.999335,     
122.         7.999433, 9.002257, 10.003888, 11.004076, 12.001602, 13.003390, 14.001623, 15.003034,    
123.         16.002561, 17.003010, 18.003897, 19.002563, 20.003530};  
124.     double yy[] = {-7.620000, -2.460000, 108.760000, 625.020000, 2170.500000, 5814.580000,  
125.         13191.840000, 26622.060000, 49230.220000, 85066.500000, 139226.280000, 217970.140000, 328843.860000,  
126.         480798.420000, 684310.000000, 951499.980000, 1296254.940000, 1734346.660000, 2283552.120000, 2963773.500000};  
127.       
128.     sizenum = sizeof(xx)/ sizeof(xx[0]);    //  拟合数据的维数  
129.     polyfit(sizenum, xx, yy, dimension, P);  
130.       
131.     printf("拟合系数, 按升序排列如下:\n");  
132.     for (i=0;i<dimension+1;i++)              //这里是升序排列，Matlab是降序排列  
133.     {  
134.         printf("P[%d]=%lf\n",i,P[i]);  
135.     }  
136. }  

  拟合结果如下所示：

  拟合系数, 按升序排列如下:
 P[0]=-18.544118
 P[1]=6.725933
 P[2]=0.236626
 P[3]=-0.529331
 P[4]=-1.450258
 P[5]=0.999157