R: вероятность/числовой интеграл двумерной (или многомерной) плотности ядра

Я использую пакет ks для оценки плотности ядра. Вот простой пример:

n <- 70
x <- rnorm(n)

library(ks)
f_kde <- kde(x)

На самом деле меня интересуют соответствующие вероятности превышения моих входных данных, которые могут быть легко возвращены ks с f_kde:

p_kde <- pkde(x, f_kde)

Это делается в ks с числовым интегрированием по правилу Симпсона. К сожалению, они реализовали это только для случая 1d. В двумерном случае в ks нет реализации какого-либо метода для возврата вероятностей:

y <- rnorm(n)
f_kde <- kde(data.frame(x,y))
# does not work, but it's what I am looking for:
p_kde <- pkde(data.frane(x,y), f_kde)

Я не смог найти какой-либо пакет или помощь в поиске в stackoverflow для решения этой проблемы в R (существуют некоторые предложения для Python, но я хотел бы сохранить его в R). Приветствуется любая строка кода или рекомендация по пакету. Несмотря на то, что меня в основном интересует двумерный случай, любые идеи для многомерного случая также приветствуются.


r kde cdf
person Felix Phl    schedule 18.07.2020    source источник
comment
Разве вы не можете просто сделать двойной интеграл, как в этот предыдущий вопрос   -  person G5W    schedule 18.07.2020
comment
@ G5W Я уже понял, что pracma::simpson2d делает то, что я ищу, но для этого требуется функция в качестве входных данных, и я не знаю, как самостоятельно кодировать двумерную плотность ядра, а также понятия не имею, с чего начать.   -  person Felix Phl    schedule 19.07.2020

Ответы (1)


arrow_upward
3
arrow_downward

kde позволяет использовать многомерную оценку ядра, поэтому мы можем использовать kde для вычисления pkde.
Для этого мы вычисляем kde на достаточно малых dx и dy шагах, используя параметр eval.points: это дает нам оценку локальной плотности на dx*dy квадрате.
Проверяем, что сумма оценок перемножается по поверхности квадратов почти равно 1:

library(ks)
set.seed(1)
n <- 10000
x <- rnorm(n)
y <- rnorm(n)
xy <- cbind(x,y)

xmin <- -10
xmax <- 10
dx <- .1

ymin <- -10
ymax <- 10
dy <- .1

pts.x <- seq(xmin, xmax, dx)
pts.y <- seq(ymin, ymax, dy)
pts <- as.data.frame(expand.grid(x = pts.x, y = pts.y))
f_kde <- kde(xy,eval.points=pts)

pts$est <- f_kde$estimate

sum(pts$est)*dx*dy
[1] 0.9998778

Теперь вы можете запросить фрейм данных pts для кумулятивной вероятности в выбранной вами области:

library(data.table)
setDT(pts)
# cumulative density
pts[x < 1 & y < 2 , .(pkde=sum(est)*dx*dy)]
        pkde
1: 0.7951228

# average density around a point
tolerance <-.1
pts[pmin(abs(x-1))<tolerance & pmin(abs(y-2))<tolerance, .(kde = mean(est))]
          kde
1: 0.01465478
person Waldi    schedule 20.07.2020
comment
Спасибо за ваш ответ. У меня есть два вопроса: 1. Я тестировал код со своими данными (n = 69), и сумма p составляет всего около 0,9, стоит ли мне об этом беспокоиться? 2. У меня возникли проблемы с преобразованием последней строки для моей цели, не могли бы вы показать, как оценить плотность вероятности xy[1,] (= x[1] и y[1]), спасибо! - person Felix Phl; 21.07.2020
comment
Что касается вопроса 1) kde только с 69 точками для растра 100 * 100 все еще очень неточен, поэтому я не удивлен, что вы получили сумму около 0,9. Для 2) pkde - это кумулятивная плотность, поэтому я сделал сумму двух неравенств. Если вам нужна плотность в одной точке, вы можете использовать pts[ x==1 & y==1] при условии, что запрашиваемая точка находится в растре. - person Waldi; 21.07.2020