标签：曲线拟合 高斯 double cout new 拟合 原理 Matrix
转：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');
标签：曲线拟合,高斯,double,cout,new,拟合,原理,Matrix 来源： https://www.cnblogs.com/sggggr/p/14244678.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。