Home > @ahrbf > private > gmem.m

gmem

PURPOSE ^

GMEM 1-d Gaussian-Mixture EM

SYNOPSIS ^

function [beta,g,P2,P0]=gmem(x,f,beta0,varargin);

DESCRIPTION ^

 GMEM 1-d Gaussian-Mixture EM

 [beta,g]=GMEM(x,f,beta0)
 min_{beta} {|g(x)-f(x;beta)|}
 where beta=[A,m,r]; 
 g=sum_i{1/sqrt(2*pi*r[i])*A[i]*exp(-(x-m[i])^2/r[i])};

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [beta,g,P2,P0]=gmem(x,f,beta0,varargin);
0002 % GMEM 1-d Gaussian-Mixture EM
0003 %
0004 % [beta,g]=GMEM(x,f,beta0)
0005 % min_{beta} {|g(x)-f(x;beta)|}
0006 % where beta=[A,m,r];
0007 % g=sum_i{1/sqrt(2*pi*r[i])*A[i]*exp(-(x-m[i])^2/r[i])};
0008 %
0009 
0010 % 10/00
0011 % alfnie@bu.edu
0012 %
0013 
0014 TOL=.005; % percentage change
0015 options=strvcat(varargin{:});
0016 sf=sum(f);
0017 f=f(:)/sf;
0018 x=x(:);
0019 A=beta0(:,1);
0020 m=beta0(:,2);
0021 r=beta0(:,3);
0022 Ni=length(A);
0023 Nx=length(x);
0024 dx=(x(end)-x(1))/(Nx-1);
0025 minr=dx.^2;
0026 
0027 beta=beta0; err=inf; 
0028 while err>TOL,
0029     % E-Step
0030     % From model -> P(x|i)
0031     xm=(x(:,ones(1,Ni)).'-m(:,ones(1,Nx))).^2;
0032     P0=exp(-xm./(2*r(:,ones(1,Nx))))./sqrt(2*pi*r(:,ones(1,Nx)))*dx;
0033     P1=A(:,ones(1,Nx)).*P0;
0034     
0035     % From P(x|i) to P(i|x)
0036     g=sum(P1,1); 
0037     P2=P1./(eps+g(ones(1,Ni),:));
0038     
0039     % M-Step
0040     % P(i)
0041     A=eps+P2*f;
0042     % <x|i>
0043     m=(P2*(f.*x))./A;
0044     % <(x-m)^2|i>
0045     if strmatch('comvar',options),     % common variance
0046         r=sum((P2.*xm)*f)*ones(Ni,1);
0047     else                                         % default
0048         r=((P2.*xm)*f)./A;
0049     end
0050     r=max(minr,r);
0051     
0052     betanew=[A,m,r];
0053     err=abs(betanew-beta)./(eps+beta); err=max(err(:));
0054     beta=betanew;
0055 end
0056 %plot(x,f,'b',x,g,'r'); hold on; stem(m,A*dx./sqrt(2*pi*r)); hold off; drawnow
0057 beta(:,1)=beta(:,1)*sf;
0058 g=g*sf;
0059

Generated on Tue 27-Mar-2007 12:06:24 by m2html © 2003