wolffd@0: function [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin)
wolffd@0: % MSTEP_COND_GAUSS Compute MLEs for mixture of Gaussians given expected sufficient statistics
wolffd@0: % function [mu, Sigma] = Mstep_cond_gauss(w, Y, YY, YTY, varargin)
wolffd@0: %
wolffd@0: % We assume P(Y|Q=i) = N(Y; mu_i, Sigma_i)
wolffd@0: % and w(i,t) = p(Q(t)=i|y(t)) = posterior responsibility
wolffd@0: % See www.ai.mit.edu/~murphyk/Papers/learncg.pdf.
wolffd@0: %
wolffd@0: % INPUTS:
wolffd@0: % w(i) = sum_t w(i,t) = responsibilities for each mixture component
wolffd@0: %  If there is only one mixture component (i.e., Q does not exist),
wolffd@0: %  then w(i) = N = nsamples,  and 
wolffd@0: %  all references to i can be replaced by 1.
wolffd@0: % YY(:,:,i) = sum_t w(i,t) y(:,t) y(:,t)' = weighted outer product
wolffd@0: % Y(:,i) = sum_t w(i,t) y(:,t) = weighted observations
wolffd@0: % YTY(i) = sum_t w(i,t) y(:,t)' y(:,t) = weighted inner product
wolffd@0: %   You only need to pass in YTY if Sigma is to be estimated as spherical.
wolffd@0: %
wolffd@0: % Optional parameters may be passed as 'param_name', param_value pairs.
wolffd@0: % Parameter names are shown below; default values in [] - if none, argument is mandatory.
wolffd@0: %
wolffd@0: % 'cov_type' - 'full', 'diag' or 'spherical' ['full']
wolffd@0: % 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0]
wolffd@0: % 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ]
wolffd@0: % 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ]
wolffd@0: % 'cov_prior' - Lambda_i, added to YY(:,:,i) [0.01*eye(d,d,Q)]
wolffd@0: %
wolffd@0: % If covariance is tied, Sigma has size d*d.
wolffd@0: % But diagonal and spherical covariances are represented in full size.
wolffd@0: 
wolffd@0: [cov_type, tied_cov,  clamped_cov, clamped_mean, cov_prior, other] = ...
wolffd@0:     process_options(varargin,...
wolffd@0: 		    'cov_type', 'full', 'tied_cov', 0,  'clamped_cov', [], 'clamped_mean', [], ...
wolffd@0: 		    'cov_prior', []);
wolffd@0: 
wolffd@0: [Ysz Q] = size(Y);
wolffd@0: N = sum(w);
wolffd@0: if isempty(cov_prior)
wolffd@0:   %cov_prior = zeros(Ysz, Ysz, Q);
wolffd@0:   %for q=1:Q
wolffd@0:   %  cov_prior(:,:,q) = 0.01*cov(Y(:,q)');
wolffd@0:   %end
wolffd@0:   cov_prior = repmat(0.01*eye(Ysz,Ysz), [1 1 Q]);
wolffd@0: end
wolffd@0: %YY = reshape(YY, [Ysz Ysz Q]) + cov_prior; % regularize the scatter matrix
wolffd@0: YY = reshape(YY, [Ysz Ysz Q]);
wolffd@0: 
wolffd@0: % Set any zero weights to one before dividing
wolffd@0: % This is valid because w(i)=0 => Y(:,i)=0, etc
wolffd@0: w = w + (w==0);
wolffd@0: 		    
wolffd@0: if ~isempty(clamped_mean)
wolffd@0:   mu = clamped_mean;
wolffd@0: else
wolffd@0:   % eqn 6
wolffd@0:   %mu = Y ./ repmat(w(:)', [Ysz 1]);% Y may have a funny size
wolffd@0:   mu = zeros(Ysz, Q);
wolffd@0:   for i=1:Q
wolffd@0:     mu(:,i) = Y(:,i) / w(i);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: if ~isempty(clamped_cov)
wolffd@0:   Sigma = clamped_cov;
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: if ~tied_cov
wolffd@0:   Sigma = zeros(Ysz,Ysz,Q);
wolffd@0:   for i=1:Q
wolffd@0:     if cov_type(1) == 's'
wolffd@0:       % eqn 17
wolffd@0:       s2 = (1/Ysz)*( (YTY(i)/w(i)) - mu(:,i)'*mu(:,i) );
wolffd@0:       Sigma(:,:,i) = s2 * eye(Ysz);
wolffd@0:     else
wolffd@0:       % eqn 12
wolffd@0:       SS = YY(:,:,i)/w(i)  - mu(:,i)*mu(:,i)';
wolffd@0:       if cov_type(1)=='d'
wolffd@0: 	SS = diag(diag(SS));
wolffd@0:       end
wolffd@0:       Sigma(:,:,i) = SS;
wolffd@0:     end
wolffd@0:   end
wolffd@0: else % tied cov
wolffd@0:   if cov_type(1) == 's'
wolffd@0:     % eqn 19
wolffd@0:     s2 = (1/(N*Ysz))*(sum(YTY,2) + sum(diag(mu'*mu) .* w));
wolffd@0:     Sigma = s2*eye(Ysz);
wolffd@0:   else
wolffd@0:     SS = zeros(Ysz, Ysz);
wolffd@0:     % eqn 15
wolffd@0:     for i=1:Q % probably could vectorize this...
wolffd@0:       SS = SS + YY(:,:,i)/N - mu(:,i)*mu(:,i)';
wolffd@0:     end
wolffd@0:     if cov_type(1) == 'd'
wolffd@0:       Sigma = diag(diag(SS));
wolffd@0:     else
wolffd@0:       Sigma = SS;
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: if tied_cov
wolffd@0:   Sigma =  repmat(Sigma, [1 1 Q]);
wolffd@0: end
wolffd@0: Sigma = Sigma + cov_prior;