V_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.
0001 function [mm,mc]=v_gaussmixm(m,v,w,z) 0002 %V_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: v_gaussmixm.m 10865 2018-09-21 17:22:45Z 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