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