In MATLAB, I've written two snippets of code that compute the PDF of a multivariate normal distribution. However there's a difference in the values these two methods produce and I can't figure out why. I've narrowed the problem down to something having to do with computing the inverse of the covariance matrix.
function p = mvnpdf_inacc(X, mu, sigma) xc = bsxfun(@minus, X, mu); [n, k] = size(xc); twopic = (2 * pi) ^ (-k / 2); sqrtdetsig = sqrt(det(sigma)) ^ -1; c = twopic * sqrtdetsig; p = zeros(n, 1); for i = 1:n xci = xc(i, :); p(i) = c * exp(-0.5 * (xci / sigma * xci')); end end
function p = mvnpdf_acc(X, mu, sigma) [R, err] = cholcov(sigma, 0); if err error('%s', 'sigma is not both symmetric and positive definite'); end X0 = bsxfun(@minus, X, mu) / R; d = min(size(X)); slogdet = sum(log(diag(R))); p = exp(-0.5 * sum(X0 .^ 2, 2) - slogdet - 0.5 * d * log(2 * pi)); end
function iseq_func = test_mvnpdf(n) x = linspace(-2, 2, n); y = x; [X, Y] = meshgrid(x, y); XY = [X(:), Y(:)]; mu = [0, 0]; sigma = [1.0, 0.5; 0.5, 0.4]; p_inacc = mvnpdf_inacc(XY, mu, sigma); p_acc = mvnpdf_acc(XY, mu, sigma); p_diff = abs(p_inacc - p_acc); iseq_func = nnz(p_diff) == 0; end
I get a value of false from running
iseq_func(25). What the heck is going on here? Thanks!