混合高斯模型的曲线拟合matlab源代码

suffer

1. function [u,sig,t,iter] = fit_mix_gaussian( X,M )
   
2. clear all
   
3. %
   
4. % fit_mix_gaussian - fit parameters for a mixed-gaussian distribution using EM algorithm
   
5. %
   
6. % format:   [u,sig,t,iter] = fit_mix_gaussian( X,M )
   
7. %
   
8. % input:    X   - input samples, N*1 vector
   
9. %           M   - number of gaussians which are assumed to compose the distribution
   
10. %
    
11. % output:   u   - fitted mean for each gaussian
    
12. %           sig - fitted standard deviation for each gaussian
    
13. %           t   - probability of each gaussian in the complete distribution
    
14. %           iter- number of iterations done by the function
    
15. %
    
16. 
    
17. 
    
18. % run with default values
    
19. % if ~nargin
    
20. %    % M   = round(rand*5)+1;
    
21. %    M=1;
    
22. %     sig = rand(1,M)*3;
    
23. %     u   = randn(1,M)*8;
    
24. %     prob= rand(1,M);
    
25. %     [u,sig,t,iter] = fit_mix_gaussian( build_mix_gaussian( u,sig,prob,1000*M ),M );
    
26. %     return
    
27. % end
    
28. %%%%%%%%added by me%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
29.    M=2;
    
30.    sig = 3;
    
31.    u   = 8;
    
32.    prob= 1;
    
33. [X,N]=build_mix_gaussian( u,sig,prob,256*M );
    
34. figure;
    
35. plot(X);
    
36. 
    
37. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
38. 
    
39. % initialize and initial guesses
    
40. N           = length( X );
    
41. Z           = ones(N,M) * 1/M;                  % indicators vector
    
42. P           = zeros(N,M);                       % probabilities vector for each sample and each model
    
43. t           = ones(1,M) * 1/M;                  % distribution of the gaussian models in the samples
    
44. u           = linspace(min(X),max(X),M);        % mean vector
    
45. sig2        = ones(1,M) * var(X) / sqrt(M);     % variance vector
    
46. C           = 1/sqrt(2*pi);                     % just a constant
    
47. Ic          = ones(N,1);                        % - enable a row replication by the * operator
    
48. Ir          = ones(1,M);                        % - enable a column replication by the * operator
    
49. Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solution
    
50. thresh      = 1e-3;
    
51. step        = N;
    
52. last_step   = inf;
    
53. iter        = 0;
    
54. min_iter    = 10;
    
55. 
    
56. % main convergence loop, assume gaussians are 1D
    
57. while ((( abs((step/last_step)-1) > thresh) & (step>(N*eps)) ) | (iter<min_iter) )
    
58.    
    
59.     % E step
    
60.     % ========
    
61.     Q   = Z;
    
62.     P   = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
    
63.     %%%%%%1/N*1X1*M*exp((N*1XM-N*1X1*M)/N*1*1*M)
    
64.    
    
65.    
    
66.     for m = 1:M
    
67.         Z(:,m)  = (P(:,m)*t(m))./(P*t(:));
    
68.     end
    
69.         
    
70.     % estimate convergence step size and update iteration number
    
71.     prog_text   = sprintf(repmat( '\b',1,(iter>0)*12+ceil(log10(iter+1)) ));
    
72.     iter        = iter + 1;
    
73.     last_step   = step * (1 + eps) + eps;
    
74.     step        = sum(sum(abs(Q-Z)));
    
75.     fprintf( '%s%d iterations\n',prog_text,iter );
    
76. 
    
77.     % M step
    
78.     % ========
    
79.     Zm              = sum(Z);               % sum each column
    
80.     Zm(find(Zm==0)) = eps;                  % avoid devision by zero
    
81.     u               = (X')*Z ./ Zm;
    
82.    
    
83.     sig2            = sum(((X*Ir - Ic*u).^2).*Z) ./ Zm;
    
84.     t               = Zm/N;
    
85. end
    
86. 
    
87. % plot the fitted distribution
    
88. % =============================
    
89. sig     = sqrt( sig2 );
    
90. plot_mix_gaussian( u,sig,t );

复制代码

aweyiun

代码不全啊    build_mix_gaussian函数没有

褚青青

请问一下这个可以用于信号分类吗？

mayaview

是这个吗？http://llwebprod2.ll.mit.edu/mis ... _Biometrics-GMM.pdf 感觉好老的算法了。