GAUSSMIXD marginal and conditional Gaussian mixture densities
Inputs: Input mixture: k mixtures, p dimensions
Output mixture: k mixtures, r dimensions
Conditioning data: n data values, q dimensions
Y(n,q) = conditioning input data: x*a'+ b'= y
M(k,p) = mixture means for x(p)
V(k,p) or V(p,p,k) variances (diagonal or full)
W(k,1)
A(q,p) = conditioning transformation: y=x*a'+ b' (where y and x are row vectors)
if A is omitted, it is assumed to be the first q rows of the
identity matrix.
B(q,1) = part of condition transform [default B=zeros(q,1)]
F(r,p) = output transformation z = x*f'+g' (where z and x are row vectors)
if F is omitted, it is assumed to be the last r=p-q rows of the
identity matrix.
G(r,1) = part of output transformation
Outputs (if 1 or 2 outputs specified):
MZ(n,r) = Global mean of z for each y
VZ(r,r,n) Global full covariances of z (now dependent on y)
Outputs (if 3 outputs specified):
MZ(k,r,n) = mixture means of z for each y
VZ(k,r) or VZ(r,r,k) variances of z (diagonal or full) independent of y
WZ(n,k)
The output mixture covariances will be diagonal if either r=1 or else the following three
conditions are all true:
(a) the input mixture covariances are diagonal and
(b) matrix A has at most one non-zero element in any row or column and
(c) matrix F has at most one non-zero element in any column
Several of the output variables can be squeezed if r=1 but this is not done automatically.
This routine can be used for inference: if you train a GMM on p-dimensional data
then, if y(n,q) contains observations of the first q components, then z=gaussmixd(y,m,v,w)
will return the estimated values of the remaining p-q components.0001 function [mz,vz,wz]=gaussmixd(y,m,v,w,a,b,f,g) 0002 %GAUSSMIXD marginal and conditional Gaussian mixture densities 0003 % 0004 % Inputs: Input mixture: k mixtures, p dimensions 0005 % Output mixture: k mixtures, r dimensions 0006 % Conditioning data: n data values, q dimensions 0007 % 0008 % Y(n,q) = conditioning input data: x*a'+ b'= y 0009 % M(k,p) = mixture means for x(p) 0010 % V(k,p) or V(p,p,k) variances (diagonal or full) 0011 % W(k,1) 0012 % A(q,p) = conditioning transformation: y=x*a'+ b' (where y and x are row vectors) 0013 % if A is omitted, it is assumed to be the first q rows of the 0014 % identity matrix. 0015 % B(q,1) = part of condition transform [default B=zeros(q,1)] 0016 % F(r,p) = output transformation z = x*f'+g' (where z and x are row vectors) 0017 % if F is omitted, it is assumed to be the last r=p-q rows of the 0018 % identity matrix. 0019 % G(r,1) = part of output transformation 0020 % 0021 % Outputs (if 1 or 2 outputs specified): 0022 % 0023 % MZ(n,r) = Global mean of z for each y 0024 % VZ(r,r,n) Global full covariances of z (now dependent on y) 0025 % 0026 % Outputs (if 3 outputs specified): 0027 % 0028 % MZ(k,r,n) = mixture means of z for each y 0029 % VZ(k,r) or VZ(r,r,k) variances of z (diagonal or full) independent of y 0030 % WZ(n,k) 0031 % 0032 % The output mixture covariances will be diagonal if either r=1 or else the following three 0033 % conditions are all true: 0034 % (a) the input mixture covariances are diagonal and 0035 % (b) matrix A has at most one non-zero element in any row or column and 0036 % (c) matrix F has at most one non-zero element in any column 0037 % 0038 % Several of the output variables can be squeezed if r=1 but this is not done automatically. 0039 % 0040 % This routine can be used for inference: if you train a GMM on p-dimensional data 0041 % then, if y(n,q) contains observations of the first q components, then z=gaussmixd(y,m,v,w) 0042 % will return the estimated values of the remaining p-q components. 0043 0044 % Copyright (C) Mike Brookes 2000-2009 0045 % Version: $Id: gaussmixd.m 713 2011-10-16 14:45:43Z dmb $ 0046 % 0047 % VOICEBOX is a MATLAB toolbox for speech processing. 0048 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0049 % 0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0051 % This program is free software; you can redistribute it and/or modify 0052 % it under the terms of the GNU General Public License as published by 0053 % the Free Software Foundation; either version 2 of the License, or 0054 % (at your option) any later version. 0055 % 0056 % This program is distributed in the hope that it will be useful, 0057 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0058 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0059 % GNU General Public License for more details. 0060 % 0061 % You can obtain a copy of the GNU General Public License from 0062 % http://www.gnu.org/copyleft/gpl.html or by writing to 0063 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0065 0066 [n,q]=size(y); 0067 [k,p]=size(m); 0068 % We set fv=1 if input V contains full covariance matrices with dimensions V(r,r,k). 0069 % This is true either if V has >2 dimensions or V=V(r,r) and r>k 0070 fv=ndims(v)>2 || size(v,1)>k; % full covariance matrix is supplied 0071 if nargin<7 || isempty(f) 0072 f=eye(p); 0073 f=f(p+1:end,:); 0074 end 0075 if nargin<5 || isempty(a) 0076 a=eye(p); 0077 a=a(1:p,:); 0078 end 0079 r=size(f,1); 0080 if nargin<8 || isempty(g) 0081 g=zeros(r,1); 0082 end 0083 if nargin<6 || isempty(b) 0084 b=zeros(q,1); 0085 end 0086 yb=y-repmat(b,n,1); % remove the b term 0087 [lp,wz]=gaussmixp(yb,m,v,w,a); % find mixture weights 0088 mz=zeros(n,r,k); 0089 ma=a~=0; % check for sparse a and f matrices (common case) 0090 mf=f~=0; 0091 % We set dvo=1 if the output mixture covariances are structurally 0092 % diagonal. This is the case if eith 0093 % (1) r=1 since in this case they are scalar values, or else 0094 % (2) the following three conditions are all true: 0095 % (a) the input mixture covariances are diagonal and 0096 % (b) matrix A has at most one non-zero element in any row or column and 0097 % (c) matrix F has at most one non-zero element in any column 0098 dvo=r==1 || (~fv && all(sum(ma,1)<=1) && all(sum(ma,2)<=1) && all(sum(mf,1)<=1)); 0099 if dvo % structurally diagonal output covariances 0100 vz=zeros(k,r); 0101 for i=1:k 0102 if fv 0103 vi=v(:,:,i); 0104 else 0105 vi=diag(v(i,:)); % convert to full covariance matrices 0106 end 0107 hi=vi*a'/(a*vi*a'); 0108 vz(i,:)=diag(f*(vi-hi*a*vi)*f')'; 0109 mi=m(i,:); 0110 m0=(mi-mi*a'*hi')+g; 0111 mz(:,:,i)=(repmat(m0,n,1)+yb*hi')*f'; 0112 end 0113 else 0114 vz=zeros(r,r,k); 0115 for i=1:k 0116 if fv 0117 vi=v(:,:,i); 0118 else 0119 vi=diag(v(i,:)); % convert to full covariance matrices 0120 end 0121 hi=vi*a'/(a*vi*a'); 0122 vz(:,:,i)=f*(vi-hi*a*vi)*f'; 0123 mi=m(i,:); 0124 m0=(mi-mi*a'*hi')+g; 0125 mz(:,:,i)=(repmat(m0,n,1)+yb*hi')*f'; 0126 end 0127 end 0128 if nargout<3 0129 mt=reshape(sum(reshape(mz,n*r,k).*repmat(wz,r,1),2),n,r); % global mean 0130 if nargout>1 % need to calculate global varaince as well 0131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0132 % We calculate the global covariance by adding up the weighted mixture 0133 % covariances corrected for the fact that the mixture mean does not equal 0134 % the global mean. 0135 % To save calculations, we calculate only the lower triangle of the 0136 % symmetric covariance matrix and then expand it at the end into a full matrix. 0137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0138 rl=r*(r+1)/2; % number of elements in lower triangular covariance matrix 0139 lix=1:r^2; 0140 cix=repmat(1:r,r,1); 0141 rix=cix'; 0142 lix(cix>rix)=[]; % index of lower triangular elements 0143 rlix=rix(lix); 0144 clix=cix(lix); 0145 lixi=zeros(r,r); 0146 lixi(lix)=1:rl; 0147 lixi=lixi'; 0148 lixi(lix)=1:rl; % reverse index to build full matrices 0149 0150 vt=zeros(n,rl); % reserve space for lower triangular output covariances 0151 for i=1:k 0152 if dvo 0153 vi=diag(vz(i,:)); 0154 else 0155 vi=vz(:,:,i); 0156 end 0157 mzt=mz(:,:,i)-mt; 0158 vt=vt+repmat(wz(:,i),1,r^2).*(repmat(vi(lix),n,1)+mzt(:,rlix).*mzt(:,clix)); 0159 end 0160 mz=mt; 0161 vz=permute(reshape(vt(:,lixi),[n,r,r]),[2 3 1]); 0162 end 0163 else 0164 mz=permute(mz,3:-1:1); 0165 end