Prawdopodobieństwa pokrycia podstawowego przedziału ufności ładowania początkowego

Mam następujące pytanie do kursu, nad którym pracuję:

Przeprowadź badanie Monte Carlo, aby oszacować prawdopodobieństwo pokrycia standardowego normalnego przedziału ufności bootstrapu i podstawowego przedziału ufności bootstrapu. Próbka z normalnej populacji i sprawdź empiryczne wskaźniki pokrycia dla średniej próby.

Prawdopodobieństwa pokrycia dla standardowego normalnego CI bootstrap są łatwe:

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

Z tego, czego nauczyłem się na tym kursie, podstawowy przedział ufności ładowania początkowego można obliczyć w następujący sposób:

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

To ma sens. Nie rozumiem, jak obliczyć prawdopodobieństwo pokrycia dla podstawowego CI bootstrap. Rozumiem, że prawdopodobieństwo pokrycia reprezentowałoby, ile razy CI zawiera prawdziwą wartość (w tym przypadku mu). Czy po prostu uruchamiam tę bootfunkcję wiele razy?

Jak mogę podejść do tego pytania inaczej?

Terminologia prawdopodobnie nie jest stosowana konsekwentnie, dlatego poniżej rozumiem tylko oryginalne pytanie. Z mojego rozumienia, obliczone przez ciebie normalne CI nie są tym, o co cię proszono. Każdy zestaw replik ładowania początkowego daje jeden przedział ufności, nie wiele. Sposób obliczania różnych typów CI na podstawie wyników zestawu replik ładowania początkowego jest następujący:

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

Ponieważ chcę porównać obliczenia z wynikami z pakietu boot, najpierw definiuję funkcję, która będzie wywoływana dla każdej replikacji. Jego argumentami są oryginalna próbka oraz wektor indeksu określający przypadki dla pojedynczej repliki. Zwraca , oszacowanie wtyczki dla , a także , oszacowanie wtyczki dla wariancji średniej . Ten ostatni będzie wymagany tylko w przypadku bootstrap -CI. μ S 2 ⋆ M σ 2 M$M^{\star}$$\mu$$S_{M}^{2\star}$$\sigma_{M}^{2}$$t$

> 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

Bez użycia pakietu bootmożesz po prostu użyć, replicate()aby uzyskać zestaw replik ładowania początkowego.

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

Ale trzymajmy się wyników z, boot.ci()aby mieć referencję.

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

Podstawowy, percentyl i CI polegają na empirycznym rozkładzie oszacowań bootstrap. Aby uzyskać kwantyle i , znajdujemy odpowiadające indeksy posortowanemu wektorowi oszacowań bootstrap (zauważ, że wykona bardziej skomplikowaną interpolację w celu znalezienia kwantyli empirycznych, gdy indeksy nie są liczbami naturalnymi) .α / 2 1 - α /$t$$\alpha/2$$1 - \alpha/2$boot.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

Dla -CI potrzebujemy oszacowań bootstrap, aby obliczyć krytyczne wartości . Dla standardowego normalnego CI, wartość krytyczna będzie po prostu być -value od rozkładu normalnego.t ⋆ t $t$$t^{\star}$$t$$z$

# 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

Aby oszacować prawdopodobieństwo pokrycia tych typów CI, będziesz musiał uruchomić tę symulację wiele razy. Po prostu zawiń kod do funkcji, zwróć listę z wynikami CI i uruchom go z replicate()podobną demonstracją w tej liście .