Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMstats/mixgauss_prob.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat, unit_norm)
wolffd@0 2 % EVAL_PDF_COND_MOG Evaluate the pdf of a conditional mixture of Gaussians
wolffd@0 3 % function [B, B2] = eval_pdf_cond_mog(data, mu, Sigma, mixmat, unit_norm)
wolffd@0 4 %
wolffd@0 5 % Notation: Y is observation, M is mixture component, and both may be conditioned on Q.
wolffd@0 6 % If Q does not exist, ignore references to Q=j below.
wolffd@0 7 % Alternatively, you may ignore M if this is a conditional Gaussian.
wolffd@0 8 %
wolffd@0 9 % INPUTS:
wolffd@0 10 % data(:,t) = t'th observation vector
wolffd@0 11 %
wolffd@0 12 % mu(:,k) = E[Y(t) | M(t)=k]
wolffd@0 13 % or mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k]
wolffd@0 14 %
wolffd@0 15 % Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k]
wolffd@0 16 % or there are various faster, special cases:
wolffd@0 17 % Sigma() - scalar, spherical covariance independent of M,Q.
wolffd@0 18 % Sigma(:,:) diag or full, tied params independent of M,Q.
wolffd@0 19 % Sigma(:,:,j) tied params independent of M.
wolffd@0 20 %
wolffd@0 21 % mixmat(k) = Pr(M(t)=k) = prior
wolffd@0 22 % or mixmat(j,k) = Pr(M(t)=k | Q(t)=j)
wolffd@0 23 % Not needed if M is not defined.
wolffd@0 24 %
wolffd@0 25 % unit_norm - optional; if 1, means data(:,i) AND mu(:,i) each have unit norm (slightly faster)
wolffd@0 26 %
wolffd@0 27 % OUTPUT:
wolffd@0 28 % B(t) = Pr(y(t))
wolffd@0 29 % or
wolffd@0 30 % B(i,t) = Pr(y(t) | Q(t)=i)
wolffd@0 31 % B2(i,k,t) = Pr(y(t) | Q(t)=i, M(t)=k)
wolffd@0 32 %
wolffd@0 33 % If the number of mixture components differs depending on Q, just set the trailing
wolffd@0 34 % entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2,
wolffd@0 35 % then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0.
wolffd@0 36
wolffd@0 37
wolffd@0 38
wolffd@0 39
wolffd@0 40 if isvectorBNT(mu) & size(mu,2)==1
wolffd@0 41 d = length(mu);
wolffd@0 42 Q = 1; M = 1;
wolffd@0 43 elseif ndims(mu)==2
wolffd@0 44 [d Q] = size(mu);
wolffd@0 45 M = 1;
wolffd@0 46 else
wolffd@0 47 [d Q M] = size(mu);
wolffd@0 48 end
wolffd@0 49 [d T] = size(data);
wolffd@0 50
wolffd@0 51 if nargin < 4, mixmat = ones(Q,1); end
wolffd@0 52 if nargin < 5, unit_norm = 0; end
wolffd@0 53
wolffd@0 54 %B2 = zeros(Q,M,T); % ATB: not needed allways
wolffd@0 55 %B = zeros(Q,T);
wolffd@0 56
wolffd@0 57 if isscalarBNT(Sigma)
wolffd@0 58 mu = reshape(mu, [d Q*M]);
wolffd@0 59 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 60 %avoid an expensive repmat
wolffd@0 61 disp('unit norm')
wolffd@0 62 %tic; D = 2 -2*(data'*mu)'; toc
wolffd@0 63 D = 2 - 2*(mu'*data);
wolffd@0 64 tic; D2 = sqdist(data, mu)'; toc
wolffd@0 65 assert(approxeq(D,D2))
wolffd@0 66 else
wolffd@0 67 D = sqdist(data, mu)';
wolffd@0 68 end
wolffd@0 69 clear mu data % ATB: clear big old data
wolffd@0 70 % D(qm,t) = sq dist between data(:,t) and mu(:,qm)
wolffd@0 71 logB2 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d
wolffd@0 72 B2 = reshape(exp(logB2), [Q M T]);
wolffd@0 73 clear logB2 % ATB: clear big old data
wolffd@0 74
wolffd@0 75 elseif ndims(Sigma)==2 % tied full
wolffd@0 76 mu = reshape(mu, [d Q*M]);
wolffd@0 77 D = sqdist(data, mu, inv(Sigma))';
wolffd@0 78 % D(qm,t) = sq dist between data(:,t) and mu(:,qm)
wolffd@0 79 logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D;
wolffd@0 80 %denom = sqrt(det(2*pi*Sigma));
wolffd@0 81 %numer = exp(-0.5 * D);
wolffd@0 82 %B2 = numer/denom;
wolffd@0 83 B2 = reshape(exp(logB2), [Q M T]);
wolffd@0 84
wolffd@0 85 elseif ndims(Sigma)==3 % tied across M
wolffd@0 86 B2 = zeros(Q,M,T);
wolffd@0 87 for j=1:Q
wolffd@0 88 % D(m,t) = sq dist between data(:,t) and mu(:,j,m)
wolffd@0 89 if isposdef(Sigma(:,:,j))
wolffd@0 90 D = sqdist(data, permute(mu(:,j,:), [1 3 2]), inv(Sigma(:,:,j)))';
wolffd@0 91 logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma(:,:,j)) - 0.5*D;
wolffd@0 92 B2(j,:,:) = exp(logB2);
wolffd@0 93 else
wolffd@0 94 error(sprintf('mixgauss_prob: Sigma(:,:,q=%d) not psd\n', j));
wolffd@0 95 end
wolffd@0 96 end
wolffd@0 97
wolffd@0 98 else % general case
wolffd@0 99 B2 = zeros(Q,M,T);
wolffd@0 100 for j=1:Q
wolffd@0 101 for k=1:M
wolffd@0 102 %if mixmat(j,k) > 0
wolffd@0 103 B2(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma(:,:,j,k));
wolffd@0 104 %end
wolffd@0 105 end
wolffd@0 106 end
wolffd@0 107 end
wolffd@0 108
wolffd@0 109 % B(j,t) = sum_k B2(j,k,t) * Pr(M(t)=k | Q(t)=j)
wolffd@0 110
wolffd@0 111 % The repmat is actually slower than the for-loop, because it uses too much memory
wolffd@0 112 % (this is true even for small T).
wolffd@0 113
wolffd@0 114 %B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2));
wolffd@0 115 %B = reshape(B, [Q T]); % undo effect of squeeze in case Q = 1
wolffd@0 116
wolffd@0 117 B = zeros(Q,T);
wolffd@0 118 if Q < T
wolffd@0 119 for q=1:Q
wolffd@0 120 %B(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % squeeze chnages order if M=1
wolffd@0 121 B(q,:) = mixmat(q,:) * permute(B2(q,:,:), [2 3 1]); % vector * matrix sums over m
wolffd@0 122 end
wolffd@0 123 else
wolffd@0 124 for t=1:T
wolffd@0 125 B(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m
wolffd@0 126 end
wolffd@0 127 end
wolffd@0 128 %t=toc;fprintf('%5.3f\n', t)
wolffd@0 129
wolffd@0 130 %tic
wolffd@0 131 %A = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2));
wolffd@0 132 %t=toc;fprintf('%5.3f\n', t)
wolffd@0 133 %assert(approxeq(A,B)) % may be false because of round off error