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