c++ - 如何在密度函数的特定值处有效地绘制二维正态分布的等值线？

Question

我有一个二维正态分布，由其均值和协方差矩阵表示。现在我想在密度函数
所在的每个点周围画一条线
超过一定的阈值。 (省略了归一化项，因此阈值可以应用于所有分布，而不管它们的大小。)

当然，这可以通过遍历整个图像并检查每个像素的公式来完成(参见 drawVariant2())。然而，这非常慢，我希望绘图机制或多或少具有实时能力。

我发现的另一种方法计算协方差矩阵的特征值和 - vector ，并使用它们来变换圆的图像(参见 drawVariant1())。这样做的问题是，我没有表示绘制的线与平均值的距离在密度函数方面意味着什么。

有什么方法可以在密度函数的某些值处绘制我的分布等值线？

这是我目前掌握的两种绘制方法的代码。它使用特征库进行矩阵和 vector 运算。

#include "Eigen/Core"
#include "Eigen/Eigenvalues"

unsigned int width = 500, height = 500;
Eigen::Matrix<double, 2, 1> mean;
Eigen::Matrix<double, 2, 2> cov;

const double C_PI = 3.14159265358979323846;

void drawVariant1(unsigned char * img)
{
    const int isolineRadius = 3;
    const double baseCircleSteps = 100;
    const double circleLength = 2 * C_PI * isolineRadius * baseCircleSteps;

    Eigen::EigenSolver<Eigen::Matrix<double, 2, 2>> eigenSolver(cov);

    eigenSolver.compute(cov);

    const double eVal1 = eigenSolver.eigenvalues().real()(0);
    const double eVal2 = eigenSolver.eigenvalues().real()(1);

    const Vector2 eVec1 = eigenSolver.eigenvectors().real().col(0);
    const Vector2 eVec2 = eigenSolver.eigenvectors().real().col(1);

    for (double phi = 0; phi < 2 * C_PI; phi += 2 * C_PI / circleLength)
    {
        const double x = isolineRadius * std::cos(phi);
        const double y = isolineRadius * std::sin(phi);

        const Vector2 posProjected = x * std::sqrt(eVal1)  * eVec1  + y * std::sqrt(eVal2) * eVec2 + mean;

        const int xP = (int)posProjected(0);
        const int yP = (int)posProjected(1);

        if (xP >= 0 && xP < width && yP >= 0 && yP < height)
        {
            img[yP * width + xP] = 255;
        }
    }
}


void drawVariant2(unsigned char * img)
{
    const double threshold = 0.5;
    for(int x = 0; x < width; ++x)
    {
        for(int y = 0; y < height; ++y)
        {
            const Eigen::Matrix<double, 2, 1> point((double)x, (double)y);
            const double likelihood = std::exp( -0.5 * (point - mean).transpose() * cov.inverse() * (point - mean) );

            if(likelihood >= threshold)
                img[y * width + x] = 128;
        }
    }
}


void main()
{
    mean << 100, 100;
    cov << 50, 0, 0, 20;

    unsigned char * img = new unsigned char[width * height];
    memset(img, 0, width*height);

    drawVariant1(img);
    drawVariant2(img);

    writePGM("test.pgm", img, width, height);

    delete [] img;
};

很抱歉，代码比较长，但我希望提供一个完整的示例可能有助于回答这个问题。

Answer 1

嗯，经过今天的投资，我想我终于找到了合适的解决方案。

想法是计算 x 轴上具有所需似然 w.r.t 的点。到标准正态分布。这可以通过二分法来完成。然后我可以用这个点的距离作为圆的半径，通过均值和协方差矩阵转换。

完整的源代码在这里。如果有人遇到类似的问题，我希望这会有所帮助。如果有人发现一些错误，我很乐意纠正它们。

#include "Eigen/Core"
#include "Eigen/Eigenvalues"

#include "PgmIO.h"

typedef Eigen::Matrix<double, 2, 1> Vector;
typedef Eigen::Matrix<double, 2, 2> Matrix;

const double C_PI = 3.14159265358979323846;

Vector mean;
Matrix cov;

int width = 500, height = 500;

void draw(unsigned char * img, double threshold)
{
    const double epsilon = 0.001;

    Matrix unity;
    Vector base;

    base << 1, 0;
    unity << 1, 0, 0, 1;

    double drawingDistance = 0;

    // Search for the Point on the X-Axis that has the desired likelihood with a bisection-method
    for(double a = 0, b = 5, likelihood = 0; std::abs(likelihood - threshold) > epsilon; )
    {
        const double currentDistance = (a + b) / 2;
        const Vector evalPoint = currentDistance * base;

        likelihood = std::exp(-0.5 * evalPoint.transpose() * unity * evalPoint); //unitiy.inverse() == unity

        // Suitable point found
        if(std::abs(likelihood - threshold) < epsilon)
        {
            drawingDistance = currentDistance;
            break;
        }

        if(likelihood > threshold)
        {
            a = currentDistance; // If the likelihood is too large search further away from the origin
        }
        else
        {
            b = currentDistance; // If the likelihood is too small search closer to the origin
        }
    }

    Eigen::EigenSolver<Matrix> eigenSolver(unity);

    eigenSolver.compute(cov);

    const double eVal1 = eigenSolver.eigenvalues().real()(0);
    const double eVal2 = eigenSolver.eigenvalues().real()(1);

    const Vector eVec1 = eigenSolver.eigenvectors().real().col(0);
    const Vector eVec2 = eigenSolver.eigenvectors().real().col(1);

    const double baseCircleSteps = 100;
    const double circleLength = 2 * C_PI * drawingDistance * baseCircleSteps;

    for (double phi = 0; phi < 2 * C_PI; phi += 2 * C_PI / circleLength)
    {
        // Compute the points on a circle within drawingDistance range
        const double x = drawingDistance * std::cos(phi);
        const double y = drawingDistance * std::sin(phi);

        // Project point to the eqivalent point on the isoline
        const Vector posProjected = x * std::sqrt(eVal1)  * eVec1  + y * std::sqrt(eVal2) * eVec2 + mean;

        const int xP = (int)posProjected(0);
        const int yP = (int)posProjected(1);

        // Set point in the image
        if (xP >= 0 && xP < width && yP >= 0 && yP < height)
        {
            img[yP * width + xP] = 255;
        }
    }
}


void main()
{
    mean << 100, 100;
    cov << 800, 100, 100, 500;

    unsigned char * img = new unsigned char[width * height];
    memset(img, 0, width*height);

    draw(img, 0.01);

    writePGM("img.pgm", img, width, height);

    delete [] img;
};