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)
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 General Public License as published by 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