计算内核双变量密度估计图下的体积
Calculate the volume under a plot of kernel bivariate density estimation
我需要计算一种称为互信息的度量。首先，我需要计算另一个度量，称为熵，例如x和y的联合熵：
-∬p(x,y)·log p(x,y)dxdy
因此，为了计算 p(x,y)，我使用了核密度估计器（以这种方式，函数 kde2d，它返回 Z 值（x 和 y 的概率window).
同样，到目前为止，我有一个 Z 值 [1x100] x [1x100] 的矩阵，这等于我的 p(x,y)。但我必须通过发现表面下的体积（双积分）来整合它。但我没有找到办法做到这一点。函数 quad2d 计算二重积分没有用，因为我只积分了一个数值矩阵 p(x,y)，它给了我一个常数....
有人知道如何找到 volume/calculate 二重积分吗？
来自persp3d的剧情图片：
谢谢大家!!!!
获得 kde2d 的结果后，计算数值积分非常简单。下面的示例会话概述了如何操作。
如您所知，数值二重积分只是二维求和。 kde2d默认以range(x)和range(y)作为二维域。我看到你有一个 100*100 的矩阵，所以我认为你在使用 kde2d 时设置了 n = 100。现在，kde$x、kde$y 定义了一个 100 * 100 的网格，den$z 给出了每个网格单元的密度。很容易计算每个网格单元的大小（它们都是相等的），然后我们做三个步骤：
- 找到归一化常数；虽然我们知道理论上密度总和（或积分）为1，但经过计算机离散化后，它只近似于1。所以我们先计算这个归一化常数，以便稍后重新缩放；
- 熵的被积函数是
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
好匹配！
我需要计算一种称为互信息的度量。首先，我需要计算另一个度量，称为熵，例如x和y的联合熵：
-∬p(x,y)·log p(x,y)dxdy
因此，为了计算 p(x,y)，我使用了核密度估计器（以这种方式，函数 kde2d，它返回 Z 值（x 和 y 的概率window).
同样，到目前为止，我有一个 Z 值 [1x100] x [1x100] 的矩阵，这等于我的 p(x,y)。但我必须通过发现表面下的体积（双积分）来整合它。但我没有找到办法做到这一点。函数 quad2d 计算二重积分没有用，因为我只积分了一个数值矩阵 p(x,y)，它给了我一个常数....
有人知道如何找到 volume/calculate 二重积分吗？
来自persp3d的剧情图片：
谢谢大家!!!!
获得 kde2d 的结果后，计算数值积分非常简单。下面的示例会话概述了如何操作。
如您所知，数值二重积分只是二维求和。 kde2d默认以range(x)和range(y)作为二维域。我看到你有一个 100*100 的矩阵，所以我认为你在使用 kde2d 时设置了 n = 100。现在，kde$x、kde$y 定义了一个 100 * 100 的网格，den$z 给出了每个网格单元的密度。很容易计算每个网格单元的大小（它们都是相等的），然后我们做三个步骤：
- 找到归一化常数；虽然我们知道理论上密度总和（或积分）为1，但经过计算机离散化后，它只近似于1。所以我们先计算这个归一化常数，以便稍后重新缩放；
- 熵的被积函数是
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
好匹配！