Warum hat mein Bootstrap-Intervall eine schreckliche Abdeckung?

29

Ich wollte eine Klassendemonstration durchführen, bei der ich ein t-Intervall mit einem Bootstrap-Intervall vergleiche und die Überdeckungswahrscheinlichkeit für beide berechne. Ich wollte, dass die Daten aus einer verzerrten Verteilung stammen, also habe ich mich dafür entschieden, die Daten als exp(rnorm(10, 0, 2)) + 1eine Stichprobe der Größe 10 aus einem verschobenen Lognormal zu generieren . Ich habe ein Skript zum Zeichnen von 1000 Samples geschrieben und für jedes Sample ein 95% t-Intervall und ein 95% Bootstrap-Perzentilintervall basierend auf 1000 Replikaten berechnet.

Wenn ich das Skript ausführe, geben beide Methoden sehr ähnliche Intervalle an und beide haben eine Abdeckungswahrscheinlichkeit von 50-60%. Ich war überrascht, weil ich dachte, das Bootstrap-Intervall wäre besser.

Meine Frage ist, habe ich

Fehler im Code gemacht?
Fehler bei der Berechnung der Intervalle gemacht?
einen Fehler gemacht, indem erwartet wurde, dass das Bootstrap-Intervall bessere Abdeckungseigenschaften hat?

Gibt es in dieser Situation auch eine Möglichkeit, ein zuverlässigeres CI zu erstellen?

 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 Flunder
quelle

3

Die Leute vergessen oft eine andere Verwendung des Bootstraps: das Erkennen und Korrigieren von Voreingenommenheit . Ich vermute, wenn Sie eine Bias-Korrektur in Ihr Bootstrapping einbeziehen, könnte das CI eine viel bessere Leistung bringen.

Whuber

@whuber: schöner Punkt, +1. Soweit ich mich erinnere, bieten Bootstrap-Methoden und ihre Anwendungen von Davison & Hinkley eine schöne und leicht zugängliche Einführung in die Bias-Korrektur und andere Verbesserungen des Bootstraps.

S. Kolassa - Wiedereinsetzung von Monica am

1

Es lohnt sich, die anderen Bootstrap-Varianten auszuprobieren, insbesondere den Basis-Bootstrap.

Frank Harrell

3

Bootstrapping ist ein Verfahren mit umfangreichen Beispielen. ist nicht groß, insbesondere für lognormale Daten .

n = 10

$n = 10$

Cliff AB

16

Bootstrap-Diagnosen und Abhilfemaßnahmen von Canto, Davison, Hinkley & Ventura (2006) scheinen ein logischer Ausgangspunkt zu sein. Sie diskutieren verschiedene Möglichkeiten, wie der Bootstrap ausfallen kann, und bieten - was noch wichtiger ist - Diagnosen und mögliche Abhilfemaßnahmen an:

Ausreißer
Falsches Resampling-Modell
Nichtpivotalität
Inkonsistenz der Bootstrap-Methode

Ich sehe in dieser Situation kein Problem mit 1, 2 und 4. Schauen wir uns 3 an. Wie @Ben Ogorek feststellt (obwohl ich mit @Glen_b übereinstimme, dass die Normalitätsdiskussion ein roter Hering sein könnte), hängt die Gültigkeit des Bootstraps von der Dreh- und Angelpunktzahl der Statistik ab, an der wir interessiert sind.

Abschnitt 4 in Canty et al. schlägt Resampling-in-Resamples vor, um ein Maß für die Abweichung und Varianz für die Parameterschätzung in jedem Bootstrap-Resample zu erhalten . Hier ist Code, um die Formeln von p zu replizieren. 15 des Artikels:

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)

Bootstrap-Diagnose

Beachten Sie die Protokollskalen - ohne Protokolle ist dies noch offensichtlicher. Wir sehen gut, wie die Varianz der Bootstrap-Mittelwertschätzung mit dem Mittelwert der Bootstrap-Stichprobe steigt. Das scheint mir eine rauchende Waffe genug zu sein, um der Nicht-Pivotalität als Schuld an der Abdeckung des niedrigen Vertrauensintervalls die Schuld zu geben.

Ich gebe jedoch gerne zu, dass man auf viele Arten weiterverfolgen kann. Beispielsweise könnten wir untersuchen, ob das Konfidenzintervall eines bestimmten Bootstrap-Replikats den wahren Mittelwert enthält, der vom Mittelwert des jeweiligen Replikats abhängt.

In Bezug auf Abhilfemaßnahmen haben Canty et al. Hier werden Transformationen besprochen und Logarithmen in den Sinn kommen (z. B. Bootstrap und Erstellen von Konfidenzintervallen nicht für den Mittelwert, sondern für den Mittelwert der protokollierten Daten), aber ich konnte es nicht wirklich zum Laufen bringen.

Canty et al. Fahren Sie mit der Diskussion fort, wie Sie sowohl die Anzahl der inneren Bootstraps als auch das verbleibende Rauschen durch wichtiges Abtasten und Glätten reduzieren und den Pivot-Plots Vertrauensbereiche hinzufügen können.

Dies könnte ein unterhaltsames Diplomarbeitsprojekt für einen klugen Studenten sein. Ich würde mich über Hinweise freuen, auf die ich mich geirrt habe, sowie auf jede andere Literatur. Und ich erlaube mir, das diagnosticTag zu dieser Frage hinzuzufügen .

S. Kolassa - Setzen Sie Monica wieder ein
quelle

13

\hat{μ} - μ

$\hat{\mu} - \mu$

\hat{μ}

$\hat{\mu}$

t

$t$

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

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

Dann habe ich ein bisschen mehr über das ganze Setup nachgedacht. Ist es dann mit nur 10 Beobachtungen und einer extrem verzerrten Verteilung nicht grundsätzlich unmöglich, den Mittelwert nichtparametrisch zu schätzen, geschweige denn Konfidenzintervalle mit der richtigen Abdeckung zu konstruieren?

$e^2 + 1 = 8.39$ $P(X \leq 2) = 0.84$ $X \sim \mathcal{N}(0,4)$ $0.84$ $0.84^{10} = 0.178$ . In etwas weniger als 18% der Fälle ist die größte Beobachtung kleiner als der Mittelwert. Um eine Abdeckung von mehr als 0,82 zu erhalten, benötigen wir die Konstruktion eines Konfidenzintervalls für den Mittelwert, der über die größte Beobachtung hinausgeht. Es fällt mir schwer, mir vorzustellen, wie eine solche Konstruktion hergestellt (und gerechtfertigt) werden kann, ohne vorher anzunehmen, dass die Verteilung extrem schief ist. Aber ich freue mich über Vorschläge.

NRH
quelle

Ich stimme mit Ihnen ein. Ich wollte wirklich darüber nachdenken, aus der Sicht von jemandem, der eine Stichprobe aus dieser Distribution hat. Woher weiß ich, dass es in diesem Fall unsicher ist, den Bootstrap munter zu benutzen? Der einzige Gedanke, an den ich denken kann, ist, dass ich möglicherweise Protokolle erstellt habe, bevor ich die Analyse durchgeführt habe, aber einer der anderen Befragten sagt, dass dies nicht wirklich hilfreich ist.

Flunder

1

Sie werden anhand der 10 Datenpunkte allein nicht wissen, ob es sicher oder unsicher ist. Wenn Sie vermuten, dass es sich um eine Schiefe oder einen starken Schwanz handelt, kann die Lösung darin bestehen, sich auf einen anderen Parameter als den Mittelwert zu konzentrieren. Zum Beispiel der log-mean oder der Median. Auf diese Weise erhalten Sie keine Schätzung (oder kein Konfidenzintervall) für den Mittelwert, es sei denn, Sie treffen zusätzliche Annahmen. Es ist jedoch möglicherweise eine bessere Idee, sich auf einen Parameter zu konzentrieren, der weniger empfindlich für die Schwänze der Verteilung ist.

NRH

6

Die Berechnungen stimmten, ich habe mit dem bekannten Paket- Boot abgeglichen . Zusätzlich habe ich das BCa-Intervall (von Efron) hinzugefügt, eine verzerrte Version des Perzentil-Bootstrap-Intervalls:

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

Ich gehe davon aus, dass die Intervalle viel besser wären, wenn die ursprüngliche Stichprobengröße größer als 10 wäre, z. B. 20 oder 50.

Darüber hinaus führt die Bootstrap-t- Methode in der Regel zu besseren Ergebnissen für verzerrte Statistiken. Es benötigt jedoch eine verschachtelte Schleife und daher mehr als 20 Mal mehr Rechenzeit.

Für das Testen von Hypothesen ist es auch sehr wichtig, dass die einseitige Abdeckung gut ist. Daher kann es oft irreführend sein, nur die doppelseitigen Abdeckungen zu betrachten.

lambruscoAcido
quelle

1

n < 100

$n<100$

5

Ich war auch darüber verwirrt und verbrachte viel Zeit mit den Bootstrap Confidence Intervals von 1996 von DiCiccio und Efron , ohne viel dafür zu beweisen .

Es hat mich tatsächlich dazu gebracht, weniger an den Bootstrap als eine Allzweckmethode zu denken. Früher stellte ich es mir als etwas vor, das Sie aus einem Stau zog, wenn Sie wirklich feststeckten. Aber ich habe sein schmutziges kleines Geheimnis gelernt: Die Bootstrap-Konfidenzintervalle basieren alle auf der einen oder anderen Normalität. Lassen Sie mich das erklären.

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$

Wenn Sie darüber nachdenken, was die Beziehung der Perzentile der Normalverteilung zu Konfidenzintervallen rechtfertigt, dann basiert sie vollständig auf dieser geeigneten zentralen Größe. Bei einer beliebigen Verteilung gibt es keinen theoretischen Zusammenhang zwischen den Perzentilen der Stichprobenverteilung und den Konfidenzintervallen , und wenn rohe Anteile der Bootstrap-Stichprobenverteilung verwendet werden, wird dies nicht beeinträchtigt.

Daher verwenden Efrons BCa-Intervalle (Bias-korrigiert) Transformationen, um die Normalität zu approximieren, und Bootstrap-t-Methoden setzen voraus, dass die resultierenden t-Statistiken ungefähr von zentraler Bedeutung sind. Jetzt kann der Bootstrap die Hölle von Momenten abschätzen, und Sie können immer von Normalität ausgehen und den Standard +/- 2 * SE verwenden. Aber in Anbetracht der ganzen Arbeit, die mit dem Bootstrap anfing, nicht parametrisch zu werden, scheint es nicht ganz fair zu sein, oder?

Ben Ogorek
quelle

2

Es ist möglich, dass ich etwas verpasst habe, aber die Tatsache, dass Bootstrapping mit schwenkbaren oder nahezu schwenkbaren Größen verbunden ist, impliziert an sich keine Assoziation mit Normalität. Gelenkmengen können unter bestimmten Umständen alle Arten von Verteilungen aufweisen. Ich sehe auch nicht, wie der kursive Satz in Ihrem vorletzten Absatz folgt.

Glen_b -Reinstate Monica

1

Wie folgt dann die Behauptung zur Normalität?

Glen_b -Reinstate Monica

1

F

$F$

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

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

2

F

$F$

2

Hinzufügen zu @Glen_b: Die Umwandlung in eine Normalverteilung muss nur vorhanden sein, um die Richtigkeit der Methode zu beweisen. Sie müssen es nicht finden, um die Methode zu verwenden. Wenn Sie Normalverteilungen nicht mögen, können Sie den gesamten Proof mit einer anderen symmetrischen, kontinuierlichen Verteilung umschreiben. Die Verwendung von Normalverteilungen ist technisch sinnvoll, aber nicht unbedingt erforderlich. Sie sagt nichts über die Datenquelle oder den Stichprobenmittelwert aus.

Peter

0

Lesen Sie Tim Hesterbergs Artikel in The American Statistician unter http://www.timhesterberg.net/bootstrap#TOC-What-Teachers-Should-Know-about-the-Bootstrap:-Resampling-in-the-Undergraduate-Statistics-Curriculum .

Im Wesentlichen hat das Bootstrap-Perzentilintervall keine hohe Abdeckungswahrscheinlichkeit für verzerrte Daten, es sei denn, n ist groß.

Gast
quelle

Warum hat mein Bootstrap-Intervall eine schreckliche Abdeckung?

Antworten: