2016-05-20 10 views
1

Мне нужно рассчитать меру, называемую взаимной информацией. Прежде всего, мне нужно вычислить другую меру, называемую энтропией, например, совместную энтропию х и у:Рассчитать объем по графику оценки двумерной плотности ядра

-∬p(x,y)·log p(x,y)dxdy 

Таким образом, чтобы вычислить p(x,y), я использовал оценщик плотности ядра (в этом случае, функция kde2d, и она вернула Z значения (вероятность наличия х и у в этом окне).

Опять же, теперь у меня есть матрица Z значений [1x100] x [1x100], это равно мой p(x,y). Но я должен интегрировать его, открыв объем под поверхностью (интеграл doble), но я не нашел способа сделать это. Функция quad2d, чтобы вычислить двойной квадратурной не работает, потому что я только интегрированный числовой матрицы p(x,y), и это дает мне постоянный ....

Каждый знает, что-то найти, что объем/вычислить двойной интеграл?

Изображение сюжета из persp3d:

density estimate

Спасибо всем !!!!

ответ

3

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

Как вы знаете, числовой двойной интеграл представляет собой просто суммирование по 2D. 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. Теперь, 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 

Хороший матч!

+0

Было бы проще, если бы я дискретировал _a priori_ и подсчитал частоту на ширине бункеров? Будет ли это производить гораздо больше ошибок? Благодарю. –

+0

Я тестировал ваш метод против гистограммы, и что-то неправильное может произойти. Когда я тестирую ваш метод на 1000 случайных выборок (как вы объясните выше), энтропия дает мне значения около нуля. Но, по сравнению с теорией, это неверно, потому что это должно дать мне log (n). Я использовал какой-то другой пакет энтропии, и он дает правильный результат ~ log (n). Вы знаете, что происходит? Благодаря! –

+0

Спасибо! Это справедливо для гауссовских нормальных ядер, но как насчет неопределенных плотностей? Тест на нормальную плотность должен давать хорошие результаты, но в моем случае при тестировании на финансовых временных рядах (не стохастическом) эти результаты значительно изменились, особенно из-за полосы пропускания (defaut - гауссовский). Я думаю, что у меня проблемы с хвостами дистрибутива. Проблемы, связанные с ошибками, состоят в том, что у меня есть обширная матрица для вычисления в попарно (временные ряды), поэтому я буду рассматривать большую визуализацию. Еще раз спасибо за помощь! –