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