Home > voicebox > gaussmixm.m

gaussmixm

PURPOSE ^

GAUSSMIXM estimate mean and variance of the magnitude of a GMM

SYNOPSIS ^

function [mm,mc]=gaussmixm(m,v,w,z)

DESCRIPTION ^

 GAUSSMIXM estimate mean and variance of the magnitude of a GMM

  Inputs:  M(K,P)   is the mean vectors (one row per mixture)
           V(K,P)   are diagonal covariances (one row per mixture) [ones(K,P)]
        or V(P,P,K) are full covariance matrices (one per mixture)
           W(K)     are the mixture weights [ones(K,1)/K]
           Z(T,P)   each row is an origin to measure magnitude from [zeros(1,P)]

 Outputs:  MM(T,1)  mean of sqrt((x-z(t))'*(x-z(t))) where x is the GMM random variate
           MC(T,T)  covariance matrix of sqrt((x-z(t))'*(x-z(t)))

 We approximate the magnitude of each mixture as a Nakagami-m distribution and
 match the moments of its squared variate. We approximate the normalized 
 correlation matrix of |x| (i.e. with unit diagonal) to be the same as that
 of |x|^2 which we can calculate exactly.
 Answers are exact for P=1 or when all M=0. Accuracy improves otherwise for
 large P or M.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [mm,mc]=gaussmixm(m,v,w,z)
0002 % GAUSSMIXM estimate mean and variance of the magnitude of a GMM
0003 %
0004 %  Inputs:  M(K,P)   is the mean vectors (one row per mixture)
0005 %           V(K,P)   are diagonal covariances (one row per mixture) [ones(K,P)]
0006 %        or V(P,P,K) are full covariance matrices (one per mixture)
0007 %           W(K)     are the mixture weights [ones(K,1)/K]
0008 %           Z(T,P)   each row is an origin to measure magnitude from [zeros(1,P)]
0009 %
0010 % Outputs:  MM(T,1)  mean of sqrt((x-z(t))'*(x-z(t))) where x is the GMM random variate
0011 %           MC(T,T)  covariance matrix of sqrt((x-z(t))'*(x-z(t)))
0012 %
0013 % We approximate the magnitude of each mixture as a Nakagami-m distribution and
0014 % match the moments of its squared variate. We approximate the normalized
0015 % correlation matrix of |x| (i.e. with unit diagonal) to be the same as that
0016 % of |x|^2 which we can calculate exactly.
0017 % Answers are exact for P=1 or when all M=0. Accuracy improves otherwise for
0018 % large P or M.
0019 
0020 %      Copyright (C) Mike Brookes 2015
0021 %      Version: $Id: gaussmixm.m 5957 2015-03-26 16:59:23Z dmb $
0022 %
0023 %   VOICEBOX is a MATLAB toolbox for speech processing.
0024 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0025 %
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 %   This program is free software; you can redistribute it and/or modify
0028 %   it under the terms of the GNU General Public License as published by
0029 %   the Free Software Foundation; either version 2 of the License, or
0030 %   (at your option) any later version.
0031 %
0032 %   This program is distributed in the hope that it will be useful,
0033 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0034 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0035 %   GNU General Public License for more details.
0036 %
0037 %   You can obtain a copy of the GNU General Public License from
0038 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0039 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0041 [k,p]=size(m);                          % k = # mixtures, p = vector dimension
0042 if nargin<4 || ~numel(z)
0043     z=zeros(1,p);
0044 end
0045 if nargin<3 || ~numel(w)
0046     w=ones(k,1);
0047 end% default to uniform weights
0048 if nargin<2 || ~numel(v)
0049     v=ones(k,p);                    % default to unit variance
0050 end
0051 [t,p1]=size(z);
0052 if p~=p1 || (p>2 && t>1 && nargout>2)
0053     error('unable to calculate a covariance matrix');
0054 end
0055 w=w(:)/sum(w);                          % make the weights sum to 1 and force column vector
0056 if p==1                                 % treat 1D case specially since we have an exact formula
0057     s=sqrt(v(:));                       % normalize mixture means by the std dev
0058     mt=m(:,ones(1,t))-z(:,ones(1,k))';  % shifted mean for each mixture x origin
0059     mts=mt./s(:,ones(1,t));             % shifted mean for each mixture x origin normalized to unit SD
0060     ncdf=normcdf(-mts);
0061     npdf=normpdf(-mts);
0062     mm=(mts.*(1-2*ncdf)+2*npdf)'*(s.*w); % exact mean of |X|
0063     mc=diag((mts.^2+1)'*(v(:).*w)); % diagonal variance elements in case t==1
0064     if nargout>1 && t>1                 % we need to calculate a covariance matrix
0065         for it=1:t
0066             for jt=1:it-1
0067                 mc(it,jt)=w.'*((v(:)+mt(:,it).*mt(:,jt)).*(1-2*abs(ncdf(:,it)-ncdf(:,jt)))+ ...
0068                     2*s.*sign(mt(:,jt)-mt(:,it)).*(npdf(:,it).*mt(:,jt)-npdf(:,jt).*mt(:,it)));
0069                 mc(jt,it)=mc(it,jt);
0070             end
0071         end   
0072     end
0073     mc=mc-mm*mm';                       % convert to variance
0074 else                                        % p>1 case
0075     fullv=ndims(v)>2 || size(v,1)>k;          % full covariance matrix is supplied
0076     if fullv                                % full covariance matrix is supplied
0077         ms=repmat(sum(m.^2+v(repmat((1:p+1:p^2)',1,k)+repmat(0:p^2:(k-1)*p^2,p,1))',2),1,t)-2*m*z'+repmat(sum(z.^2,2)',k,1);  % mean of |x|^2 for each mixture x origin
0078         vsf=zeros(t,t,k);
0079         vs=zeros(k,t);
0080         for i=1:k
0081             zmi=repmat(m(i,:),t,1)-z;
0082             si=v(:,:,i);
0083             vsi=2*trace(si^2)+4*zmi*si*zmi';
0084             vsf(:,:,i)=vsi;
0085             vs(i,:)=diag(vsi);
0086         end
0087     else                                        % else diagonal covariance matrix supplied
0088         ms=repmat(sum(m.^2,2)+sum(v,2),1,t)-2*m*z'+repmat(sum(z.^2,2)',k,1);  % mean of |x|^2 for each mixture x origin
0089         vsc=sum(v.*(2*v+4*m.^2),2);             % origin-independent part of  |x|^2 variance for each mixture
0090         vmz=4*(v.*m)*z';                        % origin-linear part of  |x|^2 variannce for each mixture x origin
0091         vs=repmat(vsc,1,t)-2*vmz+4*v*z.^2';     % variance of |x|^2 for each mixture x origin
0092     end
0093     nm=ms.^2./vs;                               % Nakagami-m parameter per mixture x origin
0094     mmk=(exp(gammaln(nm+0.5)-gammaln(nm)).*sqrt(ms./nm)); % mean of Nakagami-m distrbution per mixture x origin
0095     mm=mmk'*w;                                  % mean per origin
0096     mc = ms'*w-mm.^2 ;                          % variance in case t==1
0097     if nargout>1 && t>1                         % we need to calculate the t x t covariance matrix w.r.t. origins
0098         mvksd=sqrt(ms-mmk.^2);                  % std deviation of |x| per mixture x origin
0099         mc=zeros(t,t);
0100         if fullv                                % full covariance matrix is supplied
0101             for i=1:k
0102                 vsisd=sqrt(vs(i,:));
0103                 mcm=(mvksd(i,:)'*mvksd(i,:)).*vsf(:,:,i)./(vsisd'*vsisd)+ mmk(i,:)'*mmk(i,:); % estimated |x| correlation matrix for mixture i
0104                 mc=mc+w(i)*mcm;             % add onto total correlation matrix
0105             end
0106         else                                % else diagonal covariance matrix supplied
0107             for i=1:k                       % loop for each mixture component
0108                 iit=repmat(i,t,1);
0109                 vsi=vsc(i)-vmz(iit,:)-vmz(iit,:)'+4*z.*v(iit,:)*z'; % covariance matriz for |x|^2 for mixture i
0110                 vsisd=sqrt(vs(i,:));
0111                 mcm=(mvksd(i,:)'*mvksd(i,:)).*vsi./(vsisd'*vsisd)+ mmk(i,:)'*mmk(i,:); % estimated |x| correlation matrix for mixture i
0112                 mc=mc+w(i)*mcm;   % add onto total correlation matrix
0113             end
0114         end
0115         mc=mc-mm*mm'; % convert to covariance matrix
0116     end
0117 end

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003