wolffd@0: function [new_mu, new_Sigma, new_Sigma2] = collapse_mog(mu, Sigma, coefs)
wolffd@0: % COLLAPSE_MOG Collapse a mixture of Gaussians to a single Gaussian by moment matching
wolffd@0: % [new_mu, new_Sigma] = collapse_mog(mu, Sigma, coefs)
wolffd@0: %
wolffd@0: % coefs(i) - weight of i'th mixture component
wolffd@0: % mu(:,i), Sigma(:,:,i) - params of i'th mixture component
wolffd@0: 
wolffd@0: % S = sum_c w_c (S_c + m_c m_c' + m m' - 2 m_c m') 
wolffd@0: %   = sum_c w_c (S_c + m_c m_c') + m m' - 2 (sum_c m_c) m'
wolffd@0: %   = sum_c w_c (S_c + m_c m_c') - m m'
wolffd@0: 
wolffd@0: new_mu = sum(mu * diag(coefs), 2); % weighted sum of columns
wolffd@0: 
wolffd@0: n = length(new_mu);
wolffd@0: new_Sigma = zeros(n,n);
wolffd@0: new_Sigma2 = zeros(n,n);
wolffd@0: for j=1:length(coefs)
wolffd@0:   m = mu(:,j) - new_mu;
wolffd@0:   new_Sigma = new_Sigma + coefs(j) * (Sigma(:,:,j) + m*m');
wolffd@0:   new_Sigma2 = new_Sigma2 + coefs(j) * (Sigma(:,:,j) + mu(:,j)*mu(:,j)');
wolffd@0: end
wolffd@0: %assert(approxeq(new_Sigma, new_Sigma2 - new_mu*new_mu'))