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 values0001 function l=gmmlpdf(x,m,v,w) 0002 %GMMPDF calculated the pdf of a mixture of gaussians p=(x,m,v,w) 0003 % 0004 % Inputs: n data values, k mixtures, p parameters 0005 % 0006 % X(n,p) Input data vectors, one per row. 0007 % M(k,p) mixture means, one row per mixture. 0008 % V(k,p) mixture variances, one row per mixture (singlton dimensions will be replicated as required) 0009 % or else V(p,p,k) for full mixture covariance matrixes 0010 % W(k,1) mixture weights, one per mixture. The weights will be normalized by their sum. [default: all equal] 0011 % 0012 % Outputs: (Note that M, V and W are omitted if L==0) 0013 % 0014 % L(n,1) log PDF values 0015 0016 % Bugs/Suggestions 0017 % (1) Sort out full covariance maatrices 0018 % (2) Improve plotting 0019 % (3) Add an extra arument for plotting control 0020 0021 % Copyright (C) Mike Brookes 2000-2006 0022 % Version: $Id: gmmlpdf.m 713 2011-10-16 14:45:43Z dmb $ 0023 % 0024 % VOICEBOX is a MATLAB toolbox for speech processing. 0025 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0026 % 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 % This program is free software; you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published by 0030 % the Free Software Foundation; either version 2 of the License, or 0031 % (at your option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You can obtain a copy of the GNU General Public License from 0039 % http://www.gnu.org/copyleft/gpl.html or by writing to 0040 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 [n,p]=size(x); 0043 l=[]; % in case n=0 0044 if nargout>0 || n>0 0045 x2=x.^2; % need x^2 for variance calculation 0046 k=size(m,1); % number of mixtures 0047 if size(m,2)~=p 0048 error('x and m must have the same number of columns'); 0049 end 0050 if nargin<4 0051 if nargin<3 0052 v=1; 0053 end 0054 w=ones(k,1); 0055 end 0056 w=w/sum(w); % normalize the weights 0057 sv=size(v); 0058 if length(sv)>2 || k==1 && p>1 && sv(1)==p % full covariance matrices 0059 error('full covariance matrices not yet implemented'); 0060 else % diagonal (or constant) covariance matrices 0061 if sv(1)==1 0062 v=v(ones(k,1),:); 0063 end 0064 if sv(2)==1 0065 v=v(:,ones(1,p)); 0066 end 0067 0068 % If data size is large then do calculations in chunks 0069 0070 memsize=voicebox('memsize'); 0071 nb=min(n,max(1,floor(memsize/(8*p*k)))); % chunk size for testing data points 0072 nl=ceil(n/nb); % number of chunks 0073 0074 im=repmat(1:k,1,nb); im=im(:); 0075 0076 lpx=zeros(n,1); % log probability of each data point 0077 wk=ones(k,1); 0078 wnb=ones(1,nb); 0079 vi=v.^(-1); % calculate quantities that depend on the variances 0080 vm=sqrt(prod(vi,2)).*w; 0081 vi=-0.5*vi; 0082 0083 % first do partial chunk 0084 0085 jx=n-(nl-1)*nb; % size of first chunk 0086 ii=1:jx; 0087 kk=repmat(ii,k,1); 0088 km=repmat(1:k,1,jx); 0089 py=reshape(sum((x(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx); 0090 mx=max(py,[],1); % find normalizing factor for each data point to prevent underflow when using exp() 0091 px=exp(py-mx(wk,:)).*vm(:,ones(1,jx)); % find normalized probability of each mixture for each datapoint 0092 lpx(ii)=log(sum(px,1))+mx; 0093 ix=jx+1; 0094 0095 for il=2:nl 0096 jx=jx+nb; % increment upper limit 0097 ii=ix:jx; 0098 kk=repmat(ii,k,1); 0099 py=reshape(sum((x(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb); 0100 mx=max(py,[],1); % find normalizing factor for each data point to prevent underflow when using exp() 0101 px=exp(py-mx(wk,:)).*vm(:,wnb); % find normalized probability of each mixture for each datapoint 0102 lpx(ii)=log(sum(px,1))+mx; 0103 ix=jx+1; 0104 end 0105 l=lpx-0.5*p*log(2*pi); % log of total probability of each data point 0106 end 0107 end 0108 if nargout==0 % attempt to plot the result 0109 if p==1 % one dimensional data 0110 plot(x,l); 0111 end 0112 end