wolffd@0: % Factor analysis
wolffd@0: % Z -> X,  Z in R^k, X in R^D, k << D (high dimensional observations explained by small source)
wolffd@0: % Z ~ N(0,I),   X|Z ~ N(L z, Psi), where Psi is diagonal.
wolffd@0: %
wolffd@0: % We compare to Zoubin Ghahramani's code.
wolffd@0: 
wolffd@0: state = 0;
wolffd@0: rand('seed', state);
wolffd@0: randn('seed', state);
wolffd@0: max_iter = 3;
wolffd@0: k = 2;
wolffd@0: D = 4;
wolffd@0: N = 10;
wolffd@0: X = randn(N, D);
wolffd@0: 
wolffd@0: % Initialize as in Zoubin's ffa (fast factor analysis)
wolffd@0: X=X-ones(N,1)*mean(X);
wolffd@0: XX=X'*X/N;
wolffd@0: diagXX=diag(XX);
wolffd@0: cX=cov(X);
wolffd@0: scale=det(cX)^(1/D);
wolffd@0: randn('seed', 0);  % must reset seed here so initial params are identical to mfa
wolffd@0: L0=randn(D,k)*sqrt(scale/k);
wolffd@0: W0 = L0;
wolffd@0: Psi0=diag(cX);
wolffd@0: 
wolffd@0: [L1, Psi1, LL1] = ffa(X,k,max_iter);
wolffd@0: 
wolffd@0: 
wolffd@0: ns = [k D];
wolffd@0: dag = zeros(2,2);
wolffd@0: dag(1,2) = 1;
wolffd@0: bnet = mk_bnet(dag, ns, 'discrete', [], 'observed', 2);
wolffd@0: bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', zeros(k,1), 'cov', eye(k), 'cov_type', 'diag', ...
wolffd@0: 			   'clamp_mean', 1, 'clamp_cov', 1);
wolffd@0: bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(D,1), 'cov', diag(Psi0), 'weights', W0, ...
wolffd@0: 			   'cov_type', 'diag', 'cov_prior_weight', 0, 'clamp_mean', 1);
wolffd@0: 
wolffd@0: engine = jtree_inf_engine(bnet);
wolffd@0: evidence = cell(2,N);
wolffd@0: evidence(2,:) = num2cell(X', 1);
wolffd@0: 
wolffd@0: [bnet2, LL2] = learn_params_em(engine, evidence, max_iter);
wolffd@0: 
wolffd@0: s = struct(bnet2.CPD{2});
wolffd@0: L2 = s.weights;
wolffd@0: Psi2 = s.cov;
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: % Compare to Zoubin's code
wolffd@0: assert(approxeq(LL2, LL1));
wolffd@0: assert(approxeq(Psi2, diag(Psi1)));
wolffd@0: assert(approxeq(L2, L1));
wolffd@0: 
wolffd@0: 
wolffd@0: