## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----echo = FALSE-------------------------------------------------------------
set.seed(142)

## -----------------------------------------------------------------------------
cor_boot <- function(data, row_numbers, method = "pearson")
{ 
    # Obtain the bootstrap sample:
    sample <- data[row_numbers,]
    
    # Compute and return the statistic for the bootstrap sample:
    return(cor(sample[[1]], sample[[2]], method = method))
}

## -----------------------------------------------------------------------------
mtcars |> plot(hp ~ drat, data = _)

## -----------------------------------------------------------------------------
library(boot)
boot_res <- mtcars |> 
              subset(select = c(hp, drat)) |> 
              boot(statistic = cor_boot,
                   R = 999)

## -----------------------------------------------------------------------------
plot(boot_res)

## -----------------------------------------------------------------------------
boot.ci(boot_res)

## -----------------------------------------------------------------------------
# Compute the bootstrap p-value based on the percentile interval:
library(boot.pval)
boot.pval(boot_res, type = "perc")

## -----------------------------------------------------------------------------
cor_boot_student <- function(data, row_numbers, method = "pearson")
{ 
    sample <- data[row_numbers,]
    
    correlation <- cor(sample[[1]], sample[[2]], method = method)
    
    inner_boot <- boot(sample, cor_boot, 100)
    variance <- var(inner_boot$t)

    return(c(correlation, variance))
}

# Run the resampling:
boot_res <- mtcars |> 
              subset(select = c(hp, drat)) |> 
              boot(cor_boot_student, 999)

# Compute the bootstrap p-value based on the studentized interval:
boot.pval(boot_res, type = "stud")

## ----message = FALSE----------------------------------------------------------
# Hypothesis test for the city data
# H0: ratio = 1.4
ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.pval(city.boot, theta_null = 1.4)

# Studentized test 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.pval(grav1.boot, type = "stud", theta_null = 0)