声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4157|回复: 3

[共享资源] 混合高斯模型的曲线拟合matlab源代码

[复制链接]
发表于 2006-11-14 18:15 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  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. % run with default values
  17. % if ~nargin
  18. %    % M   = round(rand*5)+1;
  19. %    M=1;
  20. %     sig = rand(1,M)*3;
  21. %     u   = randn(1,M)*8;
  22. %     prob= rand(1,M);
  23. %     [u,sig,t,iter] = fit_mix_gaussian( build_mix_gaussian( u,sig,prob,1000*M ),M );
  24. %     return
  25. % end
  26. %%%%%%%%added by me%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  27.    M=2;
  28.    sig = 3;
  29.    u   = 8;
  30.    prob= 1;
  31. [X,N]=build_mix_gaussian( u,sig,prob,256*M );
  32. figure;
  33. plot(X);

  34. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  35. % initialize and initial guesses
  36. N           = length( X );
  37. Z           = ones(N,M) * 1/M;                  % indicators vector
  38. P           = zeros(N,M);                       % probabilities vector for each sample and each model
  39. t           = ones(1,M) * 1/M;                  % distribution of the gaussian models in the samples
  40. u           = linspace(min(X),max(X),M);        % mean vector
  41. sig2        = ones(1,M) * var(X) / sqrt(M);     % variance vector
  42. C           = 1/sqrt(2*pi);                     % just a constant
  43. Ic          = ones(N,1);                        % - enable a row replication by the * operator
  44. Ir          = ones(1,M);                        % - enable a column replication by the * operator
  45. Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solution
  46. thresh      = 1e-3;
  47. step        = N;
  48. last_step   = inf;
  49. iter        = 0;
  50. min_iter    = 10;

  51. % main convergence loop, assume gaussians are 1D
  52. while ((( abs((step/last_step)-1) > thresh) & (step>(N*eps)) ) | (iter<min_iter) )
  53.    
  54.     % E step
  55.     % ========
  56.     Q   = Z;
  57.     P   = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
  58.     %%%%%%1/N*1X1*M*exp((N*1XM-N*1X1*M)/N*1*1*M)
  59.    
  60.    
  61.     for m = 1:M
  62.         Z(:,m)  = (P(:,m)*t(m))./(P*t(:));
  63.     end
  64.         
  65.     % estimate convergence step size and update iteration number
  66.     prog_text   = sprintf(repmat( '\b',1,(iter>0)*12+ceil(log10(iter+1)) ));
  67.     iter        = iter + 1;
  68.     last_step   = step * (1 + eps) + eps;
  69.     step        = sum(sum(abs(Q-Z)));
  70.     fprintf( '%s%d iterations\n',prog_text,iter );

  71.     % M step
  72.     % ========
  73.     Zm              = sum(Z);               % sum each column
  74.     Zm(find(Zm==0)) = eps;                  % avoid devision by zero
  75.     u               = (X')*Z ./ Zm;
  76.    
  77.     sig2            = sum(((X*Ir - Ic*u).^2).*Z) ./ Zm;
  78.     t               = Zm/N;
  79. end

  80. % plot the fitted distribution
  81. % =============================
  82. sig     = sqrt( sig2 );
  83. plot_mix_gaussian( u,sig,t );
复制代码
回复
分享到:

使用道具 举报

发表于 2013-10-22 21:12 | 显示全部楼层
代码不全啊    build_mix_gaussian函数没有
发表于 2014-2-17 20:57 | 显示全部楼层
请问一下这个可以用于信号分类吗?
发表于 2014-2-25 01:10 | 显示全部楼层
是这个吗?http://llwebprod2.ll.mit.edu/mis ... _Biometrics-GMM.pdf 感觉好老的算法了。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 13:58 , Processed in 0.067486 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表