wolffd@0: function test_eval_pdf_cond_mixgauss()
wolffd@0: 
wolffd@0: %Q = 10; M = 100; d = 20; T = 500;
wolffd@0: Q = 2; M = 3; d = 4; T = 5;
wolffd@0: 
wolffd@0: mu = rand(d,Q,M);
wolffd@0: data = randn(d,T);
wolffd@0: %mixmat = mk_stochastic(rand(Q,M));
wolffd@0: mixmat = mk_stochastic(ones(Q,M));
wolffd@0: 
wolffd@0: % tied scalar
wolffd@0: Sigma = 0.01;
wolffd@0: 
wolffd@0: mu = rand(d,M,Q);
wolffd@0: weights = mixmat';
wolffd@0: N = M*ones(1,Q);
wolffd@0: tic; [B, B2, D] = parzen(data, mu, Sigma, N, weights); toc
wolffd@0: tic; [BC, B2C, DC] = parzenC(data, mu, Sigma, N); toc
wolffd@0: approxeq(B,BC)
wolffd@0: B2C = reshape(B2C,[M Q T]);
wolffd@0: approxeq(B2,B2C)
wolffd@0: DC = reshape(DC,[M Q T]);
wolffd@0: approxeq(D,DC)
wolffd@0: 
wolffd@0: 
wolffd@0: return
wolffd@0: 
wolffd@0: tic; [B, B2] = eval_pdf_cond_mixgauss(data, mu, Sigma, mixmat); toc
wolffd@0: tic; C = eval_pdf_cond_parzen(data, mu, Sigma); toc
wolffd@0: approxeq(B,C)
wolffd@0: 
wolffd@0: return;
wolffd@0: 
wolffd@0: 
wolffd@0: mu = reshape(mu, [d Q*M]);
wolffd@0: 
wolffd@0: data = mk_unit_norm(data);
wolffd@0: mu = mk_unit_norm(mu);
wolffd@0: tic; D = 2 -2*(data'*mu); toc % avoid an expensive repmat
wolffd@0: tic; D2 = sqdist(data, mu); toc
wolffd@0: approxeq(D,D2)
wolffd@0: 
wolffd@0: 
wolffd@0: % D(t,m) = sq dist between data(:,t) and mu(:,m)
wolffd@0: mu = reshape(mu, [d Q*M]);
wolffd@0: D = dist2(data', mu');
wolffd@0: %denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
wolffd@0: denom = (2*pi*Sigma)^(d/2); % sqrt(det(2*pi*Sigma))
wolffd@0: numer = exp(-0.5/Sigma  * D');
wolffd@0: B2 = numer / denom;
wolffd@0: B2 = reshape(B2, [Q M T]);
wolffd@0: 
wolffd@0: tic; B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2)); toc
wolffd@0: 
wolffd@0: tic
wolffd@0: A = zeros(Q,T);
wolffd@0: for q=1:Q
wolffd@0:   A(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % sum over m
wolffd@0: end
wolffd@0: toc
wolffd@0: assert(approxeq(A,B))
wolffd@0: 
wolffd@0: tic
wolffd@0: A = zeros(Q,T);
wolffd@0: for t=1:T
wolffd@0:   A(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m
wolffd@0: end
wolffd@0: toc
wolffd@0: assert(approxeq(A,B))
wolffd@0: 
wolffd@0:     
wolffd@0: 
wolffd@0: 
wolffd@0: mu = reshape(mu, [d Q M]);
wolffd@0: B3 = zeros(Q,M,T);
wolffd@0: for j=1:Q
wolffd@0:   for k=1:M
wolffd@0:     B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma*eye(d));
wolffd@0:   end
wolffd@0: end
wolffd@0: assert(approxeq(B2, B3))
wolffd@0: 
wolffd@0: logB4 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d
wolffd@0: B4 = reshape(exp(logB4), [Q M T]);
wolffd@0: assert(approxeq(B4, B3))
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0:   
wolffd@0: % tied cov matrix
wolffd@0: 
wolffd@0: Sigma = rand_psd(d,d);
wolffd@0: mu = reshape(mu, [d Q*M]);
wolffd@0: D = sqdist(data, mu, inv(Sigma))';
wolffd@0: denom = sqrt(det(2*pi*Sigma));
wolffd@0: numer = exp(-0.5 * D);
wolffd@0: B2 = numer / denom;
wolffd@0: B2 = reshape(B2, [Q M T]);
wolffd@0: 
wolffd@0: mu = reshape(mu, [d Q M]);
wolffd@0: B3 = zeros(Q,M,T);
wolffd@0: for j=1:Q
wolffd@0:   for k=1:M
wolffd@0:     B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma);
wolffd@0:   end
wolffd@0: end
wolffd@0: assert(approxeq(B2, B3))
wolffd@0: 
wolffd@0: logB4 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D;
wolffd@0: B4 = reshape(exp(logB4), [Q M T]);
wolffd@0: assert(approxeq(B4, B3))