Как называется метод оценки плотности, при котором все возможные пары используются для создания нормального распределения смеси?

12

Я просто подумал о аккуратном (не обязательно хорошем) способе создания одномерных оценок плотности, и мой вопрос:

У этого метода оценки плотности есть имя? Если нет, то является ли это частным случаем какого-либо другого метода в литературе?

Вот метод: Мы имеем вектор который мы предполагаем, взят из некоторого неизвестного распределения, которое мы хотели бы оценить. Способ сделать это , чтобы принять все возможные пары значений X и для каждой пары [ х я , х J ] я J соответствуют нормальному распределению с использованием максимального правдоподобия. Результирующая оценка плотности представляет собой распределение смеси, которое состоит из всех результирующих нормалей, где каждому нормальному значению присваивается равный вес.X=[x1,x2,...,xn]X[xi,xj]ij

На рисунке ниже показано использование этого метода для вектора . Здесь кружки - точки данных, цветные нормали - распределения максимального правдоподобия, оцененные с использованием каждой возможной пары, а жирная черная линия показывает итоговую оценку плотности (то есть распределение смеси).[1.3,0.15,0.73,1.4]

введите описание изображения здесь

Кстати, действительно легко реализовать метод в R, который берет образец из полученного распределения смеси:

# Generating some "data"
x <- rnorm(30)

# Drawing from the density estimate using the method described above.
density_estimate_sample <- replicate(9999, {
  pair <- sample(x, size = 2)
  rnorm(1, mean(pair), sd(pair))
})

# Plotting the density estimate compared with 
# the "data" and the "true" density.
hist(x ,xlim=c(-5, 5), main='The "data"')
hist(density_estimate_sample, xlim=c(-5, 5), main='Estimated density')
hist(rnorm(9999), xlim=c(-5, 5), main='The "true" density')

введите описание изображения здесь

Расмус Батх
источник
5
Попробуйте свой метод, используяx <- c(rnorm(30), rnorm(30, 10))
Dason
2
@ Dason Да, в этом случае метод не работает вообще! :) Тоже не сходится с большим n.
Расмус Батх
4
Это похоже на искаженную версию оценки плотности ядра, где пропускная способность оценивается перекрестной проверкой!
Сиань
X=[x1,x2,,xn]n

Ответы:

6

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

Анализ

μσ(xi,xj)

μ^(xi,xj)=xi+xj2

и

σ^(xi,xj)=|xixj|2.

Поэтому метод, описанный в вопросе

μ^(x1,x2,,xn)=2n(n1)i>jxi+xj2=1ni=1nxi,

которая является обычной оценкой среднего значения, и

σ^(x1,x2,,xn)=2n(n1)i>j|xixj|2=1n(n1)i,j|xixj|.

E=E(|xixj|)ij

E(σ^(x1,x2,,xn))=1n(n1)i,jE(|xixj|)=E.

xixj2σ22σχ(1)2/π

E=2πσ.

2/π1.128

σ^

Выводы

  1. σ^n=20,000

    фигура

  2. i,j|xixj|O(n2)O(n)n10,000R, (На других платформах требования к ОЗУ будут намного меньше, возможно, с небольшими затратами времени на вычисления).

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

    σ^OLS=(1n1i=1n(xiμ^)2)(n1)Γ((n1)/2)2Γ(n/2).

    Rn=3n=300σ^OLSσ

позже

σ^


Код

sigma <- function(x) sum(abs(outer(x, x, '-'))) / (2*choose(length(x), 2))
#
# sigma is biased.
#
y <- rnorm(1e3) # Don't exceed 2E4 or so!
mu.hat <- mean(y)
sigma.hat <- sigma(y)

hist(y, freq=FALSE,
     main="Biased (dotted red) and Unbiased (solid blue) Versions of the Estimator",
     xlab=paste("Sample size of", length(y)))
curve(dnorm(x, mu.hat, sigma.hat), col="Red", lwd=2, lty=3, add=TRUE)
curve(dnorm(x, mu.hat, sqrt(pi/4)*sigma.hat), col="Blue", lwd=2, add=TRUE)
#
# The variance of sigma is too large.
#
N <- 1e4
n <- 10
y <- matrix(rnorm(n*N), nrow=n)
sigma.hat <- apply(y, 2, sigma) * sqrt(pi/4)
sigma.ols <- apply(y, 2, sd) / (sqrt(2/(n-1)) * exp(lgamma(n/2)-lgamma((n-1)/2)))

message("Mean of unbiased estimator is ", format(mean(sigma.hat), digits=4))
message("Mean of unbiased OLS estimator is ", format(mean(sigma.ols), digits=4))
message("Variance of unbiased estimator is ", format(var(sigma.hat), digits=4))
message("Variance of unbiased OLS estimator is ", format(var(sigma.ols), digits=4))
message("Efficiency is ", format(var(sigma.ols) / var(sigma.hat), digits=4))
Whuber
источник
Соответствующая литература восходит некоторое время назад, например, Даунтон, Ф. 1966 Линейные оценки с полиномиальными коэффициентами. Биометрика 53: 129-141 дои: 10.1093 / биомет / 53.1-2.129
Ник Кокс
Вау, я получил больше, чем рассчитывал! :)
Расм Бхат