# v_gaussmixd

## PURPOSE

V_GAUSSMIXD marginal and conditional Gaussian mixture densities

## SYNOPSIS

function [mz,vz,wz]=v_gaussmixd(y,m,v,w,a,b,f,g)

## DESCRIPTION

    V_GAUSSMIXD marginal and conditional Gaussian mixture densities

     Usage: (1) [mz,vz,wz]=v_gaussmixd(y,m,v,w); % If {m,v,w} has dimension p, then y specifies
                                               % elements 1:q and {mz,vz,wz} is a GMM for elements q+1:p
            (2) [mz,vz,wz]=v_gaussmixd(y,m,v,w,[],[1 3 4]);    % y specifies elements 1,3,4 and {mz,vz,wz}
                                                           % is a GMM for the remaining elements

     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) = mixture weights
       A(q,p) = conditioning transformation: y=x*A'+ B' (where y and x are row vectors).
       B(q,1)   If A is omitted or null, y=x*I(B,:)' where I is the identity matrix.
                If B is also omitted or null, y=x*I(1:q,:)'.
       F(r,p) = output transformation z = x*F'+G'  (where z and x are row vectors).
       G(r,1)   If F is omitted or null, z = x*f' where I is the identity matrix.
                If G is also omitted or null, z=x*I(q+1:p,:)' or, if A was also null,
                the complement of y.

     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); surprisingly it is 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=v_gaussmixd(y,m,v,w)
     will return the estimated values of the remaining p-q components.

     See also: v_gaussmix, v_gaussmixg, v_gaussmixp, v_randvec

## CROSS-REFERENCE INFORMATION

This function calls:

- v_gaussmixp V_GAUSSMIXP calculate probability densities from or plot a Gaussian mixture model

This function is called by:

- v_gaussmixg V_GAUSSMIXG global mean, variance and mode of a GMM

## SOURCE CODE

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