Dlaczego mój interwał ładowania jest tak zły?

29

Chciałem zrobić demonstrację klasową, w której porównuję przedział t z przedziałem ładowania początkowego i obliczę prawdopodobieństwo pokrycia obu. Chciałem, aby dane pochodziły z przekrzywionej dystrybucji, więc postanowiłem wygenerować dane jako exp(rnorm(10, 0, 2)) + 1próbkę o wielkości 10 z przesuniętego logarytmu normalnego. Napisałem skrypt, aby narysować 1000 próbek i dla każdej próbki obliczyć zarówno 95% przedział t, jak i 95% przedział procentowy ładowania początkowego na podstawie 1000 powtórzeń.

Kiedy uruchamiam skrypt, obie metody dają bardzo podobne odstępy czasu i obie mają prawdopodobieństwo pokrycia 50-60%. Byłem zaskoczony, ponieważ myślałem, że interwał ładowania będzie lepszy.

Moje pytanie brzmi: czy mam

popełniłeś błąd w kodzie?
popełniłeś błąd przy obliczaniu interwałów?
popełniłeś błąd, oczekując, że interwał ładowania będzie miał lepsze właściwości pokrycia?

Czy istnieje również sposób na zbudowanie bardziej niezawodnego elementu CI w tej sytuacji?

 tCI.total <- 0
 bootCI.total <- 0
 m <- 10 # sample size
 true.mean <- exp(2) + 1

for (i in 1:1000){
 samp <- exp(rnorm(m,0,2)) + 1
 tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)

 boot.means <- rep(0,1000)
 for (j in 1:1000) boot.means[j] <- mean(sample(samp,m,replace=T))
 bootCI <- sort(boot.means)[c(0.025*length(boot.means), 0.975*length(boot.means))]

 if (true.mean > min(tCI) & true.mean < max(tCI)) tCI.total <- tCI.total + 1
 if (true.mean > min(bootCI) & true.mean < max(bootCI)) bootCI.total <- bootCI.total + 1 
}
tCI.total/1000     # estimate of t interval coverage probability
bootCI.total/1000  # estimate of bootstrap interval coverage probability

bootstrap diagnostic

— Flądrarz
źródło

3

Ludzie często zapominają o innym użyciu bootstrap: w celu zidentyfikowania i skorygowania uprzedzeń . Podejrzewam, że jeśli włączysz korektę błędu w ładowaniu początkowym, możesz uzyskać znacznie lepszą wydajność z CI.

— whuber

@whuber: fajny punkt, +1. O ile pamiętam, metody Bootstrap i ich aplikacje Davisona i Hinkleya stanowią przyjemne i przystępne wprowadzenie do korekcji błędów i innych ulepszeń w bootstrapie.

— S. Kolassa - Przywróć Monikę

1

Warto wypróbować inne warianty bootstrap, zwłaszcza podstawowy bootstrap.

— Frank Harrell,

3

Bootstrapping to procedura z dużą próbką. nie jest duże, szczególnie dla danych logarytmicznych .

n = 10

$n = 10$

— Cliff AB

16

Diagnostyka bootstrap i środki zaradcze autorstwa Canto, Davisona, Hinkleya i Ventury (2006) wydają się logicznym punktem wyjścia. Dyskutują na wiele sposobów, w jaki bootstrap może się zepsuć i - co ważniejsze tutaj - oferują diagnostykę i możliwe rozwiązania:

Wartości odstające
Niepoprawny model ponownego próbkowania
Nonpivotality
Niespójność metody bootstrap

W tej sytuacji nie widzę problemu z 1, 2 i 4. Spójrzmy na 3. Jak zauważa @Ben Ogorek (choć zgadzam się z @Glen_b, że dyskusja o normalności może być czerwonym śledziem), ważność bootstrap zależy od kluczowej statystyki, którą jesteśmy zainteresowani.

Sekcja 4 w Canty i in. sugeruje ponowne próbkowanie w ramach próbek, aby uzyskać pomiar odchylenia i wariancji dla oszacowania parametru w każdej próbce ładowania początkowego . Oto kod do replikacji formuł z p. 15 artykułu:

library(boot)
m <- 10 # sample size
n.boot <- 1000
inner.boot <- 1000

set.seed(1)
samp.mean <- bias <- vars <- rep(NA,n.boot)
for ( ii in 1:n.boot ) {
    samp <- exp(rnorm(m,0,2)) + 1
    samp.mean[ii] <- mean(samp)
    foo <- boot(samp,statistic=function(xx,index)mean(xx[index]),R=inner.boot)
    bias[ii] <- mean(foo$t[,1])-foo$t0
    vars[ii] <- var(foo$t[,1])
}

opar <- par(mfrow=c(1,2))
    plot(samp.mean,bias,xlab="Sample means",ylab="Bias",
        main="Bias against sample means",pch=19,log="x")
    abline(h=0)
    plot(samp.mean,vars,xlab="Sample means",ylab="Variance",
        main="Variance against sample means",pch=19,log="xy")
par(opar)

diagnostyka bootstrap

Zwróć uwagę na skale logów - bez logów jest to jeszcze bardziej rażące. Ładnie widzimy, jak wariancja średniej wartości bootstrap idzie w górę ze średnią próbki bootstrap. To dla mnie wygląda na dość pistoletu do palenia, aby obwinić o brak wzajemności jako winowajcę pokrycia niskiej częstotliwości.

Jednak z radością przyznaję, że można śledzić na wiele sposobów. Na przykład moglibyśmy spojrzeć na to, czy przedział ufności z konkretnej repliki ładowania początkowego zawiera prawdziwą średnią zależy od średniej z konkretnej repliki.

Jeśli chodzi o środki zaradcze, Canty i in. dyskutuj o transformacjach i przychodzą mi na myśl logarytmy (np. bootstrap i buduj przedziały ufności nie dla średniej, ale dla średniej z logowanych danych), ale tak naprawdę nie mogłem tego zrobić.

Canty i in. kontynuuj dyskusję, w jaki sposób można zmniejszyć zarówno liczbę wewnętrznych ładowań początkowych, jak i pozostały hałas, ważąc próbkowanie i wygładzanie, a także dodając przedziały ufności do wykresów przestawnych.

To może być fajna praca dyplomowa dla inteligentnego studenta. Byłbym wdzięczny za wszelkie wskazówki dotyczące tego, gdzie popełniłem błąd, a także do każdej innej literatury. I pozwolę sobie dodać diagnostictag do tego pytania.

— S. Kolassa - Przywróć Monikę
źródło

13

\hat{μ} - μ

$\hat{\mu} - \mu$

\hat{μ}

$\hat{\mu}$

t

$t$

\sqrt{m} \frac{\hat{μ} - μ}{\hat{σ}}

$\sqrt{m} \frac{\hat{\mu} - \mu}{\hat{\sigma}}$

Potem pomyślałem trochę więcej o całej konfiguracji. Czy przy zaledwie 10 obserwacjach i wyjątkowo wypaczonym rozkładzie nie jest w zasadzie niemożliwe, aby nieparametrycznie oszacować średnią, a tym bardziej skonstruować przedziały ufności z odpowiednim zakresem?

$e^2 + 1 = 8.39$ $P(X \leq 2) = 0.84$ $X \sim \mathcal{N}(0,4)$ $0.84$ $0.84^{10} = 0.178$ . Tak więc w nieco mniej niż 18% przypadków największa obserwacja jest mniejsza niż średnia. Aby uzyskać zasięg większy niż 0,82, potrzebujemy konstrukcji przedziału ufności dla średniej wykraczającej poza największą obserwację. Trudno mi wyobrazić sobie, jak można wykonać (i uzasadnić) taką konstrukcję bez wcześniejszych założeń, że rozkład jest bardzo wypaczony. Ale z zadowoleniem przyjmuję wszelkie sugestie.

— NRH
źródło

Zgadzam się z Tobą. Naprawdę chciałem o tym pomyśleć z punktu widzenia kogoś, kto ma jedną próbkę z tej dystrybucji. Skąd mam wiedzieć, że beztroskie korzystanie z bootstrap w tym przypadku jest niebezpieczne? Jedyną rzeczą, o której mogę myśleć, jest to, że mogłem wziąć dzienniki przed wykonaniem analizy, ale jeden z pozostałych osób odpowiadających twierdzi, że to naprawdę nie pomaga.

— Flądrowiec

1

Nie będziesz wiedział, czy jest to bezpieczne, czy tylko z 10 punktów danych. Jeśli podejrzewasz skośność lub ciężkie ogony, rozwiązaniem może być skupienie się na innym parametrze niż średnia. Na przykład średnia logarytmiczna lub mediana. Nie da to oszacowania (lub przedziału ufności dla) średniej, chyba że podejmiesz dodatkowe założenia, ale lepszym pomysłem może być skupienie się na parametrze, który jest mniej wrażliwy na ogony rozkładu.

— NRH,

6

Obliczenia były prawidłowe, sprawdziłem krzyżowo przy użyciu dobrze znanego rozruchu pakietu . Dodatkowo dodałem przedział BCa (Efron), poprawioną przez błąd wersji procentowego przedziału ładowania początkowego:

for (i in 1:1000) {
  samp <- exp(rnorm(m, 0, 2)) + 1

  boot.out <- boot(samp, function(d, i) sum(d[i]) / m, R=999)
  ci <- boot.ci(boot.out, 0.95, type="all")

  ##tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)
  tCI <- ci$normal[2:3]
      percCI <- ci$perc[4:5]
  bcaCI <- ci$bca[4:5]
      boottCI <- ci$student[4:5]

  if (true.mean > min(tCI) && true.mean < max(tCI)) tCI.total <- tCI.total + 1
  if (true.mean > min(percCI) && true.mean < max(percCI)) percCI.total <- percCI.total + 1 
  if (true.mean > min(bcaCI) && true.mean < max(bcaCI)) bcaCI.total <- bcaCI.total + 1
}

tCI.total/1000     # estimate of t interval coverage probability
0.53
percCI.total/1000  # estimate of percentile interval coverage probability
0.55
bcaCI.total/1000  # estimate of BCa interval coverage probability
0.61

Zakładam, że interwały byłyby znacznie lepsze, gdyby pierwotna wielkość próbki była większa niż 10, powiedzmy 20 lub 50.

Ponadto metoda bootstrap-t zwykle prowadzi do lepszych wyników dla przekrzywionych statystyk. Jednak potrzebuje zagnieżdżonej pętli, a zatem ponad 20 razy więcej czasu obliczeniowego.

Do testowania hipotez bardzo ważne jest również, aby jednostronne pokrycia były dobre. Dlatego patrzenie tylko na dwustronne relacje może często wprowadzać w błąd.

— lambruscoAcido
źródło

1

n < 100

$n<100$

5

Też byłem zdezorientowany i spędziłem dużo czasu na papierowych przedziałach ufności Bootstrap DiCiccio z 1996 roku i Efron , bez wiele do pokazania.

Doprowadziło mnie to do mniejszego myślenia o bootstrapie jako metodzie ogólnego zastosowania. Kiedyś myślałem o tym jak o czymś, co wyciągnęłoby cię z dżemu, kiedy naprawdę utknąłeś. Ale nauczyłem się jego brudnej małej tajemnicy: przedziały ufności ładowania początkowego oparte są w jakiś sposób na normalności. Pozwól mi wyjaśnić.

x \sim N. (μ, σ^{2)})

$x \sim N(\mu, \sigma^2)$

σ

$\sigma$

z = \frac{x - μ}{σ} \sim N. (0, 1)

$z = \frac{x - \mu}{\sigma} \sim N(0,1)$

μ

$\mu$

Pr (- 1.96 \leq \frac{x - μ}{σ} \leq 1.96) = 0.95

$\Pr(-1.96 \le \frac{x - \mu}{\sigma} \le 1.96) = 0.95$

Kiedy myślisz o tym, co uzasadnia percentyle rozkładu normalnego związane z przedziałami ufności, jest ono całkowicie oparte na tej dogodnej kluczowej wielkości. W przypadku arbitralnego rozkładu nie ma teoretycznego związku między percentylami rozkładu próbkowania i przedziałami ufności , a przyjęcie surowych proporcji rozkładu próbkowania bootstrap go nie odcina.

Dlatego interwały BCa (skorygowane odchylenie) Efrona wykorzystują transformacje, aby zbliżyć się do normalności, a metody bootstrap-t polegają na tym, że uzyskane statystyki t są w przybliżeniu kluczowe. Teraz bootstrap potrafi oszacować piekło z chwil i zawsze możesz założyć normalność i użyć standardowego +/- 2 * SE. Ale biorąc pod uwagę całą pracę związaną z wprowadzaniem nieparametrycznej wersji z bootstrapem, nie wydaje się to całkiem uczciwe, prawda?

— Ben Ogorek
źródło

2

Możliwe, że coś przeoczyłem, ale fakt, że ładowanie jest powiązane z kluczowymi lub prawie kluczowymi wielkościami, sam w sobie nie implikuje żadnego związku z normalnością. Ilości kluczowe mogą mieć różne rozkłady w szczególnych okolicznościach. Nie rozumiem też, jak następuje zdanie zapisane kursywą w ostatnim akapicie.

— Glen_b

1

Jak zatem wynika z twierdzenia dotyczącego normalności?

— Glen_b

1

F

$F$

Φ^{- 1} [F (X)]

$\Phi^{-1}[F(X)]$

2

F

$F$

2

Aby dodać do @Glen_b: transformacja do rozkładu normalnego musi istnieć tylko w celu udowodnienia poprawności metody. Nie musisz go znaleźć, aby użyć tej metody. Dodatkowo, jeśli nie lubisz normalnych dystrybucji, możesz przepisać cały dowód za pomocą innej symetrycznej, ciągłej dystrybucji. Zastosowanie normalnych rozkładów jest technicznie użyteczne, ale nie jest absolutnie konieczne, nie mówi nic o źródle danych ani średniej próbki.

— Peter

0

Sprawdź artykuł Tima Hesterberga w The American Statistician na http://www.timhesterberg.net/bootstrap#TOC-What-Teachers-Should-Know-about-the-Bootstrap:-Resampling-in-the-Undergraduate-Statistics-Curriculum .

Zasadniczo interwał centylowy ładowania początkowego nie ma dużego prawdopodobieństwa pokrycia dla wypaczonych danych, chyba że n jest duże.

— Gość
źródło