Precisión numérica de la distribución normal multivariante

Question

En MATLAB, he escrito dos fragmentos de código que calcular el PDF de una distribución normal multivariante. Sin embargo, hay una diferencia en los valores de estos dos métodos producen y no puedo entender por qué. He reducido el problema a algo que tiene que ver con la computación de la inversa de la matriz de covarianza.

Inexacta código

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

Precisa de código

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

Código de prueba

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

Puedo obtener un valor de false de ejecución iseq_func(25). ¿Qué diablos está pasando aquí? Gracias!

Answer 1

Supongo que probó el código mediante el uso de sigma que no es positivo semidefinite. El 'precisa' de implementación (que supongo que es a partir de las estadísticas de la caja de herramientas) es computing $y^{\top}y$ donde $y = C^{-\top} \left(x - \mu\right)$ donde $C$ es el resultado de cholcov $\Sigma$ (sigma). Tenga en cuenta que $y^{\top}y$ debe ser no negativo por diseño. Si usted alimenta a en una $\Sigma$ que no está de PD, cholcov silenciosamente devuelve la versión de $C$ y el resto del código de producto (me han hecho de este un error, creo).

El código que escribió, sin embargo, es computing $\left(x - \mu\right)^{\top}\Sigma^{-1}\left(x - \mu\right)$. Si $\Sigma$ no es PD, esta cantidad puede ser negativo.

Al $\Sigma$ es positiva definida, cholcov devuelve un apropiado de la factorización de Cholesky, y los resultados son los mismos (hasta redondear).