# v_gaussmixm

## PURPOSE

V_GAUSSMIXM estimate mean and variance of the magnitude of a GMM

## SYNOPSIS

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

## DESCRIPTION

    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.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    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 Lesser General Public License as published by
    0029 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0036 %
    0037 %   You can obtain a copy of the GNU Lesser General Public License from
    0038 %   https://www.gnu.org/licenses/ .
    0039 %
    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