Wie heißt die Dichteschätzmethode, bei der alle möglichen Paare verwendet werden, um eine normale Mischungsverteilung zu erstellen?

Ich habe mir gerade eine nette (nicht unbedingt gute) Methode ausgedacht, um eindimensionale Dichteschätzungen zu erstellen, und meine Frage lautet:

Hat diese Dichteschätzmethode einen Namen? Wenn nicht, handelt es sich um einen Sonderfall einer anderen Methode in der Literatur?

Hier ist die Methode: Wir haben ein Vektor $X = [x_1,x_2,...,x_n]$ dem wir annehmen, dass es aus einer unbekannten Verteilung stammt, die wir schätzen möchten. Eine Möglichkeit, dies zu tun, besteht darin, alle möglichen Wertepaare in $X$ und für jedes Paar $[x_i,x_j]_{i \neq j}$ eine Normalverteilung mit maximaler Wahrscheinlichkeit anzupassen. Die resultierende Dichteschätzung ist dann die Mischungsverteilung, die aus allen resultierenden Normalen besteht, wobei jede Normale gleich gewichtet wird.

Die folgende Abbildung zeigt die Verwendung dieser Methode auf dem Vektor . Hier sind die Kreise die Datenpunkte, die farbigen Normalen die mit jedem möglichen Paar geschätzten maximalen Wahrscheinlichkeitsverteilungen und die dicke schwarze Linie zeigt die resultierende Dichteschätzung (dh die Mischungsverteilung).$[-1.3,0.15,0.73,1.4]$

Übrigens ist es wirklich einfach, eine Methode in R zu implementieren, die eine Probe aus der resultierenden Mischungsverteilung zieht:

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

Dies ist eine faszinierende Idee, da der Schätzer der Standardabweichung weniger empfindlich gegenüber Ausreißern zu sein scheint als die üblichen Root-Mean-Square-Ansätze. Ich bezweifle jedoch, dass dieser Schätzer veröffentlicht wurde. Es gibt drei Gründe: Es ist rechnerisch ineffizient, es ist voreingenommen, und selbst wenn die Voreingenommenheit korrigiert wird, ist es statistisch ineffizient (aber nur wenig). Diese können mit einer kleinen vorläufigen Analyse gesehen werden, also lasst uns das zuerst tun und dann die Schlussfolgerungen ziehen.

Analyse

Der ML - Schätzer der Mittelwert und eine Standardabweichung & sgr; auf Daten basieren ( x i , x j ) sind$\mu$$\sigma$$(x_i, x_j)$

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

und

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

Daher ist die in der Frage beschriebene Methode

$$\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,$$

das ist der übliche Schätzer des Mittelwertes, und

$$\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|.$$

Der erwartete Wert dieses Schätzers kann leicht durch Ausnutzen der Austauschbarkeit der Daten ermittelt werden, was impliziert$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$ ist die Vorspannung in diesem Schätzer.

$\hat\sigma$ , aber - wie wir sehen werden - es ist unwahrscheinlich , viel Interesse daran sein, also werde ich es nur schätzen , mit einer schnellen Simulation.

Schlussfolgerungen

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

   

2. $\sum_{i,j}|x_i-x_j|$$O(n^2)$$O(n)$$n$$10,000$ oder so. Für die Berechnung der vorherigen Zahl wurden beispielsweise 45 Sekunden CPU-Zeit und 8 GB RAM benötigtR. (Auf anderen Plattformen wären die RAM-Anforderungen viel geringer, möglicherweise mit geringen Kosten für die Rechenzeit.)

3. Es ist statistisch ineffizient. Um die beste Darstellung zu erzielen, betrachten wir die unvoreingenommene Version und vergleichen sie mit der neutralen Version des Schätzers für kleinste Quadrate oder maximale Wahrscheinlichkeit

   $$\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$

Nachher

$\hat\sigma$

---

Code

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