与 cumsum 的离散曲面积分答案

【问题标题】：Discrete surface integral with cumsum与 cumsum 的离散曲面积分
【发布时间】：2014-06-05 22:37:54
【问题描述】：

我有一个矩阵 z(x,y) 这是一个 NxN abitary pdf，由独特的内核密度估计构造而成（即不是普通的 pdf，它没有函数）。它是多元的，不可分离的，是离散的数据。

我不想构造一个 NxN 矩阵 (F(x,y))，它是这个 pdf 的二维累积分布函数，这样我就可以随机抽样 F(x,y) = P(x

我认为多元函数的 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] 上从这个函数随机生成样本时较高的加权（即更有可能）值以这种方式出现的频率更高，样本收敛于预期值（在多个峰值的情况下），这是粒子滤波器、神经网络等算法的理想结果。

【问题讨论】：

使用trapz 而不是cumsum 可能更准确：mathworks.com.au/help/matlab/ref/trapz.html
@David 据我所知，我将如何使用 'trapz' 生成 CDF 函数，'trapz' 将离散积分评估为数值区域（或 2D 案例的体积）？
哦，当然，抱歉，我没有用。

标签： matlab probability cdf cumsum

【解决方案1】：

首先考虑一维情况。您有一个由向量 F 表示的函数，并且想要进行数值积分。 cumsum(F) 会这样做，但它使用了较差的数值积分形式。即，它将F 视为一个阶跃函数。您可以改为使用Trapezoidal rule 或Simpson's rule 进行更准确的数值积分。

二维情况也不例外。您对cumsum(cumsum(F,1),2) 的使用再次将F 视为阶跃函数，并且该假设导致的数值误差只会随着积分维数的增加而变得更糟。存在梯形规则和辛普森规则的二维类似物。由于这里要重复的数学太多了，请看这里： http://onestopgate.com/gate-study-material/mathematics/numerical-analysis/numerical-integration/2d-trapezoidal.asp.

【讨论】：

我查看了链接，它似乎很好地描述了如何在二维情况下进行数值积分，我仍然不确定如何使用这种方法产生函数（离散曲面函数： = 矩阵 NxN F(x,y))，准确地表示表面积分，然后我可以随机采样？
@Chriso 对于最初的目的，您的cumsum(cumsum(F,1),2) 方法应该可以工作，尽管准确性欠佳。您应该使用从中获得的图片更新您的问题（以便同时有 PDF 和 CDF 的图片），然后详细说明“随机样本”的含义。如果你想要做的是随机抽样一个多元高斯分布，这不是这样做的方法。
没问题，我会用更多信息更新问题。仅供参考，我正在使用多元高斯进行测试，因为它有一个精确的解决方案，我可以检查算法的准确性，但该方法将用于由多个加权内核的卷积构造的任意 pdf，从而产生一个基本上任意形状的 pdf。这个 pdf 除了构建 cdf 和样本 F(x,y) = P(x
哦，我忘了说我对 CDF 进行采样的原因是因为更可能的值更频繁地出现，这是目标。（这是粒子过滤算法的一部分）
@Chriso 如果您所做的只是采样，那么您当前的方法可能已经足够准确了。请注意，您可能可以使用两步方法进行采样：（1）采样内核权重的离散分布以确定内核（2）样本内核分布（多元高斯）。我怀疑您是否真的需要进行数值积分。

【解决方案2】：

您无需计算概率密度函数的二维积分即可从分布中采样。如果您正在计算二维积分，那么您将错误地解决问题。

这里有两种解决抽样问题的方法。

(1) 你写道你有一个核密度估计。核密度估计是混合密度的一种特殊情况。任何混合密度都可以通过首先选择一个内核（可能不同或相同的权重，应用相同的程序），然后从该内核中采样来进行采样。（这适用于任何数量的维度。）通常，内核是一些相对简单的分布，例如高斯分布，因此很容易从中采样。

(2) 任何联合密度 P(X, Y) 等于 P(X | Y) P(Y) （等价于 P(Y | X) P(X)）。因此，您可以从 P(Y)（或 P(X)）采样，然后从 P(X | Y) 采样。为了从 P(X | Y) 中采样，您需要沿 Y = y 线对 P(X, Y) 进行积分（其中 y 是 Y 的采样值），但是（这很关键）您只需要沿着这条线整合；您不需要整合 X 和 Y 的所有值。

如果您告诉我们更多有关您的问题的信息，我可以提供详细信息。

【讨论】：