www.pudn.com > voicebox2.zip > GAUSSMIX.M

function [m,v,w,lp,f]=gaussmix(x,c,l,m0,v0,w0) 
%GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,lp]=(x,xv,l,m0,v0,w0) 
% 
% n data values, k mixtures, p parameters, l loops 
% Inputs: 
%     X(n,p)   Input data vectors, one per row. 
%     c(1,p)  Minimum variance (can be a scalar if all components are identical) 
%     L        The integer portion of l gives a maximum loop count. The fractional portion gives 
%              an optional stopping threshold. Iteration will cease if the average increase in 
%              log likelihood density per data point is less than this value. Thus l=10.001 will 
%              stop after 10 iterations or when the average increase in log likelihood falls below 
%              0.001. 
%     M0(k,p)  Initial mixture means, one row per mixture. Alternatively set M0 to the number of 
%              mixtures that are wanted and omit V0 and W0; in this case the kmeans algorithm 
%              is used for initialisation. 
%     V0(k,p)  Initial mixture variances, one row per mixture. 
%     W0(k,1)  Initial mixture weights, one per mixture. The weights should sum to unity. 
% 
% Outputs: 
%     M(k,p)   Mixture means, one row per mixture. 
%     V(k,p)   Mixture variances, one row per mixture. 
%     W(k,1)   Mixture weights, one per mixture. The weights will sum to unity. 
%     LP       Total log probability of the input data. 
%     F        Fisher's F ratio measures how well the data divides into classes. 
 
%      Copyright (C) Mike Brookes 2000 
% 
%      Last modified Mon Apr 17 17:37:57 2000 
% 
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at 
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   This program is free software; you can redistribute it and/or modify 
%   it under the terms of the GNU General Public License as published by 
%   the Free Software Foundation; either version 2 of the License, or 
%   (at your option) any later version. 
% 
%   This program is distributed in the hope that it will be useful, 
%   but WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
%   GNU General Public License for more details. 
% 
%   You can obtain a copy of the GNU General Public License from 
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to 
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
[n,p]=size(x); 
x2=x.^2; 
if length(c)>1 
   c=c(ones(k,1),:); 
end 
if nargin<6 
   k=m0; 
   [m,e,j]=kmeans(x,k); 
   v=zeros(k,p); 
   w=zeros(k,1); 
   for i=1:k 
      s=i==j; 
      v(i,:)=mean(x2(s,:))-m(i,:).^2; 
      w(i)=sum(s)/n; 
   end 
else 
   k=size(m0,1); 
   m=m0; 
   v=v0; 
   w=w0; 
end 
v=max(v,c); 
im=repmat(1:k,1,n); im=im(:); 
ix=repmat(1:n,k,1); ix=ix(:); 
th=(l-floor(l))*n; 
lp=0; 
wk=ones(k,1); 
wp=ones(1,p); 
wn=ones(1,n); 
 
% EM loop 
 
for j=1:l 
   vi=v.^(-1); 
   vm=sqrt(prod(vi,2)).*w; 
   vi=-0.5*vi; 
   % Next three lines are equivalent to the one below but with less danger of underflow 
   %          px=reshape(exp(sum((x(ix,:)-m(im,:)).^2.*vi(im,:),2)).*vm(im),k,n); 
   px=reshape(sum((x(ix,:)-m(im,:)).^2.*vi(im,:),2),k,n); 
   mx=max(px,[],1); 
   px=exp(px-mx(wk,:)).*vm(:,wn); 
   ps=sum(px,1); 
   px=px./ps(wk,:); 
   pk=sum(px,2); 
   w=pk/n; 
   m=px*x./pk(:,wp); 
   v=px*x2./pk(:,wp); 
   v=max(v-m.^2,c); 
   lp0=lp; 
   lp=sum(log(ps))+sum(mx); 
   if j>1 & lp<=lp0+th break; end 
end 
if nargout > 3 
   lp=lp-0.5*p*log(2*pi); 
   if nargout > 4 
      mm=sum(m,1)/k; 
      f=prod(sum(m.^2,1)/k-mm.^2)/prod(sum(v,1)/k); 
   end 
end