www.pudn.com > box.rar > gmmlpdf.m

function l=gmmlpdf(x,m,v,w) 
%GMMPDF calculated the pdf of a mixture of gaussians p=(x,m,v,w) 
% 
% Inputs: n data values, k mixtures, p parameters 
% 
%     X(n,p)   Input data vectors, one per row. 
%     M(k,p)   mixture means, one row per mixture. 
%     V(k,p)   mixture variances, one row per mixture (singlton dimensions will be replicated as required) 
%              or else V(p,p,k) for full mixture covariance matrixes            
%     W(k,1)   mixture weights, one per mixture. The weights will be normalized by their sum. [default: all equal] 
% 
% Outputs: (Note that M, V and W are omitted if L==0) 
% 
%     L(n,1)   log PDF values 
 
%  Bugs/Suggestions 
%     (1) Sort out full covariance maatrices 
%     (2) Improve plotting 
%     (3) Add an extra arument for plotting control 
 
%      Copyright (C) Mike Brookes 2000-2006 
%      Version: $Id: gmmlpdf.m,v 1.2 2007/05/04 07:01:38 dmb Exp $ 
% 
%   VOICEBOX is a MATLAB toolbox for speech processing. 
%   Home page: 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 
%   http://www.gnu.org/copyleft/gpl.html or by writing to 
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[n,p]=size(x); 
l=[];           % in case n=0 
if nargout>0 || n>0 
    x2=x.^2;            % need x^2 for variance calculation 
    k=size(m,1);        % number of mixtures 
    if size(m,2)~=p 
        error('x and m must have the same number of columns'); 
    end 
    if nargin<4 
        if nargin<3 
            v=1; 
        end 
        w=ones(k,1); 
    end 
    w=w/sum(w);         % normalize the weights 
    sv=size(v); 
    if length(sv)>2 || k==1 && p>1 && sv(1)==p     % full covariance matrices 
        error('full covariance matrices not yet implemented'); 
    else                            % diagonal (or constant) covariance matrices 
        if sv(1)==1 
            v=v(ones(k,1),:); 
        end 
        if sv(2)==1 
            v=v(:,ones(1,p)); 
        end 
         
        % If data size is large then do calculations in chunks 
         
        memsize=voicebox('memsize');  
        nb=min(n,max(1,floor(memsize/(8*p*k))));    % chunk size for testing data points 
        nl=ceil(n/nb);                  % number of chunks 
         
        im=repmat(1:k,1,nb); im=im(:); 
         
        lpx=zeros(n,1);             % log probability of each data point 
        wk=ones(k,1); 
        wnb=ones(1,nb); 
        vi=v.^(-1);                 % calculate quantities that depend on the variances 
        vm=sqrt(prod(vi,2)).*w; 
        vi=-0.5*vi; 
         
        % first do partial chunk 
         
        jx=n-(nl-1)*nb;                % size of first chunk 
        ii=1:jx; 
        kk=repmat(ii,k,1); 
        km=repmat(1:k,1,jx); 
        py=reshape(sum((x(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx); 
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp() 
        px=exp(py-mx(wk,:)).*vm(:,ones(1,jx));  % find normalized probability of each mixture for each datapoint 
        lpx(ii)=log(sum(px,1))+mx; 
        ix=jx+1; 
         
        for il=2:nl 
            jx=jx+nb;        % increment upper limit 
            ii=ix:jx; 
            kk=repmat(ii,k,1); 
            py=reshape(sum((x(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb); 
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp() 
            px=exp(py-mx(wk,:)).*vm(:,wnb);  % find normalized probability of each mixture for each datapoint 
            lpx(ii)=log(sum(px,1))+mx; 
            ix=jx+1; 
        end 
        l=lpx-0.5*p*log(2*pi);   % log of total probability of each data point 
    end 
end 
if nargout==0                        % attempt to plot the result 
    if p==1                            % one dimensional data           
        plot(x,l); 
    end 
end