The R Project SVN
R-packages
(root)/trunk/boot/tests/Examples/boot-Ex.Rout.save – Rev 7321
Rev 6900 | Rev 7646 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2017-02-10 r72157) -- "Unsuffered Consequences"Copyright (C) 2017 The R Foundation for Statistical ComputingPlatform: x86_64-pc-linux-gnu (64-bit)R is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.Natural language support but running in an English localeR is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> pkgname <- "boot"> source(file.path(R.home("share"), "R", "examples-header.R"))> options(warn = 1)> library('boot')>> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')> cleanEx()> nameEx("Imp.Estimates")> ### * Imp.Estimates>> flush(stderr()); flush(stdout())>> ### Name: Imp.Estimates> ### Title: Importance Sampling Estimates> ### Aliases: Imp.Estimates imp.moments imp.prob imp.quantile imp.reg> ### Keywords: htest nonparametric>> ### ** Examples>> # 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 this example> # we show how certain estimates can be found using resamples taken from> # the tilted distribution.> 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 <- as.vector(tapply(d * w, strata, sum)) # drop names+ 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.L <- empinf(data = grav1, statistic = grav.fun, stype = "w",+ strata = grav1[,2], index = 3, orig = grav.z0[1])> grav.tilt <- exp.tilt(grav.L, grav.z0[3], strata = grav1[, 2])> grav.tilt.boot <- boot(grav1, grav.fun, R = 199, stype = "w",+ strata = grav1[, 2], weights = grav.tilt$p,+ orig = grav.z0[1])> # Since the weights are needed for all calculations, we shall calculate> # them once only.> grav.w <- imp.weights(grav.tilt.boot)> grav.mom <- imp.moments(grav.tilt.boot, w = grav.w, index = 3)> grav.p <- imp.prob(grav.tilt.boot, w = grav.w, index = 3, t0 = grav.z0[3])> unlist(grav.p)t0 raw rat reg1.8401182 1.0222447 0.9778246 0.9767292> grav.q <- imp.quantile(grav.tilt.boot, w = grav.w, index = 3,+ alpha = c(0.9, 0.95, 0.975, 0.99))> as.data.frame(grav.q)alpha raw rat reg1 0.900 3.048484 1.170707 5.0560042 0.950 3.237935 1.565928 5.0560043 0.975 3.629448 1.895876 5.0560044 0.990 3.629448 2.258157 5.056004>>>> cleanEx()> nameEx("abc.ci")> ### * abc.ci>> flush(stderr()); flush(stdout())>> ### Name: abc.ci> ### Title: Nonparametric ABC Confidence Intervals> ### Aliases: abc.ci> ### Keywords: nonparametric htest>> ### ** Examples>> ## Don't show:> op <- options(digits = 5)> ## End(Don't show)> # 90% and 95% confidence intervals for the correlation> # coefficient between the columns of the bigcity data>> abc.ci(bigcity, corr, conf=c(0.90,0.95))conf[1,] 0.90 0.95815 0.99173[2,] 0.95 0.94937 0.99307>> # A 95% confidence interval for the difference between the means of> # the last two samples in gravity> mean.diff <- function(y, w)+ { gp1 <- 1:table(as.numeric(y$series))[1]+ sum(y[gp1, 1] * w[gp1]) - sum(y[-gp1, 1] * w[-gp1])+ }> grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]> abc.ci(grav1, mean.diff, strata = grav1$series)[1] 0.95000 -6.70758 -0.39394> ## Don't show:> options(op)> ## End(Don't show)>>>> cleanEx()> nameEx("boot")> ### * boot>> flush(stderr()); flush(stdout())>> ### Name: boot> ### Title: Bootstrap Resampling> ### Aliases: boot boot.return c.boot> ### Keywords: nonparametric htest>> ### ** Examples>> ## Don't show:> op <- options(digits = 5)> ## End(Don't show)> # Usual bootstrap of the ratio of means using the city data> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)> boot(city, ratio, R = 999, stype = "w")ORDINARY NONPARAMETRIC BOOTSTRAPCall:boot(data = city, statistic = ratio, R = 999, stype = "w")Bootstrap Statistics :original bias std. errort1* 1.5203 0.0396 0.21671>>> # Stratified resampling for the difference of means. In this> # example we will look at the difference of means between the final> # two series in the gravity data.> diff.means <- function(d, f)+ { n <- nrow(d)+ gp1 <- 1:table(as.numeric(d$series))[1]+ m1 <- sum(d[gp1,1] * f[gp1])/sum(f[gp1])+ m2 <- sum(d[-gp1,1] * f[-gp1])/sum(f[-gp1])+ ss1 <- sum(d[gp1,1]^2 * f[gp1]) - (m1 * m1 * sum(f[gp1]))+ ss2 <- sum(d[-gp1,1]^2 * f[-gp1]) - (m2 * m2 * sum(f[-gp1]))+ c(m1 - m2, (ss1 + ss2)/(sum(f) - 2))+ }> grav1 <- gravity[as.numeric(gravity[,2]) >= 7,]> boot(grav1, diff.means, R = 999, stype = "f", strata = grav1[,2])STRATIFIED BOOTSTRAPCall:boot(data = grav1, statistic = diff.means, R = 999, stype = "f",strata = grav1[, 2])Bootstrap Statistics :original bias std. errort1* -2.8462 0.002541 1.5467t2* 16.8462 -1.457663 6.7591>> # In this example we show the use of boot in a prediction from> # regression based on the nuclear data. This example is taken> # from Example 6.8 of Davison and Hinkley (1997). Notice also> # that two extra arguments to 'statistic' are passed through boot.> nuke <- nuclear[, c(1, 2, 5, 7, 8, 10, 11)]> nuke.lm <- glm(log(cost) ~ date+log(cap)+ne+ct+log(cum.n)+pt, data = nuke)> nuke.diag <- glm.diag(nuke.lm)> nuke.res <- nuke.diag$res * nuke.diag$sd> nuke.res <- nuke.res - mean(nuke.res)>> # We set up a new data frame with the data, the standardized> # residuals and the fitted values for use in the bootstrap.> nuke.data <- data.frame(nuke, resid = nuke.res, fit = fitted(nuke.lm))>> # Now we want a prediction of plant number 32 but at date 73.00> new.data <- data.frame(cost = 1, date = 73.00, cap = 886, ne = 0,+ ct = 0, cum.n = 11, pt = 1)> new.fit <- predict(nuke.lm, new.data)>> nuke.fun <- function(dat, inds, i.pred, fit.pred, x.pred)+ {+ lm.b <- glm(fit+resid[inds] ~ date+log(cap)+ne+ct+log(cum.n)+pt,+ data = dat)+ pred.b <- predict(lm.b, x.pred)+ c(coef(lm.b), pred.b - (fit.pred + dat$resid[i.pred]))+ }>> nuke.boot <- boot(nuke.data, nuke.fun, R = 999, m = 1,+ fit.pred = new.fit, x.pred = new.data)> # The bootstrap prediction squared error would then be found by> mean(nuke.boot$t[, 8]^2)[1] 0.088157> # Basic bootstrap prediction limits would be> new.fit - sort(nuke.boot$t[, 8])[c(975, 25)][1] 6.1603 7.2988>>> # Finally a parametric bootstrap. For this example we shall look> # at the air-conditioning data. In this example our aim is to test> # the hypothesis that the true value of the index is 1 (i.e. that> # the data come from an exponential distribution) against the> # alternative that the data come from a gamma distribution with> # index not equal to 1.> air.fun <- function(data) {+ ybar <- mean(data$hours)+ para <- c(log(ybar), mean(log(data$hours)))+ ll <- function(k) {+ if (k <= 0) 1e200 else lgamma(k)-k*(log(k)-1-para[1]+para[2])+ }+ khat <- nlm(ll, ybar^2/var(data$hours))$estimate+ c(ybar, khat)+ }>> air.rg <- function(data, mle) {+ # Function to generate random exponential variates.+ # mle will contain the mean of the original data+ out <- data+ out$hours <- rexp(nrow(out), 1/mle)+ out+ }>> air.boot <- boot(aircondit, air.fun, R = 999, sim = "parametric",+ ran.gen = air.rg, mle = mean(aircondit$hours))>> # The bootstrap p-value can then be approximated by> sum(abs(air.boot$t[,2]-1) > abs(air.boot$t0[2]-1))/(1+air.boot$R)[1] 0.461> ## Don't show:> options(op)> ## End(Don't show)>>>> cleanEx()> nameEx("boot.array")> ### * boot.array>> flush(stderr()); flush(stdout())>> ### Name: boot.array> ### Title: Bootstrap Resampling Arrays> ### Aliases: boot.array> ### Keywords: nonparametric>> ### ** Examples>> # A frequency array for a nonparametric bootstrap> city.boot <- boot(city, corr, R = 40, stype = "w")> boot.array(city.boot)[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10][1,] 1 1 3 1 1 0 1 0 1 1[2,] 0 0 1 3 1 2 1 1 1 0[3,] 0 0 0 1 3 2 0 1 2 1[4,] 0 2 2 2 0 1 0 1 0 2[5,] 1 1 1 0 1 2 0 2 2 0[6,] 0 1 2 0 2 1 0 1 2 1[7,] 3 0 0 1 1 2 0 1 1 1[8,] 0 4 1 1 1 0 1 2 0 0[9,] 0 1 4 0 1 0 2 2 0 0[10,] 1 1 1 1 1 1 1 2 0 1[11,] 0 0 3 0 2 1 1 1 0 2[12,] 2 3 1 0 0 1 0 0 2 1[13,] 1 0 0 1 2 0 2 1 1 2[14,] 0 0 2 3 0 0 2 0 2 1[15,] 2 0 1 1 1 1 0 2 1 1[16,] 1 1 0 2 2 1 1 1 1 0[17,] 1 0 1 2 2 1 2 1 0 0[18,] 0 2 0 0 2 2 0 1 1 2[19,] 1 1 0 2 0 1 2 0 1 2[20,] 0 0 0 1 2 1 1 1 0 4[21,] 0 0 2 0 1 1 3 0 0 3[22,] 1 3 2 2 0 1 1 0 0 0[23,] 0 0 3 1 2 1 2 0 1 0[24,] 0 1 1 2 1 2 0 1 0 2[25,] 0 0 2 1 0 0 3 2 1 1[26,] 0 2 4 1 1 1 0 0 0 1[27,] 2 3 1 0 2 0 0 2 0 0[28,] 1 1 0 1 1 0 0 4 2 0[29,] 2 2 0 1 1 0 1 0 1 2[30,] 0 1 0 1 1 2 0 1 3 1[31,] 1 2 0 2 2 0 1 1 0 1[32,] 1 4 0 1 0 2 0 1 1 0[33,] 0 1 0 5 1 0 0 1 1 1[34,] 0 2 0 1 3 1 1 1 0 1[35,] 0 2 3 1 1 1 0 0 1 1[36,] 1 2 1 1 1 1 2 0 1 0[37,] 1 1 1 0 0 1 2 2 2 0[38,] 2 3 0 3 0 1 0 0 1 0[39,] 0 0 1 2 3 0 0 2 1 1[40,] 0 1 1 0 3 0 2 2 0 1>> perm.cor <- function(d,i) cor(d$x,d$u[i])> city.perm <- boot(city, perm.cor, R = 40, sim = "permutation")> boot.array(city.perm, indices = TRUE)[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10][1,] 7 2 8 10 6 4 9 3 1 5[2,] 10 4 6 8 5 3 2 1 9 7[3,] 5 10 6 2 3 8 1 9 7 4[4,] 10 1 4 8 9 5 3 6 7 2[5,] 1 3 5 8 2 9 6 7 10 4[6,] 10 4 3 2 1 9 8 5 7 6[7,] 2 5 1 4 10 9 3 8 7 6[8,] 5 9 6 3 1 2 4 10 8 7[9,] 8 2 7 10 4 3 1 6 9 5[10,] 1 9 3 2 8 6 7 4 5 10[11,] 6 7 10 5 3 9 2 8 4 1[12,] 9 5 1 6 10 8 3 2 7 4[13,] 9 7 3 1 8 5 4 6 2 10[14,] 7 1 3 2 9 6 10 4 5 8[15,] 4 9 6 3 2 1 10 8 7 5[16,] 9 6 1 7 5 2 8 4 10 3[17,] 7 2 3 6 10 4 1 9 5 8[18,] 2 4 3 1 5 8 6 10 9 7[19,] 4 9 6 1 10 3 7 2 5 8[20,] 6 3 7 5 8 4 2 10 1 9[21,] 9 10 2 6 7 5 8 4 3 1[22,] 10 4 5 9 8 3 1 2 6 7[23,] 9 2 5 1 10 4 6 8 7 3[24,] 9 4 3 6 8 2 1 10 5 7[25,] 10 3 8 2 5 7 1 9 6 4[26,] 7 8 3 9 4 1 5 2 10 6[27,] 4 8 1 5 3 6 10 9 7 2[28,] 8 5 7 4 10 3 2 1 9 6[29,] 8 10 2 4 7 3 9 6 1 5[30,] 1 5 7 9 3 6 4 10 2 8[31,] 10 9 7 5 4 1 2 6 3 8[32,] 8 5 1 6 3 7 10 4 9 2[33,] 1 6 8 5 2 7 9 4 10 3[34,] 8 5 7 1 9 2 6 3 10 4[35,] 4 5 9 2 7 6 8 3 1 10[36,] 5 9 10 1 3 7 2 8 4 6[37,] 8 2 7 9 10 1 3 4 6 5[38,] 5 8 1 7 4 3 9 10 6 2[39,] 10 7 6 4 8 1 3 5 9 2[40,] 2 1 8 3 4 6 9 10 7 5>>>> cleanEx()> nameEx("boot.ci")> ### * boot.ci>> flush(stderr()); flush(stdout())>> ### Name: boot.ci> ### Title: Nonparametric Bootstrap Confidence Intervals> ### Aliases: boot.ci> ### Keywords: nonparametric htest>> ### ** Examples>> # confidence intervals for the city data> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)> city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")> boot.ci(city.boot, conf = c(0.90, 0.95),+ type = c("norm", "basic", "perc", "bca"))BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 999 bootstrap replicatesCALL :boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm","basic", "perc", "bca"))Intervals :Level Normal Basic90% ( 1.124, 1.837 ) ( 1.059, 1.740 )95% ( 1.056, 1.905 ) ( 0.932, 1.799 )Level Percentile BCa90% ( 1.301, 1.982 ) ( 1.301, 1.984 )95% ( 1.242, 2.109 ) ( 1.243, 2.110 )Calculations and Intervals on Original Scale>> # studentized confidence interval for the two sample> # difference of means problem using the final two series> # of the gravity data.> diff.means <- function(d, f)+ { n <- nrow(d)+ gp1 <- 1:table(as.numeric(d$series))[1]+ m1 <- sum(d[gp1,1] * f[gp1])/sum(f[gp1])+ m2 <- sum(d[-gp1,1] * f[-gp1])/sum(f[-gp1])+ ss1 <- sum(d[gp1,1]^2 * f[gp1]) - (m1 * m1 * sum(f[gp1]))+ ss2 <- sum(d[-gp1,1]^2 * f[-gp1]) - (m2 * m2 * sum(f[-gp1]))+ c(m1 - m2, (ss1 + ss2)/(sum(f) - 2))+ }> grav1 <- gravity[as.numeric(gravity[,2]) >= 7, ]> grav1.boot <- boot(grav1, diff.means, R = 999, stype = "f",+ strata = grav1[ ,2])> boot.ci(grav1.boot, type = c("stud", "norm"))BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 999 bootstrap replicatesCALL :boot.ci(boot.out = grav1.boot, type = c("stud", "norm"))Intervals :Level Normal Studentized95% (-5.880, 0.183 ) (-7.059, -0.101 )Calculations and Intervals on Original Scale>> # Nonparametric confidence intervals for mean failure time> # of the air-conditioning data as in Example 5.4 of Davison> # and Hinkley (1997)> mean.fun <- function(d, i)+ { m <- mean(d$hours[i])+ n <- length(i)+ v <- (n-1)*var(d$hours[i])/n^2+ c(m, v)+ }> air.boot <- boot(aircondit, mean.fun, R = 999)> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"))BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 999 bootstrap replicatesCALL :boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc","stud"))Intervals :Level Normal Basic95% ( 35.5, 181.9 ) ( 26.0, 170.6 )Level Studentized Percentile95% ( 47.9, 294.5 ) ( 45.6, 190.2 )Calculations and Intervals on Original Scale>> # Now using the log transformation> # There are two ways of doing this and they both give the> # same intervals.>> # Method 1> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"),+ h = log, hdot = function(x) 1/x)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 999 bootstrap replicatesCALL :boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc","stud"), h = log, hdot = function(x) 1/x)Intervals :Level Normal Basic95% ( 4.035, 5.469 ) ( 4.118, 5.546 )Level Studentized Percentile95% ( 3.959, 5.808 ) ( 3.820, 5.248 )Calculations and Intervals on Transformed Scale>> # Method 2> vt0 <- air.boot$t0[2]/air.boot$t0[1]^2> vt <- air.boot$t[, 2]/air.boot$t[ ,1]^2> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"),+ t0 = log(air.boot$t0[1]), t = log(air.boot$t[,1]),+ var.t0 = vt0, var.t = vt)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 999 bootstrap replicatesCALL :boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc","stud"), var.t0 = vt0, var.t = vt, t0 = log(air.boot$t0[1]),t = log(air.boot$t[, 1]))Intervals :Level Normal Basic95% ( 4.069, 5.435 ) ( 4.118, 5.546 )Level Studentized Percentile95% ( 3.959, 5.808 ) ( 3.820, 5.248 )Calculations and Intervals on Original Scale>>>> cleanEx()> nameEx("censboot")> ### * censboot>> flush(stderr()); flush(stdout())>> ### Name: censboot> ### Title: Bootstrap for Censored Data> ### Aliases: censboot cens.return> ### Keywords: survival>> ### ** Examples>> library(survival)Attaching package: 'survival'The following object is masked from 'package:boot':aml> # Example 3.9 of Davison and Hinkley (1997) does a bootstrap on some> # remission times for patients with a type of leukaemia. The patients> # were divided into those who received maintenance chemotherapy and> # those who did not. Here we are interested in the median remission> # time for the two groups.> data(aml, package = "boot") # not the version in survival.> aml.fun <- function(data) {+ surv <- survfit(Surv(time, cens) ~ group, data = data)+ out <- NULL+ st <- 1+ for (s in 1:length(surv$strata)) {+ inds <- st:(st + surv$strata[s]-1)+ md <- min(surv$time[inds[1-surv$surv[inds] >= 0.5]])+ st <- st + surv$strata[s]+ out <- c(out, md)+ }+ out+ }> aml.case <- censboot(aml, aml.fun, R = 499, strata = aml$group)Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning Inf>> # Now we will look at the same statistic using the conditional> # bootstrap and the weird bootstrap. For the conditional bootstrap> # the survival distribution is stratified but the censoring> # distribution is not.>> aml.s1 <- survfit(Surv(time, cens) ~ group, data = aml)> aml.s2 <- survfit(Surv(time-0.001*cens, 1-cens) ~ 1, data = aml)> aml.cond <- censboot(aml, aml.fun, R = 499, strata = aml$group,+ F.surv = aml.s1, G.surv = aml.s2, sim = "cond")Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning Inf>>> # For the weird bootstrap we must redefine our function slightly since> # the data will not contain the group number.> aml.fun1 <- function(data, str) {+ surv <- survfit(Surv(data[, 1], data[, 2]) ~ str)+ out <- NULL+ st <- 1+ for (s in 1:length(surv$strata)) {+ inds <- st:(st + surv$strata[s] - 1)+ md <- min(surv$time[inds[1-surv$surv[inds] >= 0.5]])+ st <- st + surv$strata[s]+ out <- c(out, md)+ }+ out+ }> aml.wei <- censboot(cbind(aml$time, aml$cens), aml.fun1, R = 499,+ strata = aml$group, F.surv = aml.s1, sim = "weird")Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning InfWarning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :no non-missing arguments to min; returning Inf>> # Now for an example where a cox regression model has been fitted> # the data we will look at the melanoma data of Example 7.6 from> # Davison and Hinkley (1997). The fitted model assumes that there> # is a different survival distribution for the ulcerated and> # non-ulcerated groups but that the thickness of the tumour has a> # common effect. We will also assume that the censoring distribution> # is different in different age groups. The statistic of interest> # is the linear predictor. This is returned as the values at a> # number of equally spaced points in the range of interest.> data(melanoma, package = "boot")> library(splines)# for ns> mel.cox <- coxph(Surv(time, status == 1) ~ ns(thickness, df=4) + strata(ulcer),+ data = melanoma)> mel.surv <- survfit(mel.cox)> agec <- cut(melanoma$age, c(0, 39, 49, 59, 69, 100))> mel.cens <- survfit(Surv(time - 0.001*(status == 1), status != 1) ~+ strata(agec), data = melanoma)> mel.fun <- function(d) {+ t1 <- ns(d$thickness, df=4)+ cox <- coxph(Surv(d$time, d$status == 1) ~ t1+strata(d$ulcer))+ ind <- !duplicated(d$thickness)+ u <- d$thickness[!ind]+ eta <- cox$linear.predictors[!ind]+ sp <- smooth.spline(u, eta, df=20)+ th <- seq(from = 0.25, to = 10, by = 0.25)+ predict(sp, th)$y+ }> mel.str <- cbind(melanoma$ulcer, agec)>> # this is slow!> mel.mod <- censboot(melanoma, mel.fun, R = 499, F.surv = mel.surv,+ G.surv = mel.cens, cox = mel.cox, strata = mel.str, sim = "model")> # To plot the original predictor and a 95% pointwise envelope for it> mel.env <- envelope(mel.mod)$point> th <- seq(0.25, 10, by = 0.25)> plot(th, mel.env[1, ], ylim = c(-2, 2),+ xlab = "thickness (mm)", ylab = "linear predictor", type = "n")> lines(th, mel.mod$t0, lty = 1)> matlines(th, t(mel.env), lty = 2)>>>> cleanEx()detaching 'package:splines', 'package:survival'> nameEx("control")> ### * control>> flush(stderr()); flush(stdout())>> ### Name: control> ### Title: Control Variate Calculations> ### Aliases: control> ### Keywords: nonparametric>> ### ** Examples>> # Use of control variates for the variance of the air-conditioning data> mean.fun <- function(d, i)+ { m <- mean(d$hours[i])+ n <- nrow(d)+ v <- (n-1)*var(d$hours[i])/n^2+ c(m, v)+ }> air.boot <- boot(aircondit, mean.fun, R = 999)> control(air.boot, index = 2, bias.adj = TRUE)[1] -6.298101> air.cont <- control(air.boot, index = 2)> # Now let us try the variance on the log scale.> air.cont1 <- control(air.boot, t0 = log(air.boot$t0[2]),+ t = log(air.boot$t[, 2]))>>>> cleanEx()> nameEx("cv.glm")> ### * cv.glm>> flush(stderr()); flush(stdout())>> ### Name: cv.glm> ### Title: Cross-validation for Generalized Linear Models> ### Aliases: cv.glm> ### Keywords: regression>> ### ** Examples>> # leave-one-out and 6-fold cross-validation prediction error for> # the mammals data set.> data(mammals, package="MASS")> mammals.glm <- glm(log(brain) ~ log(body), data = mammals)> (cv.err <- cv.glm(mammals, mammals.glm)$delta)[1] 0.4918650 0.4916571> (cv.err.6 <- cv.glm(mammals, mammals.glm, K = 6)$delta)[1] 0.4851491 0.4834641>> # As this is a linear model we could calculate the leave-one-out> # cross-validation estimate without any extra model-fitting.> muhat <- fitted(mammals.glm)> mammals.diag <- glm.diag(mammals.glm)> (cv.err <- mean((mammals.glm$y - muhat)^2/(1 - mammals.diag$h)^2))[1] 0.491865>>> # leave-one-out and 11-fold cross-validation prediction error for> # the nodal data set. Since the response is a binary variable an> # appropriate cost function is> cost <- function(r, pi = 0) mean(abs(r-pi) > 0.5)>> nodal.glm <- glm(r ~ stage+xray+acid, binomial, data = nodal)> (cv.err <- cv.glm(nodal, nodal.glm, cost, K = nrow(nodal))$delta)[1] 0.1886792 0.1886792> (cv.11.err <- cv.glm(nodal, nodal.glm, cost, K = 11)$delta)[1] 0.2264151 0.2217871>>>> cleanEx()> nameEx("empinf")> ### * empinf>> flush(stderr()); flush(stdout())>> ### Name: empinf> ### Title: Empirical Influence Values> ### Aliases: empinf> ### Keywords: nonparametric math>> ### ** Examples>> # The empirical influence values for the ratio of means in> # the city data.> ratio <- function(d, w) sum(d$x *w)/sum(d$u*w)> empinf(data = city, statistic = ratio)[1] -1.04367815 -0.58417763 -0.37092459 -0.18958996 0.03164142 0.10544878[7] 0.09236345 0.20365074 1.02178280 0.73381132> city.boot <- boot(city, ratio, 499, stype="w")> empinf(boot.out = city.boot, type = "reg")1 1 1 1 1 1-1.13619987 -0.69728210 -0.45301061 -0.27615882 0.02108999 0.148963361 1 1 10.09746429 0.20622340 1.18798956 0.90092079>> # A statistic that may be of interest in the difference of means> # problem is the t-statistic for testing equality of means. In> # the bootstrap we get replicates of the difference of means and> # the variance of that statistic and then want to use this output> # to get the empirical influence values of the t-statistic.> grav1 <- gravity[as.numeric(gravity[,2]) >= 7,]> grav.fun <- function(dat, w) {+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))+ d <- dat[, 1]+ ns <- tabulate(strata)+ w <- w/tapply(w, strata, sum)[strata]+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names+ mn2 <- tapply(d * d * w, strata, sum)+ s2hat <- sum((mn2 - mns^2)/ns)+ c(mns[2] - mns[1], s2hat)+ }>> grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w",+ strata = grav1[, 2])>> # Since the statistic of interest is a function of the bootstrap> # statistics, we must calculate the bootstrap replicates and pass> # them to empinf using the t argument.> grav.z <- (grav.boot$t[,1]-grav.boot$t0[1])/sqrt(grav.boot$t[,2])> empinf(boot.out = grav.boot, t = grav.z)1 1 1 1 1 1 1-2.9326019 -1.3760327 -2.4400720 -1.2175846 0.2795352 -0.8258764 -0.81562861 1 1 1 1 1 2-0.5573332 -1.1275252 -3.1603140 1.2840693 3.5434781 9.3458860 2.66925892 2 2 2 2 2 24.4496570 3.6948000 0.9929002 -3.0100985 -3.2237464 -2.5493305 -0.65517452 2 2 2 21.9065308 0.4980530 -1.6219628 -1.6980508 -1.4528364>>>> cleanEx()> nameEx("envelope")> ### * envelope>> flush(stderr()); flush(stdout())>> ### Name: envelope> ### Title: Confidence Envelopes for Curves> ### Aliases: envelope> ### Keywords: dplot htest>> ### ** Examples>> # Testing whether the final series of measurements of the gravity data> # may come from a normal distribution. This is done in Examples 4.7> # and 4.8 of Davison and Hinkley (1997).> grav1 <- gravity$g[gravity$series == 8]> grav.z <- (grav1 - mean(grav1))/sqrt(var(grav1))> grav.gen <- function(dat, mle) rnorm(length(dat))> grav.qqboot <- boot(grav.z, sort, R = 999, sim = "parametric",+ ran.gen = grav.gen)> grav.qq <- qqnorm(grav.z, plot.it = FALSE)> grav.qq <- lapply(grav.qq, sort)> plot(grav.qq, ylim = c(-3.5, 3.5), ylab = "Studentized Order Statistics",+ xlab = "Normal Quantiles")> grav.env <- envelope(grav.qqboot, level = 0.9)> lines(grav.qq$x, grav.env$point[1, ], lty = 4)> lines(grav.qq$x, grav.env$point[2, ], lty = 4)> lines(grav.qq$x, grav.env$overall[1, ], lty = 1)> lines(grav.qq$x, grav.env$overall[2, ], lty = 1)>>>> cleanEx()> nameEx("exp.tilt")> ### * exp.tilt>> flush(stderr()); flush(stdout())>> ### Name: exp.tilt> ### Title: Exponential Tilting> ### Aliases: exp.tilt> ### Keywords: nonparametric smooth>> ### ** Examples>> # 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. This can be achieved as follows.> 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 <- as.vector(tapply(d * w, strata, sum)) # drop names+ 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.L <- empinf(data = grav1, statistic = grav.fun, stype = "w",+ strata = grav1[,2], index = 3, orig = grav.z0[1])> grav.tilt <- exp.tilt(grav.L, grav.z0[3], strata = grav1[,2])> boot(grav1, grav.fun, R = 499, stype = "w", weights = grav.tilt$p,+ strata = grav1[,2], orig = grav.z0[1])STRATIFIED WEIGHTED BOOTSTRAPCall:boot(data = grav1, statistic = grav.fun, R = 499, stype = "w",strata = grav1[, 2], weights = grav.tilt$p, orig = grav.z0[1])Bootstrap Statistics :original bias std. error mean(t*)t1* 2.846154 -0.3661063 1.705171 5.702944t2* 2.392353 -0.3538294 1.002889 3.444050t3* 0.000000 -0.5160619 1.314298 1.473456>>>> cleanEx()> nameEx("glm.diag.plots")> ### * glm.diag.plots>> flush(stderr()); flush(stdout())>> ### Name: glm.diag.plots> ### Title: Diagnostics plots for generalized linear models> ### Aliases: glm.diag.plots> ### Keywords: regression dplot hplot>> ### ** Examples>> # In this example we look at the leukaemia data which was looked at in> # Example 7.1 of Davison and Hinkley (1997)> data(leuk, package = "MASS")> leuk.mod <- glm(time ~ ag-1+log10(wbc), family = Gamma(log), data = leuk)> leuk.diag <- glm.diag(leuk.mod)> glm.diag.plots(leuk.mod, leuk.diag)>>>> cleanEx()> nameEx("jack.after.boot")> ### * jack.after.boot>> flush(stderr()); flush(stdout())>> ### Name: jack.after.boot> ### Title: Jackknife-after-Bootstrap Plots> ### Aliases: jack.after.boot> ### Keywords: hplot nonparametric>> ### ** Examples>> # To draw the jackknife-after-bootstrap plot for the head size data as in> # Example 3.24 of Davison and Hinkley (1997)> frets.fun <- function(data, i) {+ pcorr <- function(x) {+ # Function to find the correlations and partial correlations between+ # the four measurements.+ v <- cor(x)+ v.d <- diag(var(x))+ iv <- solve(v)+ iv.d <- sqrt(diag(iv))+ iv <- - diag(1/iv.d) %*% iv %*% diag(1/iv.d)+ q <- NULL+ n <- nrow(v)+ for (i in 1:(n-1))+ q <- rbind( q, c(v[i, 1:i], iv[i,(i+1):n]) )+ q <- rbind( q, v[n, ] )+ diag(q) <- round(diag(q))+ q+ }+ d <- data[i, ]+ v <- pcorr(d)+ c(v[1,], v[2,], v[3,], v[4,])+ }> frets.boot <- boot(log(as.matrix(frets)), frets.fun, R = 999)> # we will concentrate on the partial correlation between head breadth> # for the first son and head length for the second. This is the 7th> # element in the output of frets.fun so we set index = 7> jack.after.boot(frets.boot, useJ = FALSE, stinf = FALSE, index = 7)>>>> cleanEx()> nameEx("k3.linear")> ### * k3.linear>> flush(stderr()); flush(stdout())>> ### Name: k3.linear> ### Title: Linear Skewness Estimate> ### Aliases: k3.linear> ### Keywords: nonparametric>> ### ** Examples>> # To estimate the skewness of the ratio of means for the city data.> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)> k3.linear(empinf(data = city, statistic = ratio))[1] 7.831452e-05>>>> cleanEx()> nameEx("linear.approx")> ### * linear.approx>> flush(stderr()); flush(stdout())>> ### Name: linear.approx> ### Title: Linear Approximation of Bootstrap Replicates> ### Aliases: linear.approx> ### Keywords: nonparametric>> ### ** Examples>> # Using the city data let us look at the linear approximation to the> # ratio statistic and its logarithm. We compare these with the> # corresponding plots for the bigcity data>> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)> city.boot <- boot(city, ratio, R = 499, stype = "w")> bigcity.boot <- boot(bigcity, ratio, R = 499, stype = "w")> op <- par(pty = "s", mfrow = c(2, 2))>> # The first plot is for the city data ratio statistic.> city.lin1 <- linear.approx(city.boot)> lim <- range(c(city.boot$t,city.lin1))> plot(city.boot$t, city.lin1, xlim = lim, ylim = lim,+ main = "Ratio; n=10", xlab = "t*", ylab = "tL*")> abline(0, 1)>> # Now for the log of the ratio statistic for the city data.> city.lin2 <- linear.approx(city.boot,t0 = log(city.boot$t0),+ t = log(city.boot$t))> lim <- range(c(log(city.boot$t),city.lin2))> plot(log(city.boot$t), city.lin2, xlim = lim, ylim = lim,+ main = "Log(Ratio); n=10", xlab = "t*", ylab = "tL*")> abline(0, 1)>> # The ratio statistic for the bigcity data.> bigcity.lin1 <- linear.approx(bigcity.boot)> lim <- range(c(bigcity.boot$t,bigcity.lin1))> plot(bigcity.lin1, bigcity.boot$t, xlim = lim, ylim = lim,+ main = "Ratio; n=49", xlab = "t*", ylab = "tL*")> abline(0, 1)>> # Finally the log of the ratio statistic for the bigcity data.> bigcity.lin2 <- linear.approx(bigcity.boot,t0 = log(bigcity.boot$t0),+ t = log(bigcity.boot$t))> lim <- range(c(log(bigcity.boot$t),bigcity.lin2))> plot(bigcity.lin2, log(bigcity.boot$t), xlim = lim, ylim = lim,+ main = "Log(Ratio); n=49", xlab = "t*", ylab = "tL*")> abline(0, 1)>> par(op)>>>> graphics::par(get("par.postscript", pos = 'CheckExEnv'))> cleanEx()> nameEx("lines.saddle.distn")> ### * lines.saddle.distn>> flush(stderr()); flush(stdout())>> ### Name: lines.saddle.distn> ### Title: Add a Saddlepoint Approximation to a Plot> ### Aliases: lines.saddle.distn> ### Keywords: aplot smooth nonparametric>> ### ** Examples>> # In this example we show how a plot such as that in Figure 9.9 of> # Davison and Hinkley (1997) may be produced. Note the large number of> # bootstrap replicates required in this example.> expdata <- rexp(12)> vfun <- function(d, i) {+ n <- length(d)+ (n-1)/n*var(d[i])+ }> exp.boot <- boot(expdata,vfun, R = 9999)> exp.L <- (expdata - mean(expdata))^2 - exp.boot$t0> exp.tL <- linear.approx(exp.boot, L = exp.L)> hist(exp.tL, nclass = 50, probability = TRUE)> exp.t0 <- c(0, sqrt(var(exp.boot$t)))> exp.sp <- saddle.distn(A = exp.L/12,wdist = "m", t0 = exp.t0)>> # The saddlepoint approximation in this case is to the density of> # t-t0 and so t0 must be added for the plot.> lines(exp.sp, h = function(u, t0) u+t0, J = function(u, t0) 1,+ t0 = exp.boot$t0)>>>> cleanEx()> nameEx("norm.ci")> ### * norm.ci>> flush(stderr()); flush(stdout())>> ### Name: norm.ci> ### Title: Normal Approximation Confidence Intervals> ### Aliases: norm.ci> ### Keywords: htest>> ### ** Examples>> # In Example 5.1 of Davison and Hinkley (1997), normal approximation> # confidence intervals are found for the air-conditioning data.> air.mean <- mean(aircondit$hours)> air.n <- nrow(aircondit)> air.v <- air.mean^2/air.n> norm.ci(t0 = air.mean, var.t0 = air.v)conf[1,] 0.95 46.93055 169.2361> exp(norm.ci(t0 = log(air.mean), var.t0 = 1/air.n)[2:3])[1] 61.38157 190.31782>> # Now a more complicated example - the ratio estimate for the city data.> ratio <- function(d, w)+ sum(d$x * w)/sum(d$u *w)> city.v <- var.linear(empinf(data = city, statistic = ratio))> norm.ci(t0 = ratio(city,rep(0.1,10)), var.t0 = city.v)conf[1,] 0.95 1.167046 1.873579>>>> cleanEx()> nameEx("plot.boot")> ### * plot.boot>> flush(stderr()); flush(stdout())>> ### Name: plot.boot> ### Title: Plots of the Output of a Bootstrap Simulation> ### Aliases: plot.boot> ### Keywords: hplot nonparametric>> ### ** Examples>> # We fit an exponential model to the air-conditioning data and use> # that for a parametric bootstrap. Then we look at plots of the> # resampled means.> air.rg <- function(data, mle) rexp(length(data), 1/mle)>> air.boot <- boot(aircondit$hours, mean, R = 999, sim = "parametric",+ ran.gen = air.rg, mle = mean(aircondit$hours))> plot(air.boot)>> # In the difference of means example for the last two series of the> # gravity data> grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]> grav.fun <- function(dat, w) {+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))+ d <- dat[, 1]+ ns <- tabulate(strata)+ w <- w/tapply(w, strata, sum)[strata]+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names+ mn2 <- tapply(d * d * w, strata, sum)+ s2hat <- sum((mn2 - mns^2)/ns)+ c(mns[2] - mns[1], s2hat)+ }>> grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w", strata = grav1[, 2])> plot(grav.boot)> # now suppose we want to look at the studentized differences.> grav.z <- (grav.boot$t[, 1]-grav.boot$t0[1])/sqrt(grav.boot$t[, 2])> plot(grav.boot, t = grav.z, t0 = 0)>> # In this example we look at the one of the partial correlations for the> # head dimensions in the dataset frets.> frets.fun <- function(data, i) {+ pcorr <- function(x) {+ # Function to find the correlations and partial correlations between+ # the four measurements.+ v <- cor(x)+ v.d <- diag(var(x))+ iv <- solve(v)+ iv.d <- sqrt(diag(iv))+ iv <- - diag(1/iv.d) %*% iv %*% diag(1/iv.d)+ q <- NULL+ n <- nrow(v)+ for (i in 1:(n-1))+ q <- rbind( q, c(v[i, 1:i], iv[i,(i+1):n]) )+ q <- rbind( q, v[n, ] )+ diag(q) <- round(diag(q))+ q+ }+ d <- data[i, ]+ v <- pcorr(d)+ c(v[1,], v[2,], v[3,], v[4,])+ }> frets.boot <- boot(log(as.matrix(frets)), frets.fun, R = 999)> plot(frets.boot, index = 7, jack = TRUE, stinf = FALSE, useJ = FALSE)>>>> cleanEx()> nameEx("saddle")> ### * saddle>> flush(stderr()); flush(stdout())>> ### Name: saddle> ### Title: Saddlepoint Approximations for Bootstrap Statistics> ### Aliases: saddle> ### Keywords: smooth nonparametric>> ### ** Examples>> # To evaluate the bootstrap distribution of the mean failure time of> # air-conditioning equipment at 80 hours> saddle(A = aircondit$hours/12, u = 80)$spapdf cdf0.01005866 0.24446677$zeta.hat[1] -0.02580078>> # Alternatively this can be done using a conditional poisson> saddle(A = cbind(aircondit$hours/12,1), u = c(80, 12),+ wdist = "p", type = "cond")Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091$spapdf cdf0.01005943 0.24438736$zeta.hatA1 A2-0.02580805 0.89261577$zeta2.hat[1] 0.6931472>> # To use the Lugananni-Rice approximation to this> saddle(A = cbind(aircondit$hours/12,1), u = c(80, 12),+ wdist = "p", type = "cond",+ LR = TRUE)Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091$spapdf cdf0.01005943 0.24447362$zeta.hatA1 A2-0.02580805 0.89261577$zeta2.hat[1] 0.6931472>> # Example 9.16 of Davison and Hinkley (1997) calculates saddlepoint> # approximations to the distribution of the ratio statistic for the> # city data. Since the statistic is not in itself a linear combination> # of random Variables, its distribution cannot be found directly.> # Instead the statistic is expressed as the solution to a linear> # estimating equation and hence its distribution can be found. We> # get the saddlepoint approximation to the pdf and cdf evaluated at> # t = 1.25 as follows.> jacobian <- function(dat,t,zeta)+ {+ p <- exp(zeta*(dat$x-t*dat$u))+ abs(sum(dat$u*p)/sum(p))+ }> city.sp1 <- saddle(A = city$x-1.25*city$u, u = 0)> city.sp1$spa[1] <- jacobian(city, 1.25, city.sp1$zeta.hat) * city.sp1$spa[1]> city.sp1$spapdf cdf0.05565040 0.02436306$zeta.hat[1] -0.02435547>>>> cleanEx()> nameEx("saddle.distn")> ### * saddle.distn>> flush(stderr()); flush(stdout())>> ### Name: saddle.distn> ### Title: Saddlepoint Distribution Approximations for Bootstrap Statistics> ### Aliases: saddle.distn> ### Keywords: nonparametric smooth dplot>> ### ** Examples>> # The bootstrap distribution of the mean of the air-conditioning> # failure data: fails to find value on R (and probably on S too)> air.t0 <- c(mean(aircondit$hours), sqrt(var(aircondit$hours)/12))> ## Not run: saddle.distn(A = aircondit$hours/12, t0 = air.t0)>> # alternatively using the conditional poisson> saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p",+ type = "cond", t0 = air.t0)Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.344718Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.655282Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.344718Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.655282Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.832240Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.167760Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.832240Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.167760Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.444538Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.555462Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.444538Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.555462Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.494449Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.505551Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.494449Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.505551Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.594269Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.405731Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.594269Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.405731Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.544359Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.455641Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.544359Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.455641Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.519404Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.480596Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.519404Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.480596Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.561394Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.438606Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.561394Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.438606Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.290549Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.709451Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.290549Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.709451Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.019703Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.980297Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.019703Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.980297Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.748857Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.251143Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.748857Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.251143Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.478011Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.521989Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.478011Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.521989Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.207165Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.792835Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.207165Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.792835Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.936320Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.063680Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.936320Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.063680Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.665474Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.334526Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.665474Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.334526Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.785225Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.214775Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.785225Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.214775Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.175822Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.824178Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.175822Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.824178Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.566419Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.433581Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.566419Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.433581Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.957016Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.042984Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.957016Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.042984Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.347613Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.652387Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.347613Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.652387Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.738210Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.261790Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.738210Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.261790Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.128807Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.871193Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.128807Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.871193Saddlepoint Distribution ApproximationsCall :saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p",type = "cond", t0 = air.t0)Quantiles of the Distribution0.1% 27.40.5% 35.41.0% 39.72.5% 46.75.0% 53.510.0% 62.520.0% 75.350.0% 104.580.0% 139.090.0% 158.895.0% 175.997.5% 191.299.0% 209.699.5% 222.499.9% 249.5Smoothing spline used 20 points in the range 9.8 to 304.7.>> # Distribution of the ratio of a sample of size 10 from the bigcity> # data, taken from Example 9.16 of Davison and Hinkley (1997).> ratio <- function(d, w) sum(d$x *w)/sum(d$u * w)> city.v <- var.linear(empinf(data = city, statistic = ratio))> bigcity.t0 <- c(mean(bigcity$x)/mean(bigcity$u), sqrt(city.v))> Afn <- function(t, data) cbind(data$x - t*data$u, 1)> ufn <- function(t, data) c(0,10)> saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond",+ t0 = bigcity.t0, data = bigcity)Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Warning in eval(family$initialize) :non-integer counts in a binomial glm!Saddlepoint Distribution ApproximationsCall :saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond", t0 = bigcity.t0,data = bigcity)Quantiles of the Distribution0.1% 1.0700.5% 1.0921.0% 1.1042.5% 1.1225.0% 1.13910.0% 1.15820.0% 1.18450.0% 1.23780.0% 1.30490.0% 1.34895.0% 1.39297.5% 1.43699.0% 1.49499.5% 1.53799.9% 1.636Smoothing spline used 20 points in the range 1.014 to 1.96.>> # From Example 9.16 of Davison and Hinkley (1997) again, we find the> # conditional distribution of the ratio given the sum of city$u.> Afn <- function(t, data) cbind(data$x-t*data$u, data$u, 1)> ufn <- function(t, data) c(0, sum(data$u), 10)> city.t0 <- c(mean(city$x)/mean(city$u), sqrt(city.v))> saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0,+ data = city)Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.866400Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.511350Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.622251Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.866400Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.511350Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.622251Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.210844Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.107208Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.681949Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.210844Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.107208Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.681949Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.038622Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.809279Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.152100Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.038622Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.809279Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.152100Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.452511Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.160314Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.387175Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.452511Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.160314Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.387175Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.159455Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.835832Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.004713Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.159455Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.835832Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.004713Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.205022Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.133273Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.661705Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.205022Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.133273Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.661705Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.722845Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.033105Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.244049Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.722845Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.033105Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.244049Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.787431Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.388294Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.824275Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.787431Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.388294Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.824275Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.415407Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.940756Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.643837Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.415407Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.940756Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.643837Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.043383Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.493218Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.463399Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.043383Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.493218Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.463399Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.671359Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.045680Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.282961Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.671359Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.045680Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.282961Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.299336Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.598142Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.102523Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.299336Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.598142Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.102523Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.433935Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.409277Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.156788Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.433935Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.409277Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.156788Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.360158Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.271008Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.368834Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.360158Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.271008Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.368834Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.319252Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.639889Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.040859Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.319252Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.639889Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.040859Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.278346Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.008770Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.712884Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.278346Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.008770Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.712884Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.237440Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.377650Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.384909Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.237440Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.377650Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.384909Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.196534Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.746531Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.056935Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.196534Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.746531Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.056935Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.734728Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.299661Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.965612Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.734728Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.299661Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.965612Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.228786Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.666383Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.104831Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.228786Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.666383Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.104831Saddlepoint Distribution ApproximationsCall :saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0,data = city)Quantiles of the Distribution0.1% 1.2160.5% 1.2361.0% 1.2482.5% 1.2725.0% 1.30110.0% 1.34020.0% 1.39350.0% 1.50280.0% 1.61890.0% 1.68095.0% 1.73297.5% 1.77799.0% 1.83099.5% 1.86699.9% 1.938Smoothing spline used 20 points in the range 1.182 to 2.061.>>>> cleanEx()> nameEx("simplex")> ### * simplex>> flush(stderr()); flush(stdout())>> ### Name: simplex> ### Title: Simplex Method for Linear Programming Problems> ### Aliases: simplex> ### Keywords: optimize>> ### ** Examples>> # This example is taken from Exercise 7.5 of Gill, Murray and Wright (1991).> enj <- c(200, 6000, 3000, -200)> fat <- c(800, 6000, 1000, 400)> vitx <- c(50, 3, 150, 100)> vity <- c(10, 10, 75, 100)> vitz <- c(150, 35, 75, 5)> simplex(a = enj, A1 = fat, b1 = 13800, A2 = rbind(vitx, vity, vitz),+ b2 = c(600, 300, 550), maxi = TRUE)Linear Programming ResultsCall : simplex(a = enj, A1 = fat, b1 = 13800, A2 = rbind(vitx, vity,vitz), b2 = c(600, 300, 550), maxi = TRUE)Maximization Problem with Objective Function Coefficientsx1 x2 x3 x4200 6000 3000 -200Optimal solution has the following valuesx1 x2 x3 x40.0 0.0 13.8 0.0The optimal value of the objective function is 41400.>>>> cleanEx()> nameEx("smooth.f")> ### * smooth.f>> flush(stderr()); flush(stdout())>> ### Name: smooth.f> ### Title: Smooth Distributions on Data Points> ### Aliases: smooth.f> ### Keywords: smooth nonparametric>> ### ** Examples>> # 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.> 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 <- as.vector(tapply(d * w, strata, sum)) # drop names+ 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])[1] 0.01402806> imp.prob(grav.boot2, t0 = -grav.z0[3], t = -grav.boot2$t[, 3])$t0[1] -1.840118$raw[1] 0.02163715$rat[1] 0.02099078$reg[1] 0.02174393>>> # 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[1] -0.009155757> # can give very strange results (negative probabilities).>>>> cleanEx()> nameEx("tilt.boot")> ### * tilt.boot>> flush(stderr()); flush(stdout())>> ### Name: tilt.boot> ### Title: Non-parametric Tilted Bootstrap> ### Aliases: tilt.boot> ### Keywords: nonparametric>> ### ** Examples>> # Note that these examples can take a while to run.>> # Example 9.9 of Davison and Hinkley (1997).> 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 <- as.vector(tapply(d * w, strata, sum)) # drop names+ 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)> tilt.boot(grav1, grav.fun, R = c(249, 375, 375), stype = "w",+ strata = grav1[,2], tilt = TRUE, index = 3, orig = grav.z0[1])TILTED BOOTSTRAPExponential tilting usedFirst 249 replicates untilted,Next 375 replicates tilted to -2.821,Next 375 replicates tilted to 1.636.Call:tilt.boot(data = grav1, statistic = grav.fun, R = c(249, 375,375), stype = "w", strata = grav1[, 2], tilt = TRUE, index = 3,orig = grav.z0[1])Bootstrap Statistics :original bias std. errort1* 2.846154 -0.4487564 2.500644t2* 2.392353 -0.3221155 1.187574t3* 0.000000 -0.8862944 2.208945>>> # Example 9.10 of Davison and Hinkley (1997) requires a balanced> # importance resampling bootstrap to be run. In this example we> # show how this might be run.> acme.fun <- function(data, i, bhat) {+ d <- data[i,]+ n <- nrow(d)+ d.lm <- glm(d$acme~d$market)+ beta.b <- coef(d.lm)[2]+ d.diag <- boot::glm.diag(d.lm)+ SSx <- (n-1)*var(d$market)+ tmp <- (d$market-mean(d$market))*d.diag$res*d.diag$sd+ sr <- sqrt(sum(tmp^2))/SSx+ c(beta.b, sr, (beta.b-bhat)/sr)+ }> acme.b <- acme.fun(acme, 1:nrow(acme), 0)> acme.boot1 <- tilt.boot(acme, acme.fun, R = c(499, 250, 250),+ stype = "i", sim = "balanced", alpha = c(0.05, 0.95),+ tilt = TRUE, index = 3, bhat = acme.b[1])>>>> cleanEx()> nameEx("tsboot")> ### * tsboot>> flush(stderr()); flush(stdout())>> ### Name: tsboot> ### Title: Bootstrapping of Time Series> ### Aliases: tsboot ts.return> ### Keywords: nonparametric ts>> ### ** Examples>> lynx.fun <- function(tsb) {+ ar.fit <- ar(tsb, order.max = 25)+ c(ar.fit$order, mean(tsb), tsb)+ }>> # the stationary bootstrap with mean block length 20> lynx.1 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "geom")>> # the fixed block bootstrap with length 20> lynx.2 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "fixed")>> # Now for model based resampling we need the original model> # Note that for all of the bootstraps which use the residuals as their> # data, we set orig.t to FALSE since the function applied to the residual> # time series will be meaningless.> lynx.ar <- ar(log(lynx))> lynx.model <- list(order = c(lynx.ar$order, 0, 0), ar = lynx.ar$ar)> lynx.res <- lynx.ar$resid[!is.na(lynx.ar$resid)]> lynx.res <- lynx.res - mean(lynx.res)>> lynx.sim <- function(res,n.sim, ran.args) {+ # random generation of replicate series using arima.sim+ rg1 <- function(n, res) sample(res, n, replace = TRUE)+ ts.orig <- ran.args$ts+ ts.mod <- ran.args$model+ mean(ts.orig)+ts(arima.sim(model = ts.mod, n = n.sim,+ rand.gen = rg1, res = as.vector(res)))+ }>> lynx.3 <- tsboot(lynx.res, lynx.fun, R = 99, sim = "model", n.sim = 114,+ orig.t = FALSE, ran.gen = lynx.sim,+ ran.args = list(ts = log(lynx), model = lynx.model))>> # For "post-blackening" we need to define another function> lynx.black <- function(res, n.sim, ran.args) {+ ts.orig <- ran.args$ts+ ts.mod <- ran.args$model+ mean(ts.orig) + ts(arima.sim(model = ts.mod,n = n.sim,innov = res))+ }>> # Now we can run apply the two types of block resampling again but this> # time applying post-blackening.> lynx.1b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "fixed",+ n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,+ ran.args = list(ts = log(lynx), model = lynx.model))>> lynx.2b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "geom",+ n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,+ ran.args = list(ts = log(lynx), model = lynx.model))>> # To compare the observed order of the bootstrap replicates we> # proceed as follows.> table(lynx.1$t[, 1])2 3 4 5 7 8 10 11 12 13 1416 19 38 4 6 3 1 9 1 1 1> table(lynx.1b$t[, 1])2 3 4 5 6 7 8 11 12 14 156 2 22 6 4 6 3 40 7 1 2> table(lynx.2$t[, 1])2 3 4 5 6 7 8 10 11 1312 18 51 5 2 3 1 2 4 1> table(lynx.2b$t[, 1])2 3 4 5 6 7 8 9 10 11 12 13 15 212 1 21 4 1 10 4 1 3 45 3 1 2 1> table(lynx.3$t[, 1])2 3 4 5 6 7 8 9 10 11 12 13 14 154 8 11 2 1 4 2 2 2 54 6 1 1 1> # Notice that the post-blackened and model-based bootstraps preserve> # the true order of the model (11) in many more cases than the others.>>>> cleanEx()> nameEx("var.linear")> ### * var.linear>> flush(stderr()); flush(stdout())>> ### Name: var.linear> ### Title: Linear Variance Estimate> ### Aliases: var.linear> ### Keywords: nonparametric>> ### ** Examples>> # To estimate the variance of the ratio of means for the city data.> ratio <- function(d,w) sum(d$x * w)/sum(d$u * w)> var.linear(empinf(data = city, statistic = ratio))[1] 0.03248701>>>> ### * <FOOTER>> ###> options(digits = 7L)> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")Time elapsed: 29.023 0.078 29.298 0 0> grDevices::dev.off()null device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')