Jak nazywa się metoda szacowania gęstości, w której wszystkie możliwe pary są używane do utworzenia rozkładu normalnej mieszaniny?

Właśnie pomyślałem o zgrabnym (niekoniecznie dobrym) sposobie tworzenia szacunków gęstości jednowymiarowej i moje pytanie brzmi:

Czy ta metoda szacowania gęstości ma nazwę? Jeśli nie, to czy jest to szczególny przypadek innej metody w literaturze?

Oto metoda: mamy wektor który, jak zakładamy, pochodzi z nieznanego rozkładu, który chcielibyśmy oszacować. Sposobem na zrobienie tego jest pobranie wszystkich możliwych par wartości w i dla każdej pary dopasować rozkład normalny z maksymalnym prawdopodobieństwem. Wynikowa ocena gęstości jest wówczas rozkładem mieszaniny, który składa się ze wszystkich uzyskanych normalnych, gdzie każda normalna ma taką samą wagę.X [ x i , x j ] i ≠ $X = [x_1,x_2,...,x_n]$$X$$[x_i,x_j]_{i \neq j}$

Poniższy rysunek ilustruje użycie tej metody na wektorze . Tutaj koła są punktami danych, kolorowe Normalne są maksymalnymi rozkładami prawdopodobieństwa oszacowanymi przy użyciu każdej możliwej pary, a gruba czarna linia pokazuje wynikowe oszacowanie gęstości (to znaczy rozkład mieszaniny).$[-1.3,0.15,0.73,1.4]$

Nawiasem mówiąc, bardzo łatwo jest wdrożyć metodę w R, która pobiera próbkę z wynikowego rozkładu mieszaniny:

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

Jest to intrygujący pomysł, ponieważ estymator odchylenia standardowego wydaje się być mniej wrażliwy na wartości odstające niż zwykłe podejścia średniej kwadratowej. Wątpię jednak, aby ten estymator został opublikowany. Są trzy powody, dla których: jest nieefektywny obliczeniowo, jest tendencyjny, a nawet po poprawieniu błędu systematycznego jest nieefektywny statystycznie (ale tylko trochę). Można to zobaczyć z małą wstępną analizą, więc zróbmy to najpierw, a następnie wyciągnij wnioski.

Analiza

Estymatory ML średniej i odchylenia standardowego na podstawie danych wynosząσ ( x i , x j$\mu$$\sigma$$(x_i, x_j)$

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

i

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

Dlatego metoda opisana w pytaniu to

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

który jest zwykle estymatorem średniej, oraz

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

Oczekiwaną wartość tego estymatora można łatwo znaleźć, wykorzystując wymienność danych, co oznacza, że jest niezależne od i . Skąd$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.$$

Ale ponieważ i są niezależnymi zmiennymi Normalne, ich różnica jest równa zeru Normalna z wariancją . Jego wartość bezwzględna wynosi zatem razy , którego średnia to . w konsekwencji$x_i$$x_j$$2\sigma^2$$\sqrt{2}\sigma$$\chi(1)$$\sqrt{2/\pi}$

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

Współczynnik jest odchyleniem w tym estymatorze.$2/\sqrt{\pi} \approx 1.128$

W ten sam sposób, ale przy znacznie większej ilości pracy, można było obliczyć wariancję , ale - jak zobaczymy - mało prawdopodobne jest zainteresowanie tym tematem, dlatego oszacuję to za pomocą szybkiej symulacji .$\hat\sigma$

Wnioski

1. Estymator jest stronniczy. ma znaczną stałą stronniczość około + 13%. Można to poprawić. W tym przykładzie z wielkością próby zarówno histogramy z tendencyjnością, jak i z korekcją błędu systematycznego. Widoczny jest błąd 13%.$\hat\sigma$$n=20,000$

   

2. Jest to nieefektywne obliczeniowo. Ponieważ suma wartości bezwzględnych,, nie ma algebraicznego uproszczenia, jego obliczenie wymaga wysiłku zamiast wysiłku dla prawie każdego innego estymatora. To źle się skaluje, co powoduje, że jest zbyt drogie, gdy przekroczy około . Na przykład obliczenie poprzedniej liczby wymagało 45 sekund czasu procesora i 8 GB pamięci RAM . (Na innych platformach wymagania pamięci RAM byłyby znacznie mniejsze, być może przy niewielkim koszcie czasu obliczeniowego).$\sum_{i,j}|x_i-x_j|$$O(n^2)$$O(n)$$n$$10,000$R

3. Jest to statystycznie nieefektywne. Aby uzyskać najlepszy wynik, rozważmy wersję bezstronną i porównajmy ją z bezstronną wersją estymatora najmniejszych kwadratów lub maksymalnego prawdopodobieństwa

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

   Poniższy Rkod pokazuje, że bezstronna wersja estymatora w pytaniu jest zaskakująco wydajna: w zakresie wielkości próbek od do jej wariancja jest zwykle około 1% do 2% większa niż wariancja . Oznacza to, że powinieneś zaplanować zapłacenie dodatkowych 1% do 2% więcej za próbki, aby osiągnąć dowolny poziom precyzji w szacowaniu .n = 300 σ O L S$n=3$$n=300$$\hat\sigma_{OLS}$$\sigma$

Potem

Forma przypomina solidny i odporny estymator Theil-Sen - ale zamiast używać median różnic bezwzględnych, używa ich środków. Jeśli celem jest posiadanie estymatora odpornego na wartości odstające lub odpornego na odstępstwa od założenia Normalności, wówczas użycie mediany byłoby bardziej wskazane. $\hat\sigma$

---

Kod

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