我需要计算一个称为互信息的度量。首先，我需要计算另一个度量，称为熵，例如 x 和 y 的联合熵：
-∬p(x,y)·log p(x,y)dxdy
因此，为了计算 p(x,y)，我使用了核密度估计器（这样，函数 kde2d，它返回了 Z 值（在该窗口中具有 x 和 y 的概率）。
同样，到目前为止，我有一个矩阵 Z 值 [1x100] x [1x100]，这等于我的 p(x,y)。但我必须通过发现表面下的体积（双积分）来整合它。但我没有找到办法做到这一点。计算双求积的函数quad2d 不起作用，因为我只集成了一个数值矩阵p(x,y)，它给了我一个常数....
任何人都知道找到体积/计算二重积分的东西吗？
剧情图片来自persp3d：
谢谢大家！！！！
【问题讨论】：
标签： r estimation entropy kernel-density
从kde2d 获得结果后，计算数值积分就非常简单了。下面的示例会话概述了如何做到这一点。
如您所知，数值双积分只是二维求和。 kde2d，默认将range(x) 和range(y) 作为二维域。我看到你有一个 100*100 的矩阵，所以我认为你在使用 kde2d 时设置了 n = 100。现在，kde$x、kde$y 定义了一个 100 * 100 的网格，den$z 给出了每个网格单元的密度。计算每个网格单元的大小很容易（它们都相等），然后我们执行三个步骤：
z * log(z)；因为z 是一个 100 * 100 的矩阵，所以这也是一个矩阵。您只需将它们相加，然后将其乘以单元格大小cell_size，即可得到非归一化熵；## sample data: bivariate normal, with covariance/correlation 0
set.seed(123); x <- rnorm(1000, 0, 2) ## marginal variance: 4
set.seed(456); y <- rnorm(1000, 0, 2) ## marginal variance: 4
## load MASS
library(MASS)
## domain:
xlim <- range(x)
ylim <- range(y)
## 2D Kernel Density Estimation
den <- kde2d(x, y, n = 100, lims = c(xlim, ylim))
##persp(den$x,den$y,den$z)
z <- den$z ## extract density
## den$x, den$y expands a 2D grid, with den$z being density on each grid cell
## numerical integration is straighforward, by aggregation over all cells
## the size of each grid cell (a rectangular cell) is:
cell_size <- (diff(xlim) / 100) * (diff(ylim) / 100)
## normalizing constant; ideally should be 1, but actually only close to 1 due to discretization
norm <- sum(z) * cell_size
## your integrand: z * log(z) * (-1):
integrand <- z * log(z) * (-1)
## get numerical integral by summation:
entropy <- sum(integrand) * cell_size
## self-normalization:
entropy <- entropy / norm
验证
上面的代码给出了4.230938的熵。现在，Wikipedia - Multivariate normal distribution 给出了熵公式：
(k / 2) * (1 + log(2 * pi)) + (1 / 2) * log(det(Sigma))
对于上述二元正态分布，我们有k = 2。我们有Sigma（协方差矩阵）：
4 0
0 4
其行列式为16。因此，理论值为：
(1 + log(2 * pi)) + (1 / 2) * log(16) = 4.224171
好匹配！
【讨论】：