Интегрирование оценки плотности ядра в 2D

Я исхожу из этого вопроса на случай, если кто-нибудь захочет пойти по следу.

По сути, у меня есть набор данных состоящий из объектов, где к каждому объекту прикреплено заданное количество измеренных значений (два в данном случае):$\Omega$$N$

Мне нужен способ определить вероятность того, что новый объект принадлежать поэтому в этом вопросе мне посоветовали получить плотность вероятности через оценщик плотности ядра, который, как я полагаю, уже есть.$p[x_p, y_p]$$\Omega$$\hat{f}$

Поскольку моя цель - получить вероятность того, что этот новый объект ( ) будет принадлежать этому двухмерному набору данных , мне было сказано интегрировать значения pdf over " поддержки, для которой плотность меньше той, которую вы наблюдали ". «Наблюдаемая» плотность оценивается в в новом объекте , то есть: . Поэтому мне нужно решить уравнение:$p[x_p, y_p]$$\Omega$$\hat{f}$$\hat{f}$$p$$\hat{f}(x_p, y_p)$

PDF моего 2D-набора данных (полученный через модуль python stats.gaussian_kde ) выглядит следующим образом:

где красная точка представляет новый объект нанесенный поверх PDF моего набора данных.$p[x_p, y_p]$

Итак, вопрос: как я могу вычислить вышеуказанный интеграл для пределов когда pdf выглядит так?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

добавлять

Я провел несколько тестов, чтобы понять, насколько хорошо работает метод Монте-Карло, о котором я упоминал в одном из комментариев. Вот что я получил:

Значения, кажется, изменяются немного больше для областей с более низкой плотностью, причем обе полосы пропускания показывают более или менее одинаковое изменение. Наибольшее изменение в таблице происходит для точки (x, y) = (2.4,1.5), сравнивая выборочное значение Сильвермана 2500 против 1000, что дает разницу 0.0126или ~1.3%. В моем случае это было бы в значительной степени приемлемо.

Edit : Я просто заметилчто в 2 измерении правило Скотта эквивалентно Silverman в соответствии с определениемприведенным здесь .

Простой способ - растеризовать область интегрирования и вычислить дискретное приближение к интегралу.

Есть некоторые вещи, на которые стоит обратить внимание:

1. Удостоверьтесь, что вы охватили больше, чем степень точек: вам нужно включить все места, где оценка плотности ядра будет иметь какие-либо заметные значения. Это означает, что вам нужно увеличить экстент точек в три-четыре раза по пропускной способности ядра (для ядра Гаусса).

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

Вот иллюстрация для набора данных 256 точек:

Точки показаны черными точками, наложенными на две оценки плотности ядра. Шесть больших красных точек - это «зонды», на которых оценивается алгоритм. Это было сделано для четырех полос пропускания (по умолчанию от 1,8 (по вертикали) до 3 (по горизонтали), 1/2, 1 и 5 единиц) с разрешением 1000 на 1000 ячеек. Следующая матрица диаграммы рассеяния показывает, насколько сильно результаты зависят от ширины полосы для этих шести точек измерения, которые охватывают широкий диапазон плотностей:

Изменение происходит по двум причинам. Очевидно, что оценки плотности отличаются, вводя одну форму изменения. Что еще более важно, различия в оценках плотности могут создать большие различия в любой отдельной («пробной») точке. Последний вариант наиболее велик вокруг «полос» скоплений точек средней плотности - именно в тех местах, где этот расчет, вероятно, будет использоваться чаще всего.

Это демонстрирует необходимость существенной осторожности при использовании и интерпретации результатов этих вычислений, поскольку они могут быть настолько чувствительны к относительно произвольному решению (используемой полосе пропускания).

Код R

Алгоритм содержится в полдюжины строк первой функции f. Чтобы проиллюстрировать его использование, остальная часть кода генерирует предыдущие цифры.

library(MASS)     # kde2d
library(spatstat) # im class
f <- function(xy, n, x, y, ...) {
  #
  # Estimate the total where the density does not exceed that at (x,y).
  #
  # `xy` is a 2 by ... array of points.
  # `n`  specifies the numbers of rows and columns to use.
  # `x` and `y` are coordinates of "probe" points.
  # `...` is passed on to `kde2d`.
  #
  # Returns a list:
  #   image:    a raster of the kernel density
  #   integral: the estimates at the probe points.
  #   density:  the estimated densities at the probe points.
  #
  xy.kde <- kde2d(xy[1,], xy[2,], n=n, ...)
  xy.im <- im(t(xy.kde$z), xcol=xy.kde$x, yrow=xy.kde$y) # Allows interpolation $
  z <- interp.im(xy.im, x, y)                            # Densities at the probe points
  c.0 <- sum(xy.kde$z)                                   # Normalization factor $
  i <- sapply(z, function(a) sum(xy.kde$z[xy.kde$z < a])) / c.0
  return(list(image=xy.im, integral=i, density=z))
}
#
# Generate data.
#
n <- 256
set.seed(17)
xy <- matrix(c(rnorm(k <- ceiling(2*n * 0.8), mean=c(6,3), sd=c(3/2, 1)), 
               rnorm(2*n-k, mean=c(2,6), sd=1/2)), nrow=2)
#
# Example of using `f`.
#
y.probe <- 1:6
x.probe <- rep(6, length(y.probe))
lims <- c(min(xy[1,])-15, max(xy[1,])+15, min(xy[2,])-15, max(xy[2,]+15))
ex <- f(xy, 200, x.probe, y.probe, lim=lims)
ex$density; ex$integral
#
# Compare the effects of raster resolution and bandwidth.
#
res <- c(8, 40, 200, 1000)
system.time(
  est.0 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, lims=lims)$integral))
est.0
system.time(
  est.1 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1, lims=lims)$integral))
est.1
system.time(
  est.2 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1/2, lims=lims)$integral))
est.2
system.time(
  est.3 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=5, lims=lims)$integral))
est.3
results <- data.frame(Default=est.0[,4], Hp5=est.2[,4], 
                      H1=est.1[,4], H5=est.3[,4])
#
# Compare the integrals at the highest resolution.
#
par(mfrow=c(1,1))
panel <- function(x, y, ...) {
  points(x, y)
  abline(c(0,1), col="Red")
}
pairs(results, lower.panel=panel)
#
# Display two of the density estimates, the data, and the probe points.
#
par(mfrow=c(1,2))
xy.im <- f(xy, 200, x.probe, y.probe, h=0.5)$image
plot(xy.im, main="Bandwidth=1/2", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

xy.im <- f(xy, 200, x.probe, y.probe, h=5)$image
plot(xy.im, main="Bandwidth=5", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

Если у вас есть приличное количество наблюдений, вам, возможно, не понадобится делать какую-либо интеграцию вообще. Скажите, что ваша новая точка - . Предположим, у вас есть оценщик плотности ; Суммируйте количество наблюдений для которых и разделите на размер выборки. Это дает вам приближение к требуемой вероятности.${\bf{x}}_0$$\hat f$${\bf{x}}$$\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$

Это предполагает, что не является «слишком маленьким», а размер вашей выборки достаточно велик (и достаточно разбросан), чтобы дать достойную оценку в регионах с низкой плотностью. Тем не менее, 20000 случаев кажутся достаточно большими для двумерного .$\hat f({\bf{x}}_0)$${\bf{x}}$