Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMstats/mixgauss_prob_test.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 test_eval_pdf_cond_mixgauss()
wolffd@0 2
wolffd@0 3 %Q = 10; M = 100; d = 20; T = 500;
wolffd@0 4 Q = 2; M = 3; d = 4; T = 5;
wolffd@0 5
wolffd@0 6 mu = rand(d,Q,M);
wolffd@0 7 data = randn(d,T);
wolffd@0 8 %mixmat = mk_stochastic(rand(Q,M));
wolffd@0 9 mixmat = mk_stochastic(ones(Q,M));
wolffd@0 10
wolffd@0 11 % tied scalar
wolffd@0 12 Sigma = 0.01;
wolffd@0 13
wolffd@0 14 mu = rand(d,M,Q);
wolffd@0 15 weights = mixmat';
wolffd@0 16 N = M*ones(1,Q);
wolffd@0 17 tic; [B, B2, D] = parzen(data, mu, Sigma, N, weights); toc
wolffd@0 18 tic; [BC, B2C, DC] = parzenC(data, mu, Sigma, N); toc
wolffd@0 19 approxeq(B,BC)
wolffd@0 20 B2C = reshape(B2C,[M Q T]);
wolffd@0 21 approxeq(B2,B2C)
wolffd@0 22 DC = reshape(DC,[M Q T]);
wolffd@0 23 approxeq(D,DC)
wolffd@0 24
wolffd@0 25
wolffd@0 26 return
wolffd@0 27
wolffd@0 28 tic; [B, B2] = eval_pdf_cond_mixgauss(data, mu, Sigma, mixmat); toc
wolffd@0 29 tic; C = eval_pdf_cond_parzen(data, mu, Sigma); toc
wolffd@0 30 approxeq(B,C)
wolffd@0 31
wolffd@0 32 return;
wolffd@0 33
wolffd@0 34
wolffd@0 35 mu = reshape(mu, [d Q*M]);
wolffd@0 36
wolffd@0 37 data = mk_unit_norm(data);
wolffd@0 38 mu = mk_unit_norm(mu);
wolffd@0 39 tic; D = 2 -2*(data'*mu); toc % avoid an expensive repmat
wolffd@0 40 tic; D2 = sqdist(data, mu); toc
wolffd@0 41 approxeq(D,D2)
wolffd@0 42
wolffd@0 43
wolffd@0 44 % D(t,m) = sq dist between data(:,t) and mu(:,m)
wolffd@0 45 mu = reshape(mu, [d Q*M]);
wolffd@0 46 D = dist2(data', mu');
wolffd@0 47 %denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
wolffd@0 48 denom = (2*pi*Sigma)^(d/2); % sqrt(det(2*pi*Sigma))
wolffd@0 49 numer = exp(-0.5/Sigma * D');
wolffd@0 50 B2 = numer / denom;
wolffd@0 51 B2 = reshape(B2, [Q M T]);
wolffd@0 52
wolffd@0 53 tic; B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2)); toc
wolffd@0 54
wolffd@0 55 tic
wolffd@0 56 A = zeros(Q,T);
wolffd@0 57 for q=1:Q
wolffd@0 58 A(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % sum over m
wolffd@0 59 end
wolffd@0 60 toc
wolffd@0 61 assert(approxeq(A,B))
wolffd@0 62
wolffd@0 63 tic
wolffd@0 64 A = zeros(Q,T);
wolffd@0 65 for t=1:T
wolffd@0 66 A(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m
wolffd@0 67 end
wolffd@0 68 toc
wolffd@0 69 assert(approxeq(A,B))
wolffd@0 70
wolffd@0 71
wolffd@0 72
wolffd@0 73
wolffd@0 74 mu = reshape(mu, [d Q M]);
wolffd@0 75 B3 = zeros(Q,M,T);
wolffd@0 76 for j=1:Q
wolffd@0 77 for k=1:M
wolffd@0 78 B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma*eye(d));
wolffd@0 79 end
wolffd@0 80 end
wolffd@0 81 assert(approxeq(B2, B3))
wolffd@0 82
wolffd@0 83 logB4 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d
wolffd@0 84 B4 = reshape(exp(logB4), [Q M T]);
wolffd@0 85 assert(approxeq(B4, B3))
wolffd@0 86
wolffd@0 87
wolffd@0 88
wolffd@0 89
wolffd@0 90 % tied cov matrix
wolffd@0 91
wolffd@0 92 Sigma = rand_psd(d,d);
wolffd@0 93 mu = reshape(mu, [d Q*M]);
wolffd@0 94 D = sqdist(data, mu, inv(Sigma))';
wolffd@0 95 denom = sqrt(det(2*pi*Sigma));
wolffd@0 96 numer = exp(-0.5 * D);
wolffd@0 97 B2 = numer / denom;
wolffd@0 98 B2 = reshape(B2, [Q M T]);
wolffd@0 99
wolffd@0 100 mu = reshape(mu, [d Q M]);
wolffd@0 101 B3 = zeros(Q,M,T);
wolffd@0 102 for j=1:Q
wolffd@0 103 for k=1:M
wolffd@0 104 B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma);
wolffd@0 105 end
wolffd@0 106 end
wolffd@0 107 assert(approxeq(B2, B3))
wolffd@0 108
wolffd@0 109 logB4 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D;
wolffd@0 110 B4 = reshape(exp(logB4), [Q M T]);
wolffd@0 111 assert(approxeq(B4, B3))