# Example 9.8 of Davison and Hinkley (1997) requires tilting the resampling
# distribution of the studentized statistic to be centred at the observed
# value of the test statistic 1.84. In the book exponential tilting was used
# but it is also possible to use smooth.f.
data(gravity)
grav1 <- gravity[as.numeric(gravity[,2])>=7,]
grav.fun <- function(dat, w, orig)
{ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
d <- dat[, 1]
ns <- tabulate(strata)
w <- w/tapply(w, strata, sum)[strata]
mns <- tapply(d * w, strata, sum)
mn2 <- tapply(d * d * w, strata, sum)
s2hat <- sum((mn2 - mns^2)/ns)
c(mns[2]-mns[1],s2hat,(mns[2]-mns[1]-orig)/sqrt(s2hat),)
}
grav.z0 <- grav.fun(grav1,rep(1,26),0)
grav.boot <- boot(grav1, grav.fun, R=499, stype="w",
strata=grav1[,2], orig=grav.z0[1])
grav.sm <- smooth.f(grav.z0[3], grav.boot, index=3)
# Now we can run another bootstrap using these weights
grav.boot2 <- boot(grav1, grav.fun, R=499, stype="w",
strata=grav1[,2], orig=grav.z0[1],
weights=grav.sm)
# Estimated p-values can be found from these as follows
mean(grav.boot$t[,3] >= grav.z0[3])
imp.prob(grav.boot2,t0=-grav.z0[3],t=-grav.boot2$t[,3])
# Note that for the importance sampling probability we must
# multiply everything by -1 to ensure that we find the correct
# probability. Raw resampling is not reliable for probabilities
# greater than 0.5. Thus
1-imp.prob(grav.boot2,index=3,t0=grav.z0[3])$raw
# can give very strange results (negative probabilities).
Run the code above in your browser using DataLab