豆粑粑 fit Gaussian matlab

meatball1982

matlab fit Gaussian distribution
init parameter is max, std and  mean of y.
the both of the lsqcurvefit and nlinfit are applied to iterative



fun_fit_gau.m

1. function [x_new,y_new,b]=fun_fit_gau(x,y)
   
2. 
   
3. b(1)=max(y);
   
4. b(2)= std(y);
   
5. b(3)=mean(y);
   
6. 
   
7. fx=@(b,x)b(1)*exp(-b(2)*(x-b(3)).^2);
   
8. for l=1:3
   
9. b=lsqcurvefit(fx,b,x,y);
   
10. b=nlinfit(x,y,fx,b);
    
11. end
    
12. 
    
13. 
    
14. x_new=interp1([0:length(x)-1],x,[0:0.01:length(x)-1]);
    
15. y_new = b(1)*exp(-b(2)*(x_new-b(3)).^2);
    

复制代码

main.m

1. clear all
   
2. clc
   
3. clf
   
4. 
   
5. 
   
6. 
   
7. %% outline
   
8. %
   
9. 
   
10. %% main
    
11. dat1 = load('../data/deltE_QM1.dat');
    
12. 
    
13. hi_bi = linspace(-20,5,26);
    
14. 
    
15. [hi_va1]=hist(dat1(:,2),hi_bi);
    
16. hi_va1=hi_va1./sum(hi_va1);
    
17. 
    
18. [hi_bi1_new,hi_va1_new,b1]=fun_fit_gau(hi_bi,hi_va1);
    
19. 
    
20. pl_max = 1.2*max(hi_va1);
    
21. 
    
22. subplot(2,2,1)
    
23. hold on
    
24. plot(dat1(:,1),dat1(:,2),'m-','linewidth',1.5)
    
25. grid on
    
26. box on
    
27. axis([0 1000 -20 5])
    
28. title('data')
    
29. 
    
30. subplot(2,2,2)
    
31. hold on
    
32. barh(hi_bi,hi_va1,'m')
    
33. axis([0 pl_max -20 5])
    
34. grid on
    
35. box on
    
36. title('hist')
    
37. 
    
38. subplot(2,2,3)
    
39. hold on
    
40. plot(hi_bi,hi_va1,'mo')
    
41. plot(hi_bi1_new,hi_va1_new,'r-','linewidth',2)
    
42. 
    
43. box on
    
44. grid on
    
45. axis([-20 5 0 pl_max])
    
46. title('Gau fit')
    
47. 
    
48. h=gcf;
    
49. fi_na='../file_imgs/fig_gau_1dat';
    
50. fun_work_li_035_myfig_out(h,fi_na,3);
    

复制代码

meatball1982

一个较新的版本

1. function [Gau_fit,Gau_par,x_fit,y_fit,fx]=fun_mm_fit_1Gau(x,y)
   
2. %%
   
3. %
   
4. % [Gau_fit,Gau_par,x_fit,y_fit]=fun_mm_fit_1Gau(x,y)
   
5. %
   
6. % %--------------------------------------------------------------
   
7. % Inputs:
   
8. % x,y
   
9. % Outputs:
   
10. % Gau_fit : the Gaussian fit value of x
    
11. % Gau_par : the parameters in
    
12. %           b(1)*exp(-b(2)*(x-b(3)).^2)+b(4)
    
13. % x_fit   : more x points for a smoother line
    
14. % y_fit   : more y points for a smoother line
    
15. % fx      : the function to calculate the points users define
    
16. %         % x_new = [ 1 3 17 30 40];
    
17. %         % y_new = fx(Gau_par,x_new);
    
18. %
    
19. % %--------------------------------------------------------------
    
20. % example
    
21. % np = 1000;
    
22. % p = 10* randn(1,np)+3 +10;
    
23. % hi_bi = linspace(-40,60,50);
    
24. % hi_va = hist(p,hi_bi);
    
25. %
    
26. % x = hi_bi;
    
27. % y = hi_va;
    
28. % [Gau_fit,Gau_par,x_fit,y_fit,fx]=fun_mm_fit_1Gau(x,y);
    
29. %
    
30. %
    
31. % x_new = [ 1.2 3.3 17.8 30.001 40.08];
    
32. % y_new = fx(Gau_par,x_new);
    
33. %
    
34. % hold on
    
35. % stairs(hi_bi,hi_va,'b-')
    
36. % plot(x,Gau_fit,'gs')
    
37. % plot(x_fit,y_fit,'k-','linewidth',2)
    
38. % plot(x_new,y_new,'ro','markersize',10)
    
39. %
    
40. % %--------------------------------------------------------------
    
41. % type by : mm
    
42. % contact me : meatball1982@163.com
    
43. % typed on : 01-Jun-2021 14:11:21
    
44. 
    
45. %% main
    
46. % define the original function
    
47. fx=@(b,x)b(1)*exp(-b(2)*(x-b(3)).^2)+b(4);
    
48. 
    
49. % point number
    
50. n_p = length(x);
    
51. 
    
52. % smooth version of y, for find the initial b(3) with one peak
    
53. y_sm = smooth(y,ceil(n_p/10)+1);
    
54. [pks_va,pks_ind]=findpeaks(y_sm);
    
55. 
    
56. % for lsq fit, there is value < 0.
    
57. y(y==0)=-0.001;
    
58. b  = ones(1,4);
    
59. b(3)=x(pks_ind(1));
    
60. b(1)=max(y);
    
61. 
    
62. for i_ite=1:10
    
63.     % Solve nonlinear curve-fitting problems in least-squares sense.
    
64.     % just a guess.
    
65.     b=lsqcurvefit(fx,b,x,y);
    
66.     % Nonlinear regression,
    
67.     b=nlinfit(x,y,fx,b);
    
68. end
    
69. 
    
70. % more points and smoother
    
71. x_fit = linspace(min(x),max(x),n_p*5);
    
72. y_fit = fx(b,x_fit);
    
73. y_fit(y_fit<0)=0;
    
74. 
    
75. % Fit points
    
76. Gau_fit = fx(b,x);
    
77. Gau_fit(Gau_fit<0)=0;
    
78. 
    
79. % parameters
    
80. Gau_par = b;
    
81. 
    
82. 
    
83. %% logs
    
84. % mod : 01-Jun-2021 13:24:37
    
85. %
    

复制代码