|
wolffd@0
|
1 % Factor analysis
|
|
wolffd@0
|
2 % Z -> X, Z in R^k, X in R^D, k << D (high dimensional observations explained by small source)
|
|
wolffd@0
|
3 % Z ~ N(0,I), X|Z ~ N(L z, Psi), where Psi is diagonal.
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % Mixtures of FA
|
|
wolffd@0
|
6 % Now X|Z,W=i ~ N(mu(i) + L(i) Z, Psi(i))
|
|
wolffd@0
|
7 %
|
|
wolffd@0
|
8 % We compare to Zoubin Ghahramani's code.
|
|
wolffd@0
|
9
|
|
wolffd@0
|
10 randn('state', 0);
|
|
wolffd@0
|
11 max_iter = 3;
|
|
wolffd@0
|
12 M = 2;
|
|
wolffd@0
|
13 k = 3;
|
|
wolffd@0
|
14 D = 5;
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 n = 5;
|
|
wolffd@0
|
17 X1 = randn(n, D);
|
|
wolffd@0
|
18 X2 = randn(n, D) + 2; % move the mean to (2,2,2...)
|
|
wolffd@0
|
19 X = [X1; X2];
|
|
wolffd@0
|
20 N = size(X, 1);
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22 % initialise as in mfa
|
|
wolffd@0
|
23 tiny=exp(-700);
|
|
wolffd@0
|
24 mX = mean(X);
|
|
wolffd@0
|
25 cX=cov(X);
|
|
wolffd@0
|
26 scale=det(cX)^(1/D);
|
|
wolffd@0
|
27 randn('state',0); % must reset seed here so initial params are identical to mfa
|
|
wolffd@0
|
28 L0=randn(D*M,k)*sqrt(scale/k);
|
|
wolffd@0
|
29 W0 = permute(reshape(L0, [D M k]), [1 3 2]); % use D,K,M
|
|
wolffd@0
|
30 Psi0=diag(cX)+tiny;
|
|
wolffd@0
|
31 Pi0=ones(M,1)/M;
|
|
wolffd@0
|
32 Mu0=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 [Lh1, Ph1, Mu1, Pi1, LL1] = mfa(X,M,k,max_iter);
|
|
wolffd@0
|
35 Lh1 = permute(reshape(Lh1, [D M k]), [1 3 2]); % use D,K,M
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 ns = [M k D];
|
|
wolffd@0
|
39 dag = zeros(3);
|
|
wolffd@0
|
40 dag(1,3) = 1;
|
|
wolffd@0
|
41 dag(2,3) = 1;
|
|
wolffd@0
|
42 dnodes = 1;
|
|
wolffd@0
|
43 onodes = 3;
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'observed', onodes);
|
|
wolffd@0
|
46 bnet.CPD{1} = tabular_CPD(bnet, 1, Pi0);
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48 %bnet.CPD{2} = gaussian_CPD(bnet, 2, zeros(k, 1), eye(k), [], 'diag', 'untied', 'clamp_mean', 'clamp_cov');
|
|
wolffd@0
|
49
|
|
wolffd@0
|
50 bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(k, 1), 'cov', eye(k), 'cov_type', 'diag', ...
|
|
wolffd@0
|
51 'cov_prior_weight', 0, 'clamp_mean', 1, 'clamp_cov', 1);
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 %bnet.CPD{3} = gaussian_CPD(bnet, 3, Mu0', repmat(diag(Psi0), [1 1 M]), W0, 'diag', 'tied');
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', Mu0', 'cov', repmat(diag(Psi0), [1 1 M]), ...
|
|
wolffd@0
|
56 'weights', W0, 'cov_type', 'diag', 'cov_prior_weight', 0, 'tied_cov', 1);
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58 engine = jtree_inf_engine(bnet);
|
|
wolffd@0
|
59 evidence = cell(3, N);
|
|
wolffd@0
|
60 evidence(3,:) = num2cell(X', 1);
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 [bnet2, LL2, engine2] = learn_params_em(engine, evidence, max_iter);
|
|
wolffd@0
|
63
|
|
wolffd@0
|
64 s = struct(bnet2.CPD{1});
|
|
wolffd@0
|
65 Pi2 = s.CPT(:);
|
|
wolffd@0
|
66 s = struct(bnet2.CPD{3});
|
|
wolffd@0
|
67 Mu2 = s.mean;
|
|
wolffd@0
|
68 W2 = s.weights;
|
|
wolffd@0
|
69 Sigma2 = s.cov;
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 % Compare to Zoubin's code
|
|
wolffd@0
|
73 assert(approxeq(LL1,LL2));
|
|
wolffd@0
|
74 for i=1:M
|
|
wolffd@0
|
75 assert(approxeq(W2(:,:,i), Lh1(:,:,i)));
|
|
wolffd@0
|
76 assert(approxeq(Sigma2(:,:,i), diag(Ph1)));
|
|
wolffd@0
|
77 assert(approxeq(Mu2(:,i), Mu1(i,:)));
|
|
wolffd@0
|
78 assert(approxeq(Pi2(:), Pi1(:)));
|
|
wolffd@0
|
79 end
|
|
wolffd@0
|
80
|
