wolffd@0: % Fit a mixture of Gaussians using netlab and BNT
wolffd@0: 
wolffd@0: rand('state', 0);
wolffd@0: randn('state', 0);
wolffd@0: 
wolffd@0: % Q -> Y
wolffd@0: ncenters = 2; dim = 2;
wolffd@0: cov_type = 'full';
wolffd@0: 
wolffd@0: % Generate the data from a mixture of 2 Gaussians
wolffd@0: %mu = randn(dim, ncenters);
wolffd@0: mu = zeros(dim, ncenters);
wolffd@0: mu(:,1) = [-1 -1]';
wolffd@0: mu(:,1) = [1 1]';
wolffd@0: Sigma = repmat(0.1*eye(dim),[1 1 ncenters]);
wolffd@0: ndat1 = 8; ndat2 = 8;
wolffd@0: %ndat1 = 2; ndat2 = 2;
wolffd@0: ndata = ndat1+ndat2;
wolffd@0: x1 = gsamp(mu(:,1), Sigma(:,:,1), ndat1);
wolffd@0: x2 = gsamp(mu(:,2), Sigma(:,:,2), ndat2);
wolffd@0: data = [x1; x2];
wolffd@0: %plot(x1(:,1),x1(:,2),'ro', x2(:,1),x2(:,2),'bx')
wolffd@0: 
wolffd@0: % Fit using netlab
wolffd@0: max_iter = 3;
wolffd@0: mix = gmm(dim, ncenters, cov_type);
wolffd@0: options = foptions;
wolffd@0: options(1) = 1; % verbose
wolffd@0: options(14) = max_iter;
wolffd@0: 
wolffd@0: % extract initial params
wolffd@0: %mix = gmminit(mix, x, options); % Initialize with K-means
wolffd@0: mu0 = mix.centres';
wolffd@0: pi0 = mix.priors(:);
wolffd@0: Sigma0 = mix.covars; % repmat(eye(dim), [1 1 ncenters]);
wolffd@0: 
wolffd@0: [mix, options] = gmmem(mix, data, options);
wolffd@0: 
wolffd@0: % Final params
wolffd@0: ll1 = options(8);
wolffd@0: mu1 = mix.centres';
wolffd@0: pi1 = mix.priors(:);
wolffd@0: Sigma1 = mix.covars;
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: % BNT
wolffd@0: 
wolffd@0: dag = zeros(2);
wolffd@0: dag(1,2) = 1;
wolffd@0: node_sizes = [ncenters dim];
wolffd@0: discrete_nodes = 1;
wolffd@0: onodes = 2;
wolffd@0: 
wolffd@0: bnet = mk_bnet(dag, node_sizes, 'discrete', discrete_nodes, 'observed', onodes);
wolffd@0: bnet.CPD{1} = tabular_CPD(bnet, 1, pi0);
wolffd@0: bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu0, 'cov', Sigma0, 'cov_type', cov_type, ...
wolffd@0: 			   'cov_prior_weight', 0);
wolffd@0: 
wolffd@0: engine = jtree_inf_engine(bnet);
wolffd@0: 
wolffd@0: evidence = cell(2, ndata);
wolffd@0: evidence(2,:) = num2cell(data', 1);
wolffd@0: 
wolffd@0: [bnet2, LL] = learn_params_em(engine, evidence, max_iter);
wolffd@0: 
wolffd@0: ll2 = LL(end);
wolffd@0: s1 = struct(bnet2.CPD{1});
wolffd@0: pi2 = s1.CPT(:);
wolffd@0: 
wolffd@0: s2 = struct(bnet2.CPD{2});
wolffd@0: mu2 = s2.mean;
wolffd@0: Sigma2 = s2.cov;
wolffd@0: 
wolffd@0: % assert(approxeq(ll1, ll2)); % gmmem returns the value after the final M step, GMT before
wolffd@0: assert(approxeq(mu1, mu2));
wolffd@0: assert(approxeq(Sigma1, Sigma2))
wolffd@0: assert(approxeq(pi1, pi2))
wolffd@0: 
wolffd@0: