最小二乘法怎么做曲线拟合？

//最小二乘法曲线拟合
typedef CArray<double,double>CDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A)
{
//X,Y --  X,Y两轴的坐标
//M   --  结果变量组数
//N   --  采样数目
//A   --  结果参数 register long i,j,k;
double Z,D1,D2,C,P,G,Q;
CDoubleArray B,T,S;
B.SetSize(N);
T.SetSize(N);
S.SetSize(N);
if(M>N)M=N;
for(i=0;i<M;i++)
  (*A)[i]=0;
Z=0;
B[0]=1;
D1=N;
P=0;
C=0;
for(i=0;i<N;i++)
{
  P=P+(*X)[i]-Z;
  C=C+(*Y)[i];
}
C=C/D1;
P=P/D1;
(*A)[0]=C*B[0];
if(M>1)
{
  T[1]=1;
  T[0]=-P;
  D2=0;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=(*X)[i]-Z-P;
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  C=C/D2;
  P=G/D2;
  Q=D2/D1;
  D1=D2;
  (*A)[1]=C*T[1];
  (*A)[0]=C*T[0]+(*A)[0];
}
for(j=2;j<M;j++)
{
  S[j]=T[j-1];
  S[j-1]=-P*T[j-1]+T[j-2];
  if(j>=3)
  {
   for(k=j-2;k>=1;k--)
    S[k]=-P*T[k]+T[k-1]-Q*B[k];
  }
  S[0]=-P*T[0]-Q*B[0];
  D2=0;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=S[j];
   for(k=j-1;k>=0;k--)
    Q=Q*((*X)[i]-Z)+S[k];
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  C=C/D2;
  P=G/D2;
  Q=D2/D1;
  D1=D2;
  (*A)[j]=C*S[j];
  T[j]=S[j];
  for(k=j-1;k>=0;k--)
  {
   (*A)[k]=C*S[k]+(*A)[k];
   B[k]=T[k];
   T[k]=S[k];
  }
}
return TRUE;
}

http://blog.csdn.net/lost_hunter/archive/2006/09/13/1218841.aspx

直线拟合所用的函数的基是{1,x}
曲线拟合可以扩展次数比如二次曲线拟合所用的基可以是{1,x,x^2}.
有兴趣可以找本数值分析的书看，里面有详细介绍。

给出一些离散的点通过一条近似的曲线来表示这些点所代表的函数,这就是曲线拟合的目的.最小二乘法是和种简便而又实用的曲线拟合方法.       它除了能直接拟合形如y=a*x+b形的函数外.还能对1/y=a+b/x,(a>0)y=a*exp(b*x),(a>0)a=a*exp(b/x)(x>0,a>0),y=a+b*lg(x),y=a*x^b(a>0),y=1/(a+b*exp(-x))(a>0);等式子变换成线性表达式后进行拟合.相关的理论请参考相关的数值算法的书籍，我这里只给出关键的函数及主程序段，其余相关的细节就不再一一罗列了.void min2ByFixMethod1(Type* xArr,Type* yArr,int dataNum,FILE* outputFile){       Type* x_xArr,* x_yArr;/*make up data*/       Type c[2][2]={{0,0},{0,0}};/*medium data variable*/       Type d[2]={0,0};       Type tween;/*fatcor's mother*/       Type a,b;/*argu ...*/              assertF(xArr!=NULL,"in min2ByFixMethod1,xArr is NULL");       assertF(yArr!=NULL,"in min2ByFixMethod1,yArr is NULL");       assertF(outputFile!=NULL,"in min2ByFixMethod1,outputFile is NULL");       /*Memory apply*/       x_xArr=(Type*)malloc(sizeof(Type)*dataNum);       x_yArr=(Type*)malloc(sizeof(Type)*dataNum);       assertF(x_xArr!=NULL,"in min2ByFixMethod1,x_xArr is NULL");       assertF(x_yArr!=NULL,"in min2ByFixMethod1,x_yArr is NULL");       /*Other aid data make*/       arrEachBy(xArr,xArr,&x_xArr,dataNum);       arrEachBy(xArr,yArr,&x_yArr,dataNum);              c[0][0]=sumArr(x_xArr,dataNum);       c[0][1]=sumArr(xArr,dataNum);       c[1][0]=sumArr(xArr,dataNum);       c[1][1]=(Type)dataNum;              d[0]=sumArr(x_yArr,dataNum);       d[1]=sumArr(yArr,dataNum);              tween=c[0][0]*c[1][1]- c[0][1]*c[1][0];              assertF(fabs(tween)>=0.00001,"in min2ByFixMethod1,fabs(tween)<0.00001");              a=(d[0]*c[1][1]-c[0][1]*d[1])/tween;       b=(c[0][0]*d[1]-c[1][0]*d[0])/tween;              printf("the fix fun is %f+%fx.\n",b,a);       fprintf(outputFile,"the fix fun is %f+%fx.\r\n",b,a);       }

Type sumArr(Type* inArr,int len){       int i=0;/*iterator num*/       Type ans=0;       assertF(inArr!=NULL,"in sumArr,inArr is NULL\n");       for(i=0;i<len;i++)              ans+=inArr[i];       return ans;}

void arrEachBy(Type* firstArr,Type* secondArr,Type** ansArr,int len){       int i=0;/*iterator num*/       /*assertion*/       assertF(firstArr!=NULL,"in arrEachBy,firstArr is NULL\n");       assertF(secondArr!=NULL,"in arrEachBy,secondArr is NULL\n");            assertF(ansArr!=NULL,"in arrEachBy,ansArr is NULL\n");              for(i=0;i<len;i++)              (*ansArr)[i]=firstArr[i]*secondArr[i];}

/*Min2ByMethod  Algorithm test program*/#include "Global.h"#include "Ulti.h"#include "MyAssert.h"#include "Fix.h"#include <time.h>#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>



char *inFileName="inputData.txt";/*       input data specification       len,       x0,x1,...,xn;       y0,y1,...,yn;*/



char *outFileName="outputData.txt";#define DEBUG 1

void main(int argc,char* argv[]){       FILE *inputFile;/*input file*/       FILE *outputFile;/*output file*/

       double startTime,endTime,tweenTime;/*time callopsed info*/              /*The read in data*/       int len;       Type* xList,* yList;              int i;/*iterator index*/              /*input file open*/       if(argc>1)strcpy(inFileName,argv[1]);       assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");       printf("input file open success\n");              /*outpout file open*/       if(argc>2)strcpy(outFileName,argv[2]);       assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");       printf("output file open success\n");                   fscanf(inputFile,"%d,",&len);       /*Memory apply*/       xList=(Type*)malloc(sizeof(Type)*len);       yList=(Type*)malloc(sizeof(Type)*len);       /*Read info data*/       for(i=0;i<len-1;i++)                     fscanf(inputFile,"%f,",&xList[i]);       fscanf(inputFile,"%f;",&xList[i]);              for(i=0;i<len-1;i++)                     fscanf(inputFile,"%f,",&yList[i]);       fscanf(inputFile,"%f;",&yList[i]);              printf("read in data info,len:%d,\n",len);       showArrListFloat(xList,0,len);       showArrListFloat(yList,0,len);

#if  DEBUG       printf("\n*******start of test program******\n");       printf("now is runnig,please wait...\n");       startTime=(double)clock()/(double)CLOCKS_PER_SEC;       /******************Core program code*************/              min2ByFixMethod1(xList,yList,len,outputFile);       /******************End of Core program**********/       endTime=(double)clock()/(double)CLOCKS_PER_SEC;       tweenTime=endTime-startTime;/*Get the time collapsed*/       /*Time collapsed output*/       printf("the collapsed time in this algorithm implement is:%f\n",tweenTime);       fprintf(outputFile,"the collapsed time in this algorithm implement is:%f\r\n",tweenTime);         printf("\n*******end of test program******\n");#endif       free(xList);       free(yList);       printf("program end successfully,\n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer.\n");       getchar();/*Screen Delay Control*/       return;}

//测试数据17,0.6,1.3,1.64,1.8,2.1,2.3,2.44;7.05,12.2,14.4,15.2,17.4,19.6,20.2;//输出:the fix fun is 2.708288+7.150409x.

//测试数据22,0,1;3,5;the fix fun is 3.000000+2.000000x.

可以.但应该分为上下两个部分才行,毕竟拟合的二次曲线模型是y=ax^2+bx+c.

在拟合曲线前，你可以简单的在草纸上把这些点勾画一下，连接起来，看一下其基本形状，以决定在程序是把其拟合成什么样子，个人感觉如果超过二次，觉得复杂的话，可以把其拟合为三次的曲线方程。以下程序是我写的，且调试通了的，但是只是针对三次方程的曲线拟合，在早改过一个针对通用一些的程序找不到了。
///////////////////////////
#include <stdio.h>
#include <math.h>double data[2][7]={{-1.0,-0.5,0.0,0.5,1.0,1.5,2.0},{-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552}};
double G[4][4]={0},ZG_G[4][5]={0};
double b[4]={0};      //Gx=b
double result[4]={0}; //解向量
double  function(double x);
void main(void)
{
int i;
void cal_G(void);
    void print1(void);
void print2(void);
void cal_matrix(void);

double cal_wu(void);
cal_G();
printf("系数矩阵：\n");
print1();         //打印系数矩阵
printf("增广矩阵：\n");
print2();         //打印增广矩阵
cal_matrix();
printf("高斯消元后的矩阵：\n");
print2();
printf("高斯法得解为：\n");
for(i=0;i<4;i++)
printf("%lf  ",result[i]);
printf("\n");
printf("平方误差为：%lf\n",cal_wu());
}
void cal_G(void)
{
int i,j,k;
double temp1;
for(i=0;i<4;i++)
for(j=0;j<4;j++)
for(k=0;k<7;k++)
{
temp1=i+j;
G[i][j]+=pow(data[0][k],temp1);
}
//~~~~~~~~~~~~~~~~~~
for(i=0;i<4;i++)
for(j=0;j<4;j++)
ZG_G[i][j]=G[i][j];
//""""""""""""""""""
for(i=0;i<4;i++)
for(k=0;k<7;k++)
{
temp1=i;
b[i]+=pow(data[0][k],temp1)*data[1][k];
}
//////////////////////////////
for(i=0;i<4;i++)
ZG_G[i][4]=b[i];
//\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////
printf("没有消元前矩阵右侧系数b[]:\n");
for(k=0;k<4;k++)
printf("%lf ",b[k]);
printf("\n");printf("\n");
}
void cal_matrix(void)
{
int i,j,k;
double temp1;
i=0;k=0;
for(k=0;k<4;k++)    //高斯消元过程
{
for(i=k+1;i<5;i++)
{
temp1=ZG_G[i][k]/ZG_G[k][k];
for(j=0;j<5;j++)
{
ZG_G[i][j]=ZG_G[i][j]-temp1*ZG_G[k][j];
}
}
}
//------------------------------------
for(i=0;i<4;i++)
b[i]=ZG_G[i][4];
printf("\n高斯消去后的矩阵右侧系数：\n");
for(k=0;k<4;k++)
printf("%lf ",b[k]);
printf("\n");
result[3]=b[3]/ZG_G[3][3];    //迭代初值//
for(k=2;k>=0;k--)             //迭代求解//
{
for(j=k+1;j<=3;j++)
{
b[k]=b[k]-ZG_G[k][j]*result[j];
//result[b]=b[j]/ZG_G[][]
}
result[k]=b[k]/ZG_G[k][k];
}
}
double  function(double x)
{
double var;
var=result[0]+result[1]*x+result[2]*x*x+result[3]*x*x*x;
    return var;
}
//-----------------------------
double cal_wu(void)
{
double w_cha=0;
for(int i=0;i<7;i++)
{
w_cha+=(data[1][i]-function(data[0][i]))*(data[1][i]-function(data[0][i]));
printf("%lf  \n",w_cha);
printf("Hello world!");
}
return w_cha;
}
//00000000000000000000000000000
void print1(void)
{
int i,j;
for(i=0;i<4;i++)
{
for(j=0;j<4;j++)
{
printf("%lf  ",G[i][j]);
}
printf("\n");
}
printf("\n");
}
void print2(void)
{
int i,j;
for(i=0;i<4;i++)
{
for(j=0;j<5;j++)
{
printf("%4.6lf  ",ZG_G[i][j]);
}
printf("\n");
}
printf("\n");
}//运行结果
没有消元前矩阵右侧系数b[]:
0.354000         14.177500       14.448250       42.33137系数矩阵：
7.000000  3.500000  8.750000  11.375000
3.500000  8.750000  11.375000  23.187500
8.750000  11.375000  23.187500  39.593750
11.375000  23.187500  39.593750  77.421875增广矩阵：
7.000000  3.500000  8.750000  11.375000  0.354000
3.500000  8.750000  11.375000  23.187500  14.177500
8.750000  11.375000  23.187500  39.593750  14.448250
11.375000  23.187500  39.593750  77.421875  42.331375
高斯消去后的矩阵右侧系数：
0.354000         14.000500       0.005250        6.747000
高斯消元后的矩阵：
7.000000  3.500000  8.750000  11.375000  0.354000
0.000000  7.000000  7.000000  17.500000  14.000500
0.000000  0.000000  5.250000  7.875000  0.005250
0.000000  0.000000  0.000000  3.375000  6.747000高斯法得解为：
0.549119  -0.000040  -2.997667  1.999111
0.000000
Hello world!0.000004
Hello world!0.000007
Hello world!0.000010
Hello world!0.000016
Hello world!0.000021
Hello world!0.000022
Hello world!平方误差为：0.000022
Press any key to continue高斯法解得的四个系数即为三次方程的四个系数。楼主如果需要，可以自己尝试着改一下。

数学专业，计算机专业都开。很多其他工科专业也开。
数值分析英文名：numerical analysis
　　概括：研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，计算数学的主体部分。
　　运用数值分析解决问题的过程：实际问题→数学模型→数值计算方法→程序设计→上机计算求出结果
　　数值分析这门学科有如下特点：
　　1·面向计算机
　　2·有可靠的理论分析
　　3·要有好的计算复杂性
　　4·要有数值实验
　　5.要对算法进行误差分析
　　主要内容：插值法，函数逼近，曲线拟和，数值积分，数值微分，解线性方程组的直接方法，解线性方程组的迭代法，非线性方程求根，常微分方程的数值解法。个人推荐Richard L.Burden著的《Numerical Analysis》(7th)既侧重理论又注重实用，读完有种从头至尾的感觉。

cftxlin 写法不错。思路也好。

恩，cftxlin 写法不错。思路也好。

调试易

最小二乘法怎么做曲线拟合？

解决方案 »