Дискретный поверхностный интеграл с cumsum

У меня есть матрица z(x,y) Это абитарный PDF NxN, построенный из уникальная оценка плотности ядра (т.е. не обычный pdf и не имеет функции). Он многомерен и не может быть разделен и является дискретными данными.

Я не хочу строить матрицу NxN (F(x,y)), которая является кумулятивной функцией распределения в двух измерениях этого PDF-файла, чтобы затем я мог случайным образом выбрать F(x,y) = P(x ‹ X , у ‹ Y);

Аналитически я думаю, что CDF многомерной функции является поверхностным интегралом PDF.

То, что я пробовал, - это использовать функцию cumsum для вычисления поверхностного интеграла и протестировать это с многомерной нормалью по отношению к аналитическому решению, и, похоже, между ними есть некоторое несоответствие:

% multivariate parameters
delta = 100;
mu = [1 1];
Sigma = [0.25 .3; .3 1];
x1 = linspace(-2,4,delta); x2 = linspace(-2,4,delta);
[X1,X2] = meshgrid(x1,x2);
% Calculate Normal multivariate pdf
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
% My attempt at a numerical surface integral 
FN = cumsum(cumsum(F,1),2);
% Normalise the CDF
FN = FN./max(max(FN));
X = [X1(:) X2(:)];
% Analytic solution to a multivariate normal pdf
p = mvncdf(X,mu,Sigma);
p = reshape(p,delta,delta);
% Highlight the difference
dif = p - FN;
error = max(max(sqrt(dif.^2)));
% %% Plot
figure(1)
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
xlabel('x1'); ylabel('x2'); zlabel('Probability Density');
figure(2)
surf(X1,X2,FN);
xlabel('x1'); ylabel('x2');
figure(3);
surf(X1,X2,p);
xlabel('x1'); ylabel('x2');
figure(5)
surf(X1,X2,dif)
xlabel('x1'); ylabel('x2');

В частности, ошибка, по-видимому, находится в переходной области, которая является наиболее важной.

У кого-нибудь есть лучшее решение этой проблемы или я вижу, что я делаю неправильно?? Любая помощь приветствуется!

РЕДАКТИРОВАТЬ: Это желаемый результат кумулятивной интеграции. Причина, по которой эта функция представляет ценность для меня, заключается в том, что при случайном создании выборок из этой функции на закрытом интервале [0,1] более высокий взвешенный (т.е. более вероятные) значения появляются чаще, таким образом выборки сходятся к ожидаемому(ым) значению(ам) (в случае множественных пиков) это желаемый результат для таких алгоритмов, как фильтры частиц, нейронные сети и т.д.