annotate toolboxes/FullBNT-1.0.7/KPMstats/mixgauss_em.m @ 0:cc4b1211e677 tip

initial commit to HG from Changeset: 646 (e263d8a21543) added further path and more save "camirversion.m"
author Daniel Wolff
date Fri, 19 Aug 2016 13:07:06 +0200
parents
children
rev   line source
Daniel@0 1 function [mu, Sigma, prior] = mixgauss_em(Y, nc, varargin)
Daniel@0 2 % MIXGAUSS_EM Fit the parameters of a mixture of Gaussians using EM
Daniel@0 3 % function [mu, Sigma, prior] = mixgauss_em(data, nc, varargin)
Daniel@0 4 %
Daniel@0 5 % data(:, t) is the t'th data point
Daniel@0 6 % nc is the number of clusters
Daniel@0 7
Daniel@0 8 % Kevin Murphy, 13 May 2003
Daniel@0 9
Daniel@0 10 [max_iter, thresh, cov_type, mu, Sigma, method, ...
Daniel@0 11 cov_prior, verbose, prune_thresh] = process_options(...
Daniel@0 12 varargin, 'max_iter', 10, 'thresh', 1e-2, 'cov_type', 'full', ...
Daniel@0 13 'mu', [], 'Sigma', [], 'method', 'kmeans', ...
Daniel@0 14 'cov_prior', [], 'verbose', 0, 'prune_thresh', 0);
Daniel@0 15
Daniel@0 16 [ny T] = size(Y);
Daniel@0 17
Daniel@0 18 if nc==1
Daniel@0 19 % No latent variable, so there is a closed-form solution
Daniel@0 20 mu = mean(Y')';
Daniel@0 21 Sigma = cov(Y');
Daniel@0 22 if strcmp(cov_type, 'diag')
Daniel@0 23 Sigma = diag(diag(Sigma));
Daniel@0 24 end
Daniel@0 25 prior = 1;
Daniel@0 26 return;
Daniel@0 27 end
Daniel@0 28
Daniel@0 29 if isempty(mu)
Daniel@0 30 [mu, Sigma, prior] = mixgauss_init(nc, Y, cov_type, method);
Daniel@0 31 end
Daniel@0 32
Daniel@0 33 previous_loglik = -inf;
Daniel@0 34 num_iter = 1;
Daniel@0 35 converged = 0;
Daniel@0 36
Daniel@0 37 %if verbose, fprintf('starting em\n'); end
Daniel@0 38
Daniel@0 39 while (num_iter <= max_iter) & ~converged
Daniel@0 40 % E step
Daniel@0 41 probY = mixgauss_prob(Y, mu, Sigma, prior); % probY(q,t)
Daniel@0 42 [post, lik] = normalize(probY .* repmat(prior, 1, T), 1); % post(q,t)
Daniel@0 43 loglik = log(sum(lik));
Daniel@0 44
Daniel@0 45 % extract expected sufficient statistics
Daniel@0 46 w = sum(post,2); % w(c) = sum_t post(c,t)
Daniel@0 47 WYY = zeros(ny, ny, nc); % WYY(:,:,c) = sum_t post(c,t) Y(:,t) Y(:,t)'
Daniel@0 48 WY = zeros(ny, nc); % WY(:,c) = sum_t post(c,t) Y(:,t)
Daniel@0 49 WYTY = zeros(nc,1); % WYTY(c) = sum_t post(c,t) Y(:,t)' Y(:,t)
Daniel@0 50 for c=1:nc
Daniel@0 51 weights = repmat(post(c,:), ny, 1); % weights(:,t) = post(c,t)
Daniel@0 52 WYbig = Y .* weights; % WYbig(:,t) = post(c,t) * Y(:,t)
Daniel@0 53 WYY(:,:,c) = WYbig * Y';
Daniel@0 54 WY(:,c) = sum(WYbig, 2);
Daniel@0 55 WYTY(c) = sum(diag(WYbig' * Y));
Daniel@0 56 end
Daniel@0 57
Daniel@0 58 % M step
Daniel@0 59 prior = normalize(w);
Daniel@0 60 [mu, Sigma] = mixgauss_Mstep(w, WY, WYY, WYTY, 'cov_type', cov_type, 'cov_prior', cov_prior);
Daniel@0 61
Daniel@0 62 if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
Daniel@0 63 num_iter = num_iter + 1;
Daniel@0 64 converged = em_converged(loglik, previous_loglik, thresh);
Daniel@0 65 previous_loglik = loglik;
Daniel@0 66
Daniel@0 67 end
Daniel@0 68
Daniel@0 69 if prune_thresh > 0
Daniel@0 70 ndx = find(prior < prune_thresh);
Daniel@0 71 mu(:,ndx) = [];
Daniel@0 72 Sigma(:,:,ndx) = [];
Daniel@0 73 prior(ndx) = [];
Daniel@0 74 end