高斯曲線擬合原理及實現

轉：https://blog.csdn.net/c914620529/article/details/50393238/

高斯擬合(Gaussian Fitting)即使用形如：

    
          Gi(x)=Ai*exp((x-Bi)^2/Ci^2)

        的高斯函數對數據點集進行函數逼近的擬合方法。

        其實可以跟多項式擬合類比起來，不同的是多項式擬合是用冪函數系，
        而高斯擬合是用高斯函數系。

        使用高斯函數來進行擬合，優點在於計算積分十分簡單快捷。這一點
        在很多領域都有應用，特別是計算化學。著名的化學軟件Gaussian98
        就是建立在高斯基函數擬合的數學基礎上的。

 

 c#中用math.net 進行矩陣運算  實現方案

double[,] a = new double[fitDatas.Count, 3];
                        double[] b = new double[fitDatas.Count];
                        double[] X = new double[3] { 0, 0, 0 };
                        for (int i = 0; i < fitDatas.Count; i++)
                        {
                            b[i] = Math.Log(fitDatas[i].Intensity);
                            a[i, 0] = 1;
                            a[i, 1] = fitDatas[i].WaveLength;
                            a[i, 2] = a[i, 1] * a[i, 1];
                        }
                        // Matrix.Equation(datas.Count, 3, a, b, X);
                        MathNet.Numerics.LinearAlgebra.Matrix matrixA = new MathNet.Numerics.LinearAlgebra.Matrix(a);
                        MathNet.Numerics.LinearAlgebra.Matrix matrixB = new MathNet.Numerics.LinearAlgebra.Matrix(b, b.Length);
                        MathNet.Numerics.LinearAlgebra.Matrix matrixC = matrixA.Solve(matrixB);
                        X = matrixC.GetColumnVector(0);
                        double S = -1 / X[2];
                        double xMax = X[1] * S / 2.0;
                        double yMax = Math.Exp(X[0] + xMax * xMax / S);

運用c++實現方案

#include<iostream.h>
 
#include<math.h>
 
#include<stdlib.h>
 
#include <windows.h>
 
double f(int n,double x){             //f(n,x)用來返回x的n次方
 
        double y=1.0;
 
        if(n==0)return 1.0;
 
               else{
 
               for(int i=0;i<n;i++)y*=x;
 
               return y;
 
        }
 
}
 
int xianxingfangchengzu(double **a,int n,double *b,double *p,double dt)//用高斯列主元法來求解法方程組
 
 {
 
   int i,j,k,l;
 
   double c,t;
 
   for(k=1;k<=n;k++)
 
   {
 
   c=0.0;
 
   for(i=k;i<=n;i++)
 
           if(fabs(a[i-1][k-1])>fabs(c))
 
           {
 
           c=a[i-1][k-1];
 
           l=i;
 
           }if(fabs(c)<=dt)
 
                  return(0);
 
           if(l!=k)
 
           {
 
                  for(j=k;j<=n;j++)
 
                  {
 
                          t=a[k-1][j-1];
 
                          a[k-1][j-1]=a[l-1][j-1];
 
                          a[l-1][j-1]=t;
 
                  }
 
                  t=b[k-1];
 
                  b[k-1]=b[l-1];
 
                  b[l-1]=t;
 
           }
 
                  c=1/c;
 
                  for(j=k+1;j<=n;j++)
 
                  {
 
                          a[k-1][j-1]=a[k-1][j-1]*c;
 
                          for(i=k+1;i<=n;i++)
 
                                  a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
 
                  }
 
                  b[k-1]*=c;
 
                   for(i=k+1;i<=n;i++)
 
                               b[i-1]-=b[k-1]*a[i-1][k-1];
 
           }
 
           for(i=n;i>=1;i--)
 
                  for(j=i+1;j<=n;j++)
 
                          b[i-1]-=b[j-1]*a[i-1][j-1];
 
                  
 
   cout.precision(12);
 
   for(i=0;i<n;i++)p[i]=b[i];
 
}
 
double** create(int a,int b)//動態生成數組
 
{
 
        double **P=new double *[a];
 
        for(int i=0;i<b;i++)
 
               P[i]=new double[b];
 
        return P;
 
}
 
 
 
void zuixiaoerchengnihe(double x[],double y[],int n,double a[],int m)
 
{
 
        int i,j,k,l;
 
        double **A,*B;
 
        A=create(m,m);
 
        B=new double[m];
 
        for(i=0;i<m;i++)
 
               for(j=0;j<m;j++)A[i][j]=0.0;
 
               for(k=0;k<m;k++)
 
                       for(l=0;l<m;l++)
 
                               for(j=0;j<n;j++)A[k][l]+=f(k,x[j])*f(l,x[j]);//計算法方程組系數矩陣A[k][l]
 
        cout<<"法方程組的系數矩陣為:"<<endl;
 
        for(i=0;i<m;i++)
 
               for(j=0,k=1;j<m;j++,k++){
 
                       cout<<A[i][j]<<'\t';
 
                       if(k&&k%m==0)cout<<endl;
 
               }
 
    for(i=0;i<m;i++)B[i]=0.0;
 
        for(i=0;i<m;i++)
 
        for(j=0;j<n;j++)B[i]+=y[j]*f(i,x[j]);
 
        for(i=0;i<m;i++)cout<<"B["<<i<<"]="<<B[i]<<endl;//記錄B[n]
 
        xianxingfangchengzu(A,m,B,a,1e-6);
 
        delete[]A;
 
        delete B;
 
}
 
double pingfangwucha(double x[],double y[],int n,double a[],int m)//計算最小二乘解的平方誤差
 
{
 
        double deta,q=0.0,r=0.0;
 
        int i,j;
 
        double *B;
 
        B=new double[m];
 
        for(i=0;i<m;i++)B[i]=0.0;
 
        for(i=0;i<m;i++)
 
        for(j=0;j<n;j++)B[i]+=y[j]*f(i,x[j]);
 
        for(i=0;i<n;i++)q+=y[i]*y[i];
 
        for(j=0;j<m;j++)r+=a[j]*B[j];
 
        deta=fabs(q-r);
 
        return deta;
 
        delete B;
 
}       
 
void main(void){
 
        int i,n,m;
 
        double *x,*y,*a;
 
        char ch='y';
 
        do{
 
        system("cls");
 
        cout<<"請輸入所給擬合數據點的個數n=";
 
        cin>>n;
 
        cout<<"請輸入所要擬合多項式的項數m=";
 
        cin>>m;
 
        while(n<=m){
 
           cout<<"你所輸入的數據點無法確定擬合項數，請重新輸入"<<endl;
 
           Sleep(1000);
 
           system("cls");
 
           cout<<"請輸入所給擬合數據點的個數n=";
 
           cin>>n;
 
           cout<<"請輸入所要擬合多項式的項數m=";
 
           cin>>m;
 
           }
 
        x=new double[n];                             //存放數據點x
 
        y=new double[n];                             //存放數據點y
 
        a=new double[m];                             //存放擬合多項式的系數
 
        cout<<"請輸入所給定的"<<n<<"個數據x"<<endl;
 
        for(i=0;i<n;i++)
 
        {
 
               cout<<"x["<<i+1<<"]=";
 
               cin>>x[i];
 
        }
 
        cout<<"請輸入所給定的"<<n<<"個數據y"<<endl;
 
        for(i=0;i<n;i++)
 
        {
 
               cout<<"y["<<i+1<<"]=";
 
               cin>>y[i];
 
        }
 
        zuixiaoerchengnihe(x,y,n,a,m+1);
 
    cout<<endl;
 
        cout<<"擬合多項式的系數為:"<<endl;
 
        for(i=0;i<=m;i++)cout<<"a["<<i<<"]="<<a[i]<<'\t';
 
        cout<<endl;
 
        cout<<"平方誤差為："<<pingfangwucha(x,y,n,a,m+1)<<endl;
 
        delete x;   delete y;
 
        cout<<"按y繼續，按其他字符退出"<<endl;
 
        cin>>ch;
 
    }while(ch=='y'||ch=='Y');

 參考資料：https://wenku.baidu.com/view/297f4f5da100a6c30c22590102020740be1ecd09.html