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