Daniel@0: function p = gaussian_prob(x, m, C, use_log)
Daniel@0: % GAUSSIAN_PROB Evaluate a multivariate Gaussian density.
Daniel@0: % p = gaussian_prob(X, m, C)
Daniel@0: % p(i) = N(X(:,i), m, C) where C = covariance matrix and each COLUMN of x is a datavector
Daniel@0: 
Daniel@0: % p = gaussian_prob(X, m, C, 1) returns log N(X(:,i), m, C) (to prevents underflow).
Daniel@0: %
Daniel@0: % If X has size dxN, then p has size Nx1, where N = number of examples
Daniel@0: 
Daniel@0: if nargin < 4, use_log = 0; end
Daniel@0: 
Daniel@0: if length(m)==1 % scalar
Daniel@0:   x = x(:)';
Daniel@0: end
Daniel@0: [d N] = size(x);
Daniel@0: %assert(length(m)==d); % slow
Daniel@0: m = m(:);
Daniel@0: M = m*ones(1,N); % replicate the mean across columns
Daniel@0: denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
Daniel@0: mahal = sum(((x-M)'*inv(C)).*(x-M)',2);   % Chris Bregler's trick
Daniel@0: if any(mahal<0)
Daniel@0:   warning('mahal < 0 => C is not psd')
Daniel@0: end
Daniel@0: if use_log
Daniel@0:   p = -0.5*mahal - log(denom);
Daniel@0: else
Daniel@0:   p = exp(-0.5*mahal) / (denom+eps);
Daniel@0: end