高斯(Gaussian)拟合的实现

转载于:https://blog.csdn.net/c914620529/article/details/50393238

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

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

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

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




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

 


 
  1. double[,] a = new double[fitDatas.Count, 3];
  2. double[] b = new double[fitDatas.Count];
  3. double[] X = new double[ 3] { 0, 0, 0 };
  4. for ( int i = 0; i < fitDatas.Count; i++)
  5. {
  6. b[i] = Math.Log(fitDatas[i].Intensity);
  7. a[i, 0] = 1;
  8. a[i, 1] = fitDatas[i].WaveLength;
  9. a[i, 2] = a[i, 1] * a[i, 1];
  10. }
  11. // Matrix.Equation(datas.Count, 3, a, b, X);
  12. MathNet.Numerics.LinearAlgebra.Matrix matrixA = new MathNet.Numerics.LinearAlgebra.Matrix(a);
  13. MathNet.Numerics.LinearAlgebra.Matrix matrixB = new MathNet.Numerics.LinearAlgebra.Matrix(b, b.Length);
  14. MathNet.Numerics.LinearAlgebra.Matrix matrixC = matrixA.Solve(matrixB);
  15. X = matrixC.GetColumnVector( 0);
  16. double S = -1 / X[ 2];
  17. double xMax = X[ 1] * S / 2.0;
  18. double yMax = Math.Exp(X[ 0] + xMax * xMax / S);



运用c++实现方案


 
  1. #include<iostream.h>
  2. #include<math.h>
  3. #include<stdlib.h>
  4. #include <windows.h>
  5. double f(int n,double x){ //f(n,x)用来返回x的n次方
  6. double y= 1.0;
  7. if(n== 0) return 1.0;
  8. else{
  9. for( int i= 0;i<n;i++)y*=x;
  10. return y;
  11. }
  12. }
  13. int xianxingfangchengzu(double **a,int n,double *b,double *p,double dt)//用高斯列主元法来求解法方程组
  14. {
  15. int i,j,k,l;
  16. double c,t;
  17. for(k= 1;k<=n;k++)
  18. {
  19. c= 0.0;
  20. for(i=k;i<=n;i++)
  21. if( fabs(a[i -1][k -1])> fabs(c))
  22. {
  23. c=a[i -1][k -1];
  24. l=i;
  25. } if( fabs(c)<=dt)
  26. return( 0);
  27. if(l!=k)
  28. {
  29. for(j=k;j<=n;j++)
  30. {
  31. t=a[k -1][j -1];
  32. a[k -1][j -1]=a[l -1][j -1];
  33. a[l -1][j -1]=t;
  34. }
  35. t=b[k -1];
  36. b[k -1]=b[l -1];
  37. b[l -1]=t;
  38. }
  39. c= 1/c;
  40. for(j=k+ 1;j<=n;j++)
  41. {
  42. a[k -1][j -1]=a[k -1][j -1]*c;
  43. for(i=k+ 1;i<=n;i++)
  44. a[i -1][j -1]-=a[i -1][k -1]*a[k -1][j -1];
  45. }
  46. b[k -1]*=c;
  47. for(i=k+ 1;i<=n;i++)
  48. b[i -1]-=b[k -1]*a[i -1][k -1];
  49. }
  50. for(i=n;i>= 1;i--)
  51. for(j=i+ 1;j<=n;j++)
  52. b[i -1]-=b[j -1]*a[i -1][j -1];
  53. cout.precision( 12);
  54. for(i= 0;i<n;i++)p[i]=b[i];
  55. }
  56. double** create(int a,int b)//动态生成数组
  57. {
  58. double **P= new double *[a];
  59. for( int i= 0;i<b;i++)
  60. P[i]= new double[b];
  61. return P;
  62. }
  63. void zuixiaoerchengnihe(double x[],double y[],int n,double a[],int m)
  64. {
  65. int i,j,k,l;
  66. double **A,*B;
  67. A=create(m,m);
  68. B= new double[m];
  69. for(i= 0;i<m;i++)
  70. for(j= 0;j<m;j++)A[i][j]= 0.0;
  71. for(k= 0;k<m;k++)
  72. for(l= 0;l<m;l++)
  73. for(j= 0;j<n;j++)A[k][l]+=f(k,x[j])*f(l,x[j]); //计算法方程组系数矩阵A[k][l]
  74. cout<< "法方程组的系数矩阵为:"<< endl;
  75. for(i= 0;i<m;i++)
  76. for(j= 0,k= 1;j<m;j++,k++){
  77. cout<<A[i][j]<< '\t';
  78. if(k&&k%m== 0) cout<< endl;
  79. }
  80. for(i= 0;i<m;i++)B[i]= 0.0;
  81. for(i= 0;i<m;i++)
  82. for(j= 0;j<n;j++)B[i]+=y[j]*f(i,x[j]);
  83. for(i= 0;i<m;i++) cout<< "B["<<i<< "]="<<B[i]<< endl; //记录B[n]
  84. xianxingfangchengzu(A,m,B,a, 1e-6);
  85. delete[]A;
  86. delete B;
  87. }
  88. double pingfangwucha(double x[],double y[],int n,double a[],int m)//计算最小二乘解的平方误差
  89. {
  90. double deta,q= 0.0,r= 0.0;
  91. int i,j;
  92. double *B;
  93. B= new double[m];
  94. for(i= 0;i<m;i++)B[i]= 0.0;
  95. for(i= 0;i<m;i++)
  96. for(j= 0;j<n;j++)B[i]+=y[j]*f(i,x[j]);
  97. for(i= 0;i<n;i++)q+=y[i]*y[i];
  98. for(j= 0;j<m;j++)r+=a[j]*B[j];
  99. deta= fabs(q-r);
  100. return deta;
  101. delete B;
  102. }
  103. void main(void){
  104. int i,n,m;
  105. double *x,*y,*a;
  106. char ch= 'y';
  107. do{
  108. system( "cls");
  109. cout<< "请输入所给拟合数据点的个数n=";
  110. cin>>n;
  111. cout<< "请输入所要拟合多项式的项数m=";
  112. cin>>m;
  113. while(n<=m){
  114. cout<< "你所输入的数据点无法确定拟合项数,请重新输入"<< endl;
  115. Sleep( 1000);
  116. system( "cls");
  117. cout<< "请输入所给拟合数据点的个数n=";
  118. cin>>n;
  119. cout<< "请输入所要拟合多项式的项数m=";
  120. cin>>m;
  121. }
  122. x= new double[n]; //存放数据点x
  123. y= new double[n]; //存放数据点y
  124. a= new double[m]; //存放拟合多项式的系数
  125. cout<< "请输入所给定的"<<n<< "个数据x"<< endl;
  126. for(i= 0;i<n;i++)
  127. {
  128. cout<< "x["<<i+ 1<< "]=";
  129. cin>>x[i];
  130. }
  131. cout<< "请输入所给定的"<<n<< "个数据y"<< endl;
  132. for(i= 0;i<n;i++)
  133. {
  134. cout<< "y["<<i+ 1<< "]=";
  135. cin>>y[i];
  136. }
  137. zuixiaoerchengnihe(x,y,n,a,m+ 1);
  138. cout<< endl;
  139. cout<< "拟合多项式的系数为:"<< endl;
  140. for(i= 0;i<=m;i++) cout<< "a["<<i<< "]="<<a[i]<< '\t';
  141. cout<< endl;
  142. cout<< "平方误差为:"<<pingfangwucha(x,y,n,a,m+ 1)<< endl;
  143. delete x; delete y;
  144. cout<< "按y继续,按其他字符退出"<< endl;
  145. cin>>ch;
  146. } while(ch== 'y'||ch== 'Y');




  • 2
    点赞
  • 14
    收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页
评论 5
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值