Как рассчитать PDF двумерного нормального распределения?

У меня есть набор точек и я извлекаю их небольшое подмножество для вычисления двумерного нормального распределения. После этого я проверяю все остальные точки, соответствуют ли они этому распределению, вычисляя PDF для каждой точки и отклоняя точки со значением ниже некоторого порогового значения.

Так много о теории …

PDF имеет согласно википедии формулу:

σ — стандартное отклонение, а µ — среднее значение, рассчитанное следующим образом:

cv::Scalar mean;
cv::Scalar stdDev;
dataPoints = dataPoints.reshape(3); // convert 3 columns to 3 channels
cv::meanStdDev(dataPoints, mean, stdDev);
dataPoints = dataPoints.reshape(1); // convert back
meanX = mean.val[0];
meanY = mean.val[1];
sigmaX = stdDev.val[0];
sigmaY = stdDev.val[1];

dataPoints — это cv :: Mat с 3 столбцами с плавающей точкой (x, y, index).

ρ — коэффициент корреляции, который я вычисляю так:

cv::matchTemplate(dataPoints.col(0), dataPoints.col(1), rho, cv::TM_CCOEFF_NORMED);

Последний шаг вычисляет вероятность для каждой точки, используя это:

double p = (1. / (2. * M_PI * sigmaX * sigmaY * sqrt(1. - pow(rho, 2))));
double e = exp((-1. / 2.) * D(x, y, rho));
double ret = p * e;

И D () должно быть, насколько я знаю, расстояние Махаланобиса, но формула из OpenCV cv::Mahalanobis(x, y, rho) возвращает другое значение, чем когда я сам его вычисляю:

double cX = (x - meanX) / sigmaX;
double cY = (y - meanY) / sigmaY;
double a = (1. / (1. - pow(rho, 2)));
double b = (pow(cX, 2) + pow(cY, 2) - 2. * rho * cX * cY);
double ret = a * b;

Так и сейчас моя проблема:

Насколько я знаю, интеграл по PDF должен быть 1, а максимальное значение PDF должно быть в (meanX, meanY), поэтому, когда σ будет 0, PDF в среднем должен быть 1. Но с помощью приведенных выше вычислений я могу получить значения, превышающие 1. Что я ошибаюсь?