¿Cuál es el nombre del método de estimación de densidad donde se usan todos los pares posibles para crear una distribución de mezcla Normal?

Acabo de pensar en una forma ordenada (no necesariamente buena) de crear estimaciones de densidad unidimensionales y mi pregunta es:

¿Este método de estimación de densidad tiene un nombre? Si no, ¿es un caso especial de algún otro método en la literatura?

Aquí está el método: Tenemos un vector que suponemos se extrae de alguna distribución desconocida que nos gustaría estimar. Una forma de hacerlo es tomar todos los pares de valores posibles en X y para cada par [ x i , x j ] i ≠ j ajustar una distribución Normal usando la máxima verosimilitud. La estimación de densidad resultante es entonces la distribución de la mezcla que consiste en todas las normales resultantes, donde a cada normal se le asigna el mismo peso.$X = [x_1,x_2,...,x_n]$$X$$[x_i,x_j]_{i \neq j}$

La siguiente figura ilustra el uso de este método en el vector . Aquí los círculos son los puntos de datos, las normales de colores son las distribuciones de máxima probabilidad estimadas usando cada par posible y la línea negra gruesa muestra la estimación de densidad resultante (es decir, la distribución de la mezcla).$[-1.3,0.15,0.73,1.4]$

Por cierto, es realmente fácil implementar un método en R que extraiga una muestra de la distribución resultante de la mezcla:

# 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')

Esta es una idea intrigante, porque el estimador de la desviación estándar parece ser menos sensible a los valores atípicos que los enfoques habituales de raíz cuadrática media. Sin embargo, dudo que este estimador haya sido publicado. Hay tres razones por las cuales: es computacionalmente ineficiente, está sesgado e incluso cuando se corrige el sesgo, es estadísticamente ineficiente (pero solo un poco). Esto se puede ver con un pequeño análisis preliminar, así que hagamos eso primero y luego saquemos las conclusiones.

Análisis

$\mu$$\sigma$$(x_i, x_j)$

$$\hat\mu(x_i,x_j) = \frac{x_i+x_j}{2}$$

y

$$\hat\sigma(x_i,x_j) = \frac{|x_i-x_j|}{2}.$$

Por lo tanto, el método descrito en la pregunta es

$$\hat\mu(x_1, x_2, \ldots, x_n) = \frac{2}{n(n-1)} \sum_{i\gt j} \frac{x_i+x_j}{2} = \frac{1}{n}\sum_{i=1}^nx_i,$$

cuál es el estimador habitual de la media, y

$$\hat\sigma(x_1, x_2, \ldots, x_n) = \frac{2}{n(n-1)}\sum_{i\gt j}\frac{|x_i-x_j|}{2} = \frac{1}{n(n-1)}\sum_{i,j}|x_i-x_j|.$$

$E = \mathbb{E}(|x_i-x_j|)$$i$$j$

$$\mathbb{E}(\hat\sigma(x_1, x_2, \ldots, x_n)) = \frac{1}{n(n-1)}\sum_{i,j}\mathbb{E}(|x_i-x_j|) = E.$$

$x_i$$x_j$$2\sigma^2$$\sqrt{2}\sigma$$\chi(1)$$\sqrt{2/\pi}$

$$E = \frac{2}{\sqrt{\pi}} \sigma.$$

$2/\sqrt{\pi} \approx 1.128$

$\hat\sigma$

Conclusiones

1. $\hat\sigma$$n=20,000$

   

2. $\sum_{i,j}|x_i-x_j|$$O(n^2)$$O(n)$$n$$10,000$R. (En otras plataformas, los requisitos de RAM serían mucho menores, tal vez a un bajo costo en tiempo de cálculo).

3. Es estadísticamente ineficiente. Para darle el mejor resultado, consideremos la versión imparcial y compárela con la versión imparcial del estimador de mínimos cuadrados o de máxima verosimilitud

   $$\hat\sigma_{OLS} = \sqrt{\left(\frac{1}{n-1} \sum_{i=1}^n \left(x_i - \hat\mu\right)^2\right)} \frac{(n-1)\Gamma((n-1)/2)}{2\Gamma(n/2)}.$$

   R$n=3$$n=300$$\hat\sigma_{OLS}$$\sigma$

Después

$\hat\sigma$

---

Código

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))