wolffd@0: function [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat, unit_norm)
wolffd@0: % EVAL_PDF_COND_MOG Evaluate the pdf of a conditional mixture of Gaussians
wolffd@0: % function [B, B2] = eval_pdf_cond_mog(data, mu, Sigma, mixmat, unit_norm)
wolffd@0: %
wolffd@0: % Notation: Y is observation, M is mixture component, and both may be conditioned on Q.
wolffd@0: % If Q does not exist, ignore references to Q=j below.
wolffd@0: % Alternatively, you may ignore M if this is a conditional Gaussian.
wolffd@0: %
wolffd@0: % INPUTS:
wolffd@0: % data(:,t) = t'th observation vector 
wolffd@0: %
wolffd@0: % mu(:,k) = E[Y(t) | M(t)=k] 
wolffd@0: % or mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k]
wolffd@0: %
wolffd@0: % Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k]
wolffd@0: % or there are various faster, special cases:
wolffd@0: %   Sigma() - scalar, spherical covariance independent of M,Q.
wolffd@0: %   Sigma(:,:) diag or full, tied params independent of M,Q. 
wolffd@0: %   Sigma(:,:,j) tied params independent of M. 
wolffd@0: %
wolffd@0: % mixmat(k) = Pr(M(t)=k) = prior
wolffd@0: % or mixmat(j,k) = Pr(M(t)=k | Q(t)=j) 
wolffd@0: % Not needed if M is not defined.
wolffd@0: %
wolffd@0: % unit_norm - optional; if 1, means data(:,i) AND mu(:,i) each have unit norm (slightly faster)
wolffd@0: %
wolffd@0: % OUTPUT:
wolffd@0: % B(t) = Pr(y(t)) 
wolffd@0: % or
wolffd@0: % B(i,t) = Pr(y(t) | Q(t)=i) 
wolffd@0: % B2(i,k,t) = Pr(y(t) | Q(t)=i, M(t)=k) 
wolffd@0: %
wolffd@0: % If the number of mixture components differs depending on Q, just set the trailing
wolffd@0: % entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2,
wolffd@0: % then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0.
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: if isvectorBNT(mu) & size(mu,2)==1
wolffd@0:   d = length(mu);
wolffd@0:   Q = 1; M = 1;
wolffd@0: elseif ndims(mu)==2
wolffd@0:   [d Q] = size(mu);
wolffd@0:   M = 1;
wolffd@0: else
wolffd@0:   [d Q M] = size(mu);
wolffd@0: end
wolffd@0: [d T] = size(data);
wolffd@0: 
wolffd@0: if nargin < 4, mixmat = ones(Q,1); end
wolffd@0: if nargin < 5, unit_norm = 0; end
wolffd@0: 
wolffd@0: %B2 = zeros(Q,M,T); % ATB: not needed allways
wolffd@0: %B = zeros(Q,T);
wolffd@0: 
wolffd@0: if isscalarBNT(Sigma)
wolffd@0:   mu = reshape(mu, [d Q*M]);
wolffd@0:   if unit_norm % (p-q)'(p-q) = p'p + q'q - 2p'q = n+m -2p'q since p(:,i)'p(:,i)=1
wolffd@0:     %avoid an expensive repmat
wolffd@0:     disp('unit norm')
wolffd@0:     %tic; D = 2 -2*(data'*mu)'; toc 
wolffd@0:     D = 2 - 2*(mu'*data);
wolffd@0:     tic; D2 = sqdist(data, mu)'; toc
wolffd@0:     assert(approxeq(D,D2)) 
wolffd@0:   else
wolffd@0:     D = sqdist(data, mu)';
wolffd@0:   end
wolffd@0:   clear mu data % ATB: clear big old data
wolffd@0:   % D(qm,t) = sq dist between data(:,t) and mu(:,qm)
wolffd@0:   logB2 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d
wolffd@0:   B2 = reshape(exp(logB2), [Q M T]);
wolffd@0:   clear logB2 % ATB: clear big old data
wolffd@0:   
wolffd@0: elseif ndims(Sigma)==2 % tied full
wolffd@0:   mu = reshape(mu, [d Q*M]);
wolffd@0:   D = sqdist(data, mu, inv(Sigma))';
wolffd@0:   % D(qm,t) = sq dist between data(:,t) and mu(:,qm)
wolffd@0:   logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D;
wolffd@0:   %denom = sqrt(det(2*pi*Sigma));
wolffd@0:   %numer = exp(-0.5 * D);
wolffd@0:   %B2 = numer/denom;
wolffd@0:   B2 = reshape(exp(logB2), [Q M T]);
wolffd@0:   
wolffd@0: elseif ndims(Sigma)==3 % tied across M
wolffd@0:   B2 = zeros(Q,M,T);
wolffd@0:   for j=1:Q
wolffd@0:     % D(m,t) = sq dist between data(:,t) and mu(:,j,m)
wolffd@0:     if isposdef(Sigma(:,:,j))
wolffd@0:       D = sqdist(data, permute(mu(:,j,:), [1 3 2]), inv(Sigma(:,:,j)))';
wolffd@0:       logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma(:,:,j)) - 0.5*D;
wolffd@0:       B2(j,:,:) = exp(logB2);
wolffd@0:     else
wolffd@0:       error(sprintf('mixgauss_prob: Sigma(:,:,q=%d) not psd\n', j));
wolffd@0:     end
wolffd@0:   end
wolffd@0:   
wolffd@0: else % general case
wolffd@0:   B2 = zeros(Q,M,T);
wolffd@0:   for j=1:Q
wolffd@0:     for k=1:M
wolffd@0:       %if mixmat(j,k) > 0
wolffd@0:       B2(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma(:,:,j,k));
wolffd@0:       %end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % B(j,t) = sum_k B2(j,k,t) * Pr(M(t)=k | Q(t)=j) 
wolffd@0: 
wolffd@0: % The repmat is actually slower than the for-loop, because it uses too much memory
wolffd@0: % (this is true even for small T).
wolffd@0: 
wolffd@0: %B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2));
wolffd@0: %B = reshape(B, [Q T]); % undo effect of squeeze in case Q = 1
wolffd@0:   
wolffd@0: B = zeros(Q,T);
wolffd@0: if Q < T
wolffd@0:   for q=1:Q
wolffd@0:     %B(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % squeeze chnages order if M=1
wolffd@0:     B(q,:) = mixmat(q,:) * permute(B2(q,:,:), [2 3 1]); % vector * matrix sums over m
wolffd@0:   end
wolffd@0: else
wolffd@0:   for t=1:T
wolffd@0:     B(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m
wolffd@0:   end
wolffd@0: end
wolffd@0: %t=toc;fprintf('%5.3f\n', t)
wolffd@0: 
wolffd@0: %tic
wolffd@0: %A = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2));
wolffd@0: %t=toc;fprintf('%5.3f\n', t)
wolffd@0: %assert(approxeq(A,B)) % may be false because of round off error