Вычислить объем под графиком двумерной оценки плотности ядра

Получив результаты kde2d, очень просто вычислить числовой интеграл. В приведенном ниже примере сеанса показано, как это сделать.

Как вы знаете, числовой двойной интеграл — это просто двумерное суммирование. kde2d по умолчанию принимает range(x) и range(y) в качестве 2D-домена. Я вижу, что у вас есть матрица 100 * 100, поэтому я думаю, что вы установили n = 100 при использовании kde2d. Теперь kde$x, kde$y определяют сетку 100 * 100, а den$z задает плотность в каждой ячейке сетки. Легко вычислить размер каждой ячейки сетки (все они равны), тогда делаем три шага:

1. найти нормирующие константы; хотя мы знаем, что теоретически плотность дает в сумме (или интегрирует) 1, но после компьютерной дискретизации она только приближается к 1. Итак, мы сначала вычисляем эту нормализующую константу для последующего масштабирования;
2. подынтегральное выражение для энтропии равно z * log(z); так как z это матрица 100*100, то это тоже матрица. Вы просто суммируете их и умножаете на размер ячейки cell_size, тогда вы получаете ненормализованную энтропию;
3. перемасштабировать ненормализованную энтропию на нормализованную.

## 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. Теперь Википедия — Многомерное нормальное распределение дает формулу энтропии:

(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

Хорошая партия!