最小二乘法線性擬合和2次曲線擬合算法

最近由于項目要求，應用了最小二乘法線性擬合和2次曲線擬合算法，現總結如下：

最小二乘法線性擬合應用已有的采樣時間點，再現這些點所描述的線性變化，即求出一個線性方程y=ax+b(這個算法的主要問題也就是如何用給定的數據求線性方程系數a和b)

//最小二乘法線性擬合,線性方程求系數,Xval時間數據，Yval每個時間點上的值數據，n數據的個數，Aval線性方程系數a,Bval線性方程系數b
BOOL DlgDataAnalyse::TwoCurveCompose(double *Xval,double *Yval,long n,double *Aval,double *Bval)
{
 double mX,mY,mXX,mXY;
 mX=mY=mXX=mXY=0;
 for (int i=0;i  {
  mX+=Xval[i];
  mY+=Yval[i];
  mXX+=Xval[i]*Xval[i];
  mXY+=Xval[i]*Yval[i];
 }
 if(mX*mX-mXX*n==0)return FALSE;
 *Aval=(mY*mX-mXY*n)/(mX*mX-mXX*n);
 *Bval=(mY-mX*(*Aval))/n;
 return TRUE;
}
最小二乘法2次曲線擬合應用已有的采樣時間點，再現這些點所描述的2次曲線的變化，即求出一個二次曲線方程y=ax2+bx+c (這個算法的主要問題也就是如何用給定的數據求方程系數abc)

今天使用擬合的最小二乘法，求出了給定的一組坐標系上的點對最接近的直線的。
　　其具體理論如下：
   在科學實驗數據處理中，往往要根據一組給定的實驗數據，求出自變量x與因變量y的函數關系，這是為待定參數，由于觀測數據總有誤差，且待定參數ai的數量比給定數據點的數量少(即n＜m)，因此它不同于插值問題.這類問題不要求通過點，而只要求在給定點上的誤差的平方和最小.當時，即
　　　　　　　　　(5.8.1)
這里是線性無關的函數族，假定在上給出一組數據，以及對應的一組權，這里為權系數，要求使最小，其中
　　　　　　　　　　　(5.8.2)

這就是最小二乘逼近，得到的擬合曲線為y=s(x)，這種方法稱為曲線擬合的最小二乘法.
　　(5.8.2)中實際上是關于的多元函數，求I的最小值就是求多元函數I的極值，由極值必要條件，可得
　　　　　　(5.8.3)
根據內積定義(見第三章)引入相應帶權內積記號
　　　　　　　　　　　　　(5.8.4)
則(5.8.3)可改寫為
　　　　　
這是關于參數的線性方程組，用矩陣表示為
　　　　　　　　　　　　(5.8.5)
(5.8.5)稱為法方程.當線性無關，且在點集上至多只有n個不同零點，則稱在X上滿足Haar條件，此時(5.8.5)的解存在唯一(證明見［3］).記(5.8.5)的解為
　　　　　　　　 
從而得到最小二乘擬合曲線
　　　　　　　　　　(5.8.6)
可以證明對，有
　　　　　
故(5.8.6)得到的即為所求的最小二乘解.它的平方誤差為
　　　　　　　　　　　　　　　　　　(5.8.7)
均方誤差為
　　　　　　　　　　
　　在最小二乘逼近中，若取，則，表示為
　　　　　　　　　　　　　(5.8.8)
此時關于系數的法方程(5.8.5)是病態方程，通常當n≥3時都不直接取作為基。

//最小二乘法二次曲線擬合算法,Xval時間數據，Yval每個時間點上的值數據，M代表幾次曲線(如：2次的話就是3)，N數據的個數，A二次曲線方程的系數(A[2]代表a,A[1]代表b,A[0]代表c)
BOOL DlgDataAnalyse::CalculateCurveParameter(double *Xval,double *Yval,long M,long N,double *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   A[i]=0;
 Z=0;
 B[0]=1;
 D1=N;
 P=0;
 C=0;
 for(i=0;i  {
  P=P+Xval[i]-Z;
  C=C+Yval[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   {
   Q=Xval[i]-Z-P;
   D2=D2+Q*Q;
   C=Yval[i]*Q+C;
   G=(Xval[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  {
  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   {
   Q=S[j];
   for(k=j-1;k>=0;k--)
    Q=Q*(Xval[i]-Z)+S[k];
   D2=D2+Q*Q;
   C=Yval[i]*Q+C;
   G=(Xval[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;
}