机器算法验证 - 多元正态分布的数值精度 - 吾爱随笔录

在 MATLAB 中，我编写了两段代码来计算多元正态分布的 PDF。但是，这两种方法产生的值存在差异，我不知道为什么。我已经将问题缩小到与计算协方差矩阵的逆有关的问题。

不准确的代码

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

我从 running 得到一个 false 值iseq_func(25)。这到底是怎么回事？谢谢！