v_gausprod

PURPOSE

V_GAUSPROD calculates a product of gaussians [G,U,K]=(M,C)

SYNOPSIS

function [g,u,k]=v_gausprod(m,c,e)

DESCRIPTION

```V_GAUSPROD calculates a product of gaussians [G,U,K]=(M,C)
calculates the product of n d-dimensional multivariate gaussians
this product is itself a gaussian
Inputs: m(d,n) - each column is the mean of one of the gaussians
c(d,d,n) - contains the d#d covariance matrix for each gaussian
Alternatives: (i) c(d,n) if diagonal (ii) c(n) if c*I or (iii) omitted if I
not yet implemented:
e(d,d,n) - contains orthogonal eigenvalue matrices and c(d,n) contains eigenvalues.
Covariance matrix is E(:,:,k)*diag(c(:,k))*E(:,:,k)'
c(d,n) can then contain 0 and Inf entries

Outputs: g log gain (= log(integral) )
u(d,1) mean vector
k(d,d) or k(d) or k(1) = covariance matrix, diagonal or multiple of I (same form as input)```

CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:

SOURCE CODE

```0001 function [g,u,k]=v_gausprod(m,c,e)
0002 %V_GAUSPROD calculates a product of gaussians [G,U,K]=(M,C)
0003 % calculates the product of n d-dimensional multivariate gaussians
0004 % this product is itself a gaussian
0005 % Inputs: m(d,n) - each column is the mean of one of the gaussians
0006 %         c(d,d,n) - contains the d#d covariance matrix for each gaussian
0007 %                    Alternatives: (i) c(d,n) if diagonal (ii) c(n) if c*I or (iii) omitted if I
0008 % not yet implemented:
0009 %         e(d,d,n) - contains orthogonal eigenvalue matrices and c(d,n) contains eigenvalues.
0010 %                    Covariance matrix is E(:,:,k)*diag(c(:,k))*E(:,:,k)'
0011 %                    c(d,n) can then contain 0 and Inf entries
0012 %
0013 % Outputs: g log gain (= log(integral) )
0014 %          u(d,1) mean vector
0015 %          k(d,d) or k(d) or k(1) = covariance matrix, diagonal or multiple of I (same form as input)
0016 %
0017
0018 % Bugs: this version works with singular covariance matrices provided that their null spaces are distinct
0019 %       we could improve it slightly by doing the pseudo inverses locally and keeping track of null spaces
0020
0021 %       Copyright (C) Mike Brookes 2004
0022 %      Version: \$Id: v_gausprod.m 10865 2018-09-21 17:22:45Z dmb \$
0023 %
0024 %   VOICEBOX is a MATLAB toolbox for speech processing.
0026 %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 %   This program is free software; you can redistribute it and/or modify
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
0043 [d,n]=size(m);
0044 if nargin>2
0045     error('third argument not yet implemented in gausprod');
0046 end
0047 if nargin<2     % all covariance matrices = I
0048     c=ones(n,1);
0049 end
0050 if ~nargout     % save input data for plotting
0051     m0=m;
0052     c0=c;
0053 end
0054
0055 sc=size(c);
0056 if length(sc)<3
0057     if(sc(2)==1) & (n>1)    % covariance matrices are multiples of the identity
0058         jj=1;
0059         jj2=2;
0060         gj=zeros(n,1);
0061         while jj<n
0062             for j=1+jj:jj2:n        % we combine the gaussians in pairs
0063                 k=j-jj;
0064                 cjk=c(k)+c(j);
0065                 dm=m(:,k)-m(:,j);
0066                 gj(k)=gj(k)+gj(j)-0.5*(d*log(cjk)+dm'*dm/cjk);
0067                 m(:,k)=(c(k)*m(:,j)+c(j)*m(:,k))/cjk;
0068                 c(k)=c(k)*c(j)/cjk;
0069             end
0070             jj=jj2;
0071             jj2=2*jj;
0072         end
0073         g=gj(1);
0074         k=c(1);
0075         u=m(:,1);
0076     else                    % diagonal covariance matrices
0077         jj=1;
0078         jj2=2;
0079         gj=zeros(n,1);
0080         while jj<n
0081             for j=1+jj:jj2:n        % we combine the gaussians in pairs
0082                 k=j-jj;
0083                 cjk=c(:,k)+c(:,j);
0084                 dm=m(:,k)-m(:,j);
0085                 ix=cjk>d*max(cjk)*eps;      % calculate the psedo inverse directly
0086                 piv=zeros(d,1);
0087                 piv(ix)=cjk(ix).^(-1);
0088                 gj(k)=gj(k)+gj(j)-0.5*(log(prod(cjk))+piv'*dm.^2);
0089                 m(:,k)=piv.*(c(:,k).*m(:,j)+c(:,j).*m(:,k));
0090                 c(:,k)=c(:,k).*piv.*c(:,j);
0091             end
0092             jj=jj2;
0093             jj2=2*jj;
0094         end
0095         g=gj(1);
0096         k=c(:,1);
0097         u=m(:,1);
0098     end
0099 else                        % full covariance matrices
0100     jj=1;
0101     jj2=2;
0102     gj=zeros(n,1);
0103     while jj<n
0104         for j=1+jj:jj2:n        % we combine the gaussians in pairs
0105             k=j-jj;
0106             cjk=c(:,:,k)+c(:,:,j);
0107             dm=m(:,k)-m(:,j);
0108             piv=pinv(cjk);
0109             gj(k)=gj(k)+gj(j)-0.5*(log(det(cjk))+dm'*piv*dm);
0110             m(:,k)=c(:,:,k)*piv*m(:,j)+c(:,:,j)*piv*m(:,k);
0111             c(:,:,k)=c(:,:,k)*piv*c(:,:,j);
0112             c(:,:,k)=0.5*(c(:,:,k)+c(:,:,k)');  % ensure exactly symmetric
0113         end
0114         jj=jj2;
0115         jj2=2*jj;
0116     end
0117     g=gj(1);
0118     k=c(:,:,1);
0119     u=m(:,1);
0120 end
0121 g=g-0.5*(n-1)*d*log(2*pi);
0122
0123 if ~nargout                 % plot results if no output arguments
0124     if d==1                 % one-dimensional vectors
0125         x0=linspace(-3,3,100)';
0126         x=zeros(length(x0),n);
0127         y=x;
0128         for j=1:n
0129             x(:,j)=x0+m0(1,j);
0130             cj=c0(j);
0131             y(:,j)=normpdf(x0,0,sqrt(cj));
0132         end
0133         plot(x,log10(y),':',x0+u,log10(normpdf(x0,0,k)*exp(g)),'k-');
0134         ylabel('Log10(pdf)');
0135     else
0136         if length(sc)<3
0137             if(sc(2)==1) & (n>1)    % covariance matrices are multiples of the identity
0138                 sk=k*eye(d);
0139             else                    % diagonal covariance matrices
0140                 sk=diag(k);
0141             end
0142             uk=eye(d);
0143             vk=uk;
0144         else                        % full covariance matrices
0145             [uk,sk,vk]=svd(k);
0146         end
0147
0148
0149         u2=uk(:,1:2);
0150         t0=linspace(0,2*pi,100);
0151         x=zeros(length(t0),n);
0152         y=x;
0153         x0=[cos(t0); sin(t0)];
0154         for j=1:n
0155             if length(sc)<3
0156                 if(sc(2)==1) & (n>1)    % covariance matrices are multiples of the identity
0157                     cj=c0(j)*eye(2);
0158                 else                    % diagonal covariance matrices
0159                     cj=u2'*diag(c0(:,j))*u2;
0160                 end
0161             else                        % full covariance matrices
0162                 cj=u2'*c0(:,:,j)*u2;
0163             end
0164             mj=u2'*m0(:,j);
0165             v=sqrt(sum((x0'/cj).*x0',2).^(-1));
0166             x(:,j)=mj(1)+v.*x0(1,:)';
0167             y(:,j)=mj(2)+v.*x0(2,:)';
0168         end
0169
0170         if length(sc)<3
0171             if(sc(2)==1) & (n>1)    % covariance matrices are multiples of the identity
0172                 cj=k*eye(2);
0173             else                    % diagonal covariance matrices
0174                 cj=u2'*diag(k)*u2;
0175             end
0176         else                        % full covariance matrices
0177             cj=u2'*k*u2;
0178         end
0179         mj=u2'*u;
0180         v=sqrt(sum((x0'/cj).*x0',2).^(-1));
0181         plot(x,y,':',mj(1)+v.*x0(1,:)',mj(2)+v.*x0(2,:)','k-');
0182         axis equal;
0183     end
0184 end```