Abdeckungswahrscheinlichkeiten des grundlegenden Bootstrap-Konfidenzintervalls

11

Ich habe die folgende Frage für einen Kurs, an dem ich arbeite:

Führen Sie eine Monte-Carlo-Studie durch, um die Abdeckungswahrscheinlichkeiten des normalen Standard-Bootstrap-Konfidenzintervalls und des grundlegenden Bootstrap-Konfidenzintervalls abzuschätzen. Stichprobe aus einer normalen Population und Überprüfung der empirischen Abdeckungsraten für den Stichprobenmittelwert.

Die Abdeckungswahrscheinlichkeiten für das normale Standard-Bootstrap-CI sind einfach:

n = 1000;
alpha = c(0.025, 0.975);
x = rnorm(n, 0, 1);
mu = mean(x);
sqrt.n = sqrt(n);

LNorm = numeric(B);
UNorm = numeric(B);

for(j in 1:B)
{
    smpl = x[sample(1:n, size = n, replace = TRUE)];
    xbar = mean(smpl);
    s = sd(smpl);

    LNorm[j] = xbar + qnorm(alpha[1]) * (s / sqrt.n);
    UNorm[j] = xbar + qnorm(alpha[2]) * (s / sqrt.n);
}

mean(LNorm < 0 & UNorm > 0); # Approximates to 0.95
# NOTE: it is not good enough to look at overall coverage
# Must compute separately for each tail

Aus dem, was ich für diesen Kurs gelernt habe, kann das grundlegende Bootstrap- Konfidenzintervall wie folgt berechnet werden:

# Using x from previous...
R = boot(data = x, R=1000, statistic = function(x, i){ mean(x[i]); });
result = 2 * mu - quantile(R$t, alpha, type=1);

Das macht Sinn. Was ich nicht verstehe, ist die Berechnung der Abdeckungswahrscheinlichkeiten für das grundlegende Bootstrap-CI. Ich verstehe, dass die Abdeckungswahrscheinlichkeit die Häufigkeit darstellt, mit der das CI den wahren Wert enthält (in diesem Fall mu). Führe ich die bootFunktion einfach viele Male aus?

Wie kann ich diese Frage anders angehen?

TheCloudlessSky
quelle
Ist dein size=100Tippfehler? Ich glaube nicht, dass Sie die richtigen oberen und unteren Grenzen erhalten, da die implizite Stichprobengröße 1000 zu sein scheint, wenn Sie Ihre CIs in der Schleife berechnen (da Sie sie sqrt.nin der Berechnung verwenden). Warum vergleichen Sie mit muund nicht direkt mit 0 (letzteres ist der wahre Mittelwert)?
Kardinal
Auch smpl = x[sample(1:n, size = 100, replace = TRUE)]; kann vereinfacht werden , um smpl = sample(x, size=100, replace=TRUE).
Kardinal
@cardinal - Ja, es war ein Tippfehler und das Gleiche muwie 0. Das normale CI funktioniert einwandfrei. Es ist das grundlegende Bootstrap-CI, mit dem ich Schwierigkeiten habe.
TheCloudlessSky

Antworten:

16

Die Terminologie wird wahrscheinlich nicht konsistent verwendet, daher verstehe ich im Folgenden nur die ursprüngliche Frage. Nach meinem Verständnis sind die von Ihnen berechneten normalen CIs nicht das, wonach gefragt wurde. Jeder Satz von Bootstrap-Replikaten gibt Ihnen ein Konfidenzintervall, nicht viele. Die Methode zum Berechnen verschiedener CI-Typen aus den Ergebnissen einer Reihe von Bootstrap-Replikaten lautet wie folgt:

B    <- 999                  # number of replicates
muH0 <- 100                  # for generating data: true mean
sdH0 <- 40                   # for generating data: true sd
N    <- 200                  # sample size
DV   <- rnorm(N, muH0, sdH0) # simulated data: original sample

bootM.μS.M.2σM.2t

> getM <- function(orgDV, idx) {
+     bsM   <- mean(orgDV[idx])                       # M*
+     bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N    # S^2*(M)
+     c(bsM, bsS2M)
+ }

> library(boot)                                       # for boot(), boot.ci()
> bOut <- boot(DV, statistic=getM, R=B)
> boot.ci(bOut, conf=0.95, type=c("basic", "perc", "norm", "stud"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL : 
boot.ci(boot.out = bOut, conf = 0.95, type = c("basic", "perc", "norm", "stud"))

Intervals : 
Level      Normal            Basic         Studentized        Percentile    
95%   ( 95.6, 106.0 )   ( 95.7, 106.2 )  ( 95.4, 106.2 )   ( 95.4, 106.0 )  
Calculations and Intervals on Original Scale

Ohne Verwendung eines Pakets können bootSie einfach replicate()eine Reihe von Bootstrap-Replikaten abrufen.

boots <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))

Aber bleiben wir bei den Ergebnissen von boot.ci(), um eine Referenz zu haben.

boots   <- bOut$t                     # estimates from all replicates
M       <- mean(DV)                   # M from original sample
S2M     <- (((N-1)/N) * var(DV)) / N  # S^2(M) from original sample
Mstar   <- boots[ , 1]                # M* for each replicate
S2Mstar <- boots[ , 2]                # S^2*(M) for each replicate
biasM   <- mean(Mstar) - M            # bias of estimator M

tα/.21- -α/.2boot.ci()

(idx <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2)) # indices for sorted vector of estimates
[1] 25 975

> (ciBasic <- 2*M - sort(Mstar)[idx])          # basic CI
[1] 106.21826  95.65911

> (ciPerc <- sort(Mstar)[idx])                 # percentile CI
[1] 95.42188 105.98103

tttz

# standard normal CI with bias correction
> zCrit   <- qnorm(c(0.025, 0.975))   # z-quantiles from std-normal distribution
> (ciNorm <- M - biasM + zCrit * sqrt(var(Mstar)))
[1] 95.5566 106.0043

> tStar <- (Mstar-M) / sqrt(S2Mstar)  # t*
> tCrit <- sort(tStar)[idx]           # t-quantiles from empirical t* distribution
> (ciT  <- M - tCrit * sqrt(S2M))     # studentized t-CI
[1] 106.20690  95.44878

Um die Abdeckungswahrscheinlichkeiten dieser CI-Typen abzuschätzen, müssen Sie diese Simulation viele Male ausführen. Wickeln Sie einfach den Code in eine Funktion ein, geben Sie eine Liste mit den CI-Ergebnissen zurück und führen Sie sie replicate()wie in dieser Übersicht gezeigt aus .

Karakal
quelle
Beeindruckend! - Tolle Erklärung, was ich falsch gemacht habe. Auch - danke für die Code-Tipps! Das funktioniert perfekt!
TheCloudlessSky
Ok, eine letzte Frage: Wenn ich versuche, diese Informationen zu replizieren, habe ich eine Funktion erstellt computeCIsund aufgerufen results = replicate(500, computeCIs());. Am Ende computeCIskehrt es zurück c(ciBasic, ciPerc). Sollte ich dann nicht testen mean(results[1, ] < 0 & results[2, ] > 0), um alle Basis-CIs zu testen, die den wahren Mittelwert (die Deckungswahrscheinlichkeit) enthalten, um die Abdeckungswahrscheinlichkeiten zu testen? Wenn ich das mache, bekomme ich, 1wenn ich denke, ich sollte es bekommen 0.95.
TheCloudlessSky
@ TheCloudlessSky Für die vollständige Funktion und vollständige Simulation mit den erwarteten Ergebnissen in Bezug auf die Abdeckungsfrequenzen, siehe pastebin.com/qKpNKK0D
caracal
Ja, ich bin ein Idiot :) ... Ich habe beim Kopieren des Codes in R einen Tippfehler gemacht ... danke für all Ihre Hilfe! :)
TheCloudlessSky
Danke @caracal für die nette Antwort. Die Verbindung pastebin.com/qKpNKK0Dist unterbrochen. Würde mich freuen, wenn Sie es aktualisieren und die vollständige Funktion und vollständige Simulation bereitstellen. Danke
MYaseen208