Kolejny wariant, który jest nieco uproszczony, ale myślę, że dostarczam wiadomość bez jawnego korzystania z biblioteki, boot
która może mylić niektóre osoby z używaną składnią.
Mamy model liniowy: ,y= Xβ+ ϵε ~ N( 0 ,σ2))
Poniżej przedstawiono parametryczny bootstrap dla tego modelu liniowego, co oznacza, że nie próbkujemy ponownie naszych oryginalnych danych, ale w rzeczywistości generujemy nowe dane z naszego dopasowanego modelu. Dodatkowo zakładamy, że rozkład początkowy współczynnika regresji jest symetryczny, a więc niezmienny w tłumaczeniu. (Mówiąc z grubsza, że możemy przesuwać jego oś, wpływając na jej właściwości) Pomysł polega na tym, że wahania są spowodowane a zatem przy wystarczającej liczbie próbek powinny one zapewnić dobre przybliżenie prawdziwego rozkładu z . Tak jak poprzednio, ponownie testujemy i zdefiniowaliśmy nasze wartości p jakoββϵβH.0: 0 =βjot„prawdopodobieństwo, biorąc pod uwagę zerową hipotezę rozkładu prawdopodobieństwa danych, że wynik byłby tak ekstremalny, jak bardziej ekstremalny niż obserwowany wynik” (gdzie obserwowane wyniki w tym przypadku są, które otrzymaliśmy dla naszego oryginalnego modelu). Więc oto idzie:β
# Sample Size
N <- 2^12;
# Linear Model to Boostrap
Model2Boot <- lm( mpg ~ wt + disp, mtcars)
# Values of the model coefficients
Betas <- coefficients(Model2Boot)
# Number of coefficents to test against
M <- length(Betas)
# Matrix of M columns to hold Bootstraping results
BtStrpRes <- matrix( rep(0,M*N), ncol=M)
for (i in 1:N) {
# Simulate data N times from the model we assume be true
# and save the resulting coefficient in the i-th row of BtStrpRes
BtStrpRes[i,] <-coefficients(lm(unlist(simulate(Model2Boot)) ~wt + disp, mtcars))
}
#Get the p-values for coefficient
P_val1 <-mean( abs(BtStrpRes[,1] - mean(BtStrpRes[,1]) )> abs( Betas[1]))
P_val2 <-mean( abs(BtStrpRes[,2] - mean(BtStrpRes[,2]) )> abs( Betas[2]))
P_val3 <-mean( abs(BtStrpRes[,3] - mean(BtStrpRes[,3]) )> abs( Betas[3]))
#and some parametric bootstrap confidence intervals (2.5%, 97.5%)
ConfInt1 <- quantile(BtStrpRes[,1], c(.025, 0.975))
ConfInt2 <- quantile(BtStrpRes[,2], c(.025, 0.975))
ConfInt3 <- quantile(BtStrpRes[,3], c(.025, 0.975))
Jak wspomniano, cały pomysł polega na tym, że masz bootstrapped dystrybucję zbliżoną do ich prawdziwej. (Oczywiście ten kod jest zoptymalizowany pod kątem szybkości, ale pod kątem czytelności. :))β