高斯曲线拟合原理及实现

高斯拟合(Gaussian Fitting)即使用形如：
    
          Gi(x)=Ai*exp((x-Bi)^2/Ci^2)

        的高斯函数对数据点集进行函数逼近的拟合方法。

        其实可以跟多项式拟合类比起来，不同的是多项式拟合是用幂函数系，
        而高斯拟合是用高斯函数系。

        使用高斯函数来进行拟合，优点在于计算积分十分简单快捷。这一点
        在很多领域都有应用，特别是计算化学。著名的化学软件Gaussian98
        就是建立在高斯基函数拟合的数学基础上的。

c#中用mathnet 惊醒矩阵运算  实现方案

 

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');