---
title: Simulation for Inference I --- The Bootstrap
author: 36-467/667
date: 22 October 2020 (Lecture 16)
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Prob}[1]{\mathbb{P}\left[ #1 \right]}
\newcommand{\Probwrt}[2]{\mathbb{P}_{#2}\left( #1 \right)}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\]

```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small", background="white",
               cache=TRUE, autodep=TRUE,
               tidy=FALSE, warning=FALSE, message=FALSE,
               echo=TRUE)
```


## In our last episode

- General strategy for simulating
- The Monte Carlo method: replace complicated probability calculations with repeated simulation
    + To find $\Expect{f(Y)}$, do $m$ simulations $Y_1, Y_2, ldots Y_m$
    + Then $m^{-1}\sum_{1}^{m}{f(Y_i)} \approx \Expect{f(Y)}$
    + Also works for variances, probabilities, etc.


## Agenda

> - The sampling distribution of a statistic tells us about statistical uncertainty (standard errors, biases, confidence sets)
> - The **bootstrap principle**: _approximate_ the sampling distribution by _simulating_ from a good model, and treating the simulation run just like the data
> - Sometimes we simulate from the model we're estimating (model-based or "parametric" bootstrap)
> - Sometimes we simulate by re-sampling the original data (resampling or "nonparametric" bootstrap)
> - Stronger assumptions $\Rightarrow$ less uncertainty _if we're right_


## Statistical Uncertainty

> - Re-run the experiment (survey, census, ...) and get different data
> - $\therefore$ everything we calculate from data (estimates, test statistics, policies, ...) will change from trial to trial as well
> - This variability is (the source of) **statistical uncertainty**
> - Quantifying this is a way  to be honest about what we actually know


## Measures of Uncertainty


> - Standard error = standard deviation of an estimator
>     + (could use median absolute deviation instead, etc.)
> - $p$-value = Probability we'd see a signal this big if there was just noise
> - Confidence region = All the parameter values we can't reject at low error rates:

> 1. _Either_ the true parameter is in the confidence region
> 2. _or_ we are very unlucky
>    + (Something happened whose probability was very low under _all_ parameter values)
> 3. _or_ our model is wrong

> - etc., etc.

## The Sampling Distribution Is the Source of All Knowledge


- Data $X \sim P_{X,\theta_0}$, for some true $\theta_0$
- We calculate a statistic $T = \tau(X)$ so it has distribution $P_{T,\theta_0}$
- If we knew $P_{T,\theta_0}$, we could calculate
    + $\Var{T}$ (and so standard error)
    + $\Expect{T}$ (and so bias)
    + quantiles (and so confidence intervals or $p$-values), etc.

## The Problems

> 1. Most of the time, $P_{X,\theta_0}$ is very complicated
> 2. Most of the time, $\tau$ is a very complicated function
> 3. We certainly don't know $\theta_0$

> - Upshot: We don't know $P_{T,\theta_0}$ and can't use it to calculate anything

## The Solutions

> - Classically ($\approx 1900$--$\approx 1975$): Restrict the model and the
statistic until you can calculate the sampling distribution, at least for very
large $n$
>     + This is the asymptotic theory we did in [Lecture 13](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/13/lecture-13.html) and [Lecture 14](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/14/lecture-14.html)
> - Modern ($\approx 1975$--): Use complex models and statistics, but _simulate_ calculating the statistic on the model

## The Bootstrap Principle

> 1. Find a good estimate $\hat{P}$ for $P_{X,\theta_0}$
> 2. Generate a simulation $\tilde{X}$ from $\hat{P}$; set $\tilde{T} = \tau(\tilde{X})$; repeat
> 3. Use the simulated distribution of the $\tilde{T}$ to approximate $P_{T,\theta_0}$

> - First step: find a good estimate $\hat{P}$ for $P_{X,\theta_0}$

## Model-based Bootstrap

If we are using a model, our best guess at $P_{X,\theta_0}$ is
$P_{X,\hat{\theta}}$, with our best estimate $\hat{\theta}$ of the parameters

#### The Model-based Bootstrap

- Get data $X$, estimate $\hat{\theta}$ from $X$
- Repeat $m$ times:
    + Simulate $\tilde{X}$ from $P_{X,\hat{\theta}}$ (simulate data of same size/"shape" as real data)
    + Calculate $\tilde{T} = \tau(\tilde{X}$) (treat simulated data the same as real data)
- Use the simulated distribution of $\tilde{T}$ as $P_{T,\theta_0}$


## Example: lynxes in Canada

![](http://farm4.static.flickr.com/3131/2295650288_a8f7a9a51b_m.jpg)

## Example: lynxes in Canada

- Annual time series of number of lynxes trapped for the fur company in the Mackenzie River area of Canada:

```{r}
data(lynx)
plot(lynx, xlab="Year", ylab="Lynxes", lwd=2)
```

## Fitting an AR model to the example data

- This oscillates, so an AR(1) isn't appropriate --- try an AR(2)

```{r}
lynx.ar2 <- ar.ols(lynx, order.max=2, aic=FALSE, demean=FALSE, intercept=TRUE)
lynx.ar2$ar
lynx.ar2$x.intercept
```

so the estimated model is
\[
X(t)  = `r lynx.ar2$x.intercept` + `r lynx.ar2$ar[1]` X(t-1) - `r abs(lynx.ar2$ar[2])` X(t-2) + \epsilon(t)
\]

## Simulating from the fitted model

- What do some time series from this look like?
    + (We wrote `ar2.sim` last time)

```{r, include=FALSE}
# Simulate a time series of length n from a Gaussian AR(2) model with given
# coefficients, noise level and initial condition
  # Copied from lecture 17
# Inputs: length of output time series (n)
  # intercept of AR(2) (a)
  # vecotr of coefficients of AR(2) (b)
  # standard deviation of the innovations (sd.innovation)
  # initial two values of time series (x.start)
# Output: vector of length n
# Presumes: innovations should be Gaussian white noise
  # a and sd.innovation are length-1 numerical values
  # b and x.start are length-2 numerical vectors
  # b has coefficient on X(t-1) first
  # x.start has first value first
ar2.sim <- function(n, a, b, sd.innovation, x.start) {
    x <- vector(length=n)
    x[1] <- x.start[1] # Because R enumerates from 0
    x[2] <- x.start[2]
    for (t in 1:(n-2)) {
        x[t+2] <- a+b[1]*x[t+1] + b[2]*x[t] + rnorm(1, sd=sd.innovation)
    }
    return(x)
}
```

```{r}
lynx.sims <- replicate(5, ar2.sim(n=length(lynx),
                                   a=lynx.ar2$x.intercept,
                                   b=lynx.ar2$ar,
                                   sd.innovation=sd(lynx.ar2$resid, na.rm=TRUE),
                                   x.start=lynx[1:2]))
plot(lynx, xlab="Year", ylab="Lynxes caught", lwd=2)
for (i in 1:ncol(lynx.sims)) {
    lines(1821:1934, lynx.sims[,i], col="grey", lty=i)
}
abline(h=0, col="blue", lwd=0.5)
```

## Uncertainty in the coefficients {.smaller}

- Write a function which estimates the coefficients

```{r}
# Estimate the coefficients of an AR(2) model
# Inputs: a time series (x)
# Outputs: vector containing intercept and 2 AR slopes
ar2.estimator <- function(x) {
    # Treat the input exactly the same way we treated the real data
    fit <- ar.ols(x, order.max=2, aic=FALSE, demean=FALSE, intercept=TRUE)
    return(c(a=fit$x.intercept, b1=fit$ar[1], b2=fit$ar[2]))
}
```

Check that this gives the right answer on the data, and a _different_ answer
on a simulation:
```{r}
ar2.estimator(lynx)
ar2.estimator(lynx.sims[,1])
```


## Uncertainty in the coefficients {.smaller}

Now do _lots_ of simulations, and apply this function to each simulation:

```{r}
lynx.coef.sims <- replicate(1000, ar2.estimator( ar2.sim(n=length(lynx),
                                   a=lynx.ar2$x.intercept,
                                   b=lynx.ar2$ar,
                                   sd.innovation=sd(lynx.ar2$resid, na.rm=TRUE),
                                   x.start=lynx[1:2])))
dim(lynx.coef.sims)
lynx.coef.sims[1:3, 1:5]
```

## Uncertainty in the coefficients {.smaller}

Standard errors $\approx$ standard deviation across simulations:
```{r}
apply(lynx.coef.sims, 1, sd)
```

95% confidence intervals $\approx$ quantiles across simulations:
```{r}
apply(lynx.coef.sims, 1, quantile, prob=c(0.025, 0.975))
```

## For Gaussian AR$(p)$ models...

- ... the bootstrapping isn't needed, because there the asymptotic theory we worked out before is actually exact
    + Basically: Homework 7, only with even more book-keeping
- But we can use this bootstrap approach for _any_ generative time-series model
    + E.g., make $X(t)$ have a Poisson distribution rather than a Gaussian
        * No negative lynx counts, only whole numbers of lynxes...
    + All we have to do is swap out the simulation

## Prediction

- You can apply these ideas when the "statistic" is $X(t)$ for a fixed $t$
- Or when it's $X(t) | X(t-1) = x$ for a fixed $x$
    + Since the model tells us how to calculate that!
- Expected values across simulations = best guess for the point prediction
- Variances and intervals incorporate parameter uncertainty


## Resampling

- Problem: Suppose we don't have a trust-worthy model
- Resource: We do have data, which tells us a lot about the distribution
- Solution: **Resampling**, treat the sample like a whole population

## The Resampling ("Nonparameteric") Bootstrap

- Get data $X$, containing $n$ samples, find $T=\tau(X)$
- Fix a **block length** $k$ and divide $X$ into (overlapping) blocks of length $k$
    + So first block is $(X(1), \ldots X(k))$, second block is $(X(2), \ldots X(k+1))$, down to block $n-k+1$, which is $(X(n-k+1), X(n-k+2), \ldots X(n))$
- Repeat $m$ times:
    + Generate $\tilde{X}$ by drawing $n/k$ blocks from $X$ _with replacement_ (resample the data)
    + Concatenate the blocks to make $\tilde{X}$
    + Calculate $\tilde{T} = \tau(\tilde{X})$ (treat simulated data the same as real data)
- Use simulated distribution of $\tilde{T}$ as $P_{T,\theta}$

## Some useful code for block resampling {.smaller}

```{r}
# Convert a time series into a data frame of lagged values
# Inputs: a time series (ts)
  # maximum lag to use (order)
  # whether older values go on the right or the left (right.older)
# Output: a data frame with (order+1) columns, named lag0, lag1, ... , and
  # length(ts)-order rows
design.matrix.from.ts <- function(ts,order,right.older=TRUE) {
  n <- length(ts)
  x <- ts[(order+1):n]
  for (lag in 1:order) {
    if (right.older) {
      x <- cbind(x,ts[(order+1-lag):(n-lag)])
    } else {
      x <- cbind(ts[(order+1-lag):(n-lag)],x)
    }
  }
  lag.names <- c("lag0",paste("lag",1:order,sep=""))
  if (right.older) {
    colnames(x) <- lag.names
  } else {
    colnames(x) <- rev(lag.names)
  }
  return(as.data.frame(x))
}
```



## Example: Back to the lynxes {.smaller}

```{r, include=FALSE}
# Source the PPW block-length estimator
source("http://www.math.ucsd.edu/~politis/SOFT/PPW/ppw.R")
# This estimates a good block-length for two different time-series bootstraps,
# one with randomized block lengths (the "stationary bootstrap"), and one
# with fixed block lengths (the "circular bootstrap"); we're using the
# second one
```

```{r}
head(lynx)
head(design.matrix.from.ts(lynx,2, right.older=FALSE))
```

## Example: Back to the lynxes

```{r}
# Simple block bootstrap
# Inputs: time series (ts), block length, length of output
# Output: one resampled time series
# Presumes: ts is a univariate time series
rblockboot <- function(ts,block.length,len.out=length(ts)) {
  # chop up ts into blocks
  the.blocks <- as.matrix(design.matrix.from.ts(ts,block.length-1,
    right.older=FALSE))
    # look carefully at design.matrix.from.ts to see why we need the -1
  # How many blocks is that?
  blocks.in.ts <- nrow(the.blocks)
  # Sanity-check
  stopifnot(blocks.in.ts == length(ts) - block.length+1)
  # How many blocks will we need (round up)?
  blocks.needed <- ceiling(len.out/block.length)
  # Sample blocks with replacement
  picked.blocks <- sample(1:blocks.in.ts,size=blocks.needed,replace=TRUE)
  # put the blocks in the randomly-selected order
  x <- the.blocks[picked.blocks,]
  # convert from a matrix to a vector and return
    # need transpose because R goes from vectors to matrices and back column by
    # column, not row by row
  x.vec <- as.vector(t(x))
    # Discard uneeded extra observations at the end silently
  return(x.vec[1:len.out])
}
```

```{r}
head(lynx)
head(rblockboot(lynx, 3))
```

## Resampling the lynxes

- We can do everything the same as before, just changing how we simulate

```{r}
lynx.coef.sims.resamp <- replicate(1000, ar2.estimator(rblockboot(lynx, 3)))
```

- Standard errors:
```{r}
apply(lynx.coef.sims.resamp, 1, sd)
```
- 95% confidence intervals:
```{r}
apply(lynx.coef.sims.resamp, 1, quantile, prob=c(0.025, 0.975))
```



## Why should this work?

- Preserves the distribution within blocks
- Independence across blocks
- So $\tilde{X} \approx$ stationary, with the same short-range distribution as $X$ but no longer-range dependence
- Should work well if $X$ is stationary without long-range dependence
- Want to let $k$ grow with $n$, but still have a growing number of blocks, so $n/k_n \rightarrow \infty$, or $k=o(n)$
- Some (complicated) theory suggests $k \propto n^{1/3}$
    + @PolitisWhite2004 provides an automatic procedure for picking $k$
       * Estimates bias (from too-short blocks) and variance (from too-few long blocks) and balances
    + R implementation at [http://www.math.ucsd.edu/~politis/SOFT/PPW/ppw.R]
    + ...which suggests a block-length of 3 for the lynxes

## What about non-stationary time series?

- Model-based: OK if the model is non-stationary
- Otherwise:
    + Estimate the trend
    + Bootstrap the deviations from the trend (by model or resampling)
    + Add these back to the trend
    + $\tilde{X} =$ trend + bootstrapped deviations


## What about spatial data?

- Model-based: simulate the model
- Resampling: use spatial blocks (e.g., rectangles)
    + @Levina-Bickel-resampling-random-fields suggests the best blocks aren't rectangles but rectangles with a notch cut out at one corner



## Sources of Error in Bootstrapping

> - **Simulation** Using only $m$ bootstrap replicates
>     + Make this small by letting $m\rightarrow\infty$
>     + Costs computing time
>     + Diminishing returns: error is generally $O(1/\sqrt{m})$
> - **Approximation** Using $\hat{P}$ instead of $P_{X,\theta_0}$
>     + Make this small by careful statistical modeling
> - **Estimation** Only a finite number of samples $n$
>     + Make this small by being careful about what we simulate (example in backup slides)

> - Generally: for fixed $n$, resampling boostrap shows more uncertainty than model-based bootstraps, but is less at risk to modeling mistakes (yet another bias-variance tradeoff)

## Summary

- Bootstrapping is applying the Monte Carlo idea to the sampling distribution
    + Find a good approximation to the data-generating distribution
    + Simulate it many times
    + Treat each simulation run just like the original data
    + Distribution over simulations $\approx$ sampling distribution
- Choice of how to simulate:
    + Based on a model
    + By resampling blocks
- Many, many refinements, but this is the core of it


## Backup: Some General Ways of Improving the Bootstrap

- Improving the initial estimate $\hat{P}$
    + Model-based vs. resampling is a bias/variance trade-off here
- Reducing the number of simulations or speeding them up
    + E.g., if we're estimating a mean, we'd like to get pairs of _negatively_ correlated samples; there are tricks for doing this ("antithetic sampling")
- Transforming $\tau$ so the final approximation is more stable
    + The distribution of $\tilde{\theta}$ approaches that of $\hat{\theta}$ as $n$ grows
    + But often the distribution of $\tilde{\theta}-\hat{\theta}$ approaches that of $\hat{\theta}-\theta_0$ more quickly ("pivoting"; see below)
    + And the distribution of $\frac{\tilde{\theta}-\hat{\theta}}{\sqrt{\Var{\tilde{\theta}}}}$ is an even better approximation to the distribution of $\frac{\hat{\theta}-\theta_0}{\sqrt{\Var{\hat{\theta}}}}$ ("studentizing")


## Backup: The "surrogate data" method

- You think that $X$ comes from some interesting class of processes $\mathcal{H}_1$
- Your skeptical friend thinks maybe it's just some boring process in  $\mathcal{H}_0$
- Fit a model to $X$ using only the boring processes in $\mathcal{H}_0$, say $\hat{Q}$
- Now simulate $\hat{Q}$ to get many "surrogate data" sets $\tilde{X}$
- Look at the distribution of  $\tau(\tilde{X})$ over your simulations
- If $\tau(X)$ is really unlikely under that simulation, it's evidence against $\mathcal{H}_0$
    + When would it be evidence _for_ $\mathcal{H}_1$?
- This is a test of the  hypothesis $\mathcal{H}_0$ against the alternative $\mathcal{H}_1$
- The term "surrogate data" comes out of physics [@Theiler-et-al-surrogate-data]
    + testing the null of "this is a linear stochastic process" vs. "this is a deterministic, but chaotic, nonlinear dynamical system"


## Backup: Incorporating uncertainty in the initial conditions

- Simulating a dynamical model needs an initial value $X(0)$ (or $X(1)$, or maybe $X(1), X(2)$, etc.)
- What if this isn't known precisely?
- Easy case: it follows a probability distribution
    + Then draw $X(0)$ from that distribution independently on each run
    + Aggregate across runs as usual
- Less easy: $X(0)$ within known limits but no distribution otherwise
    + My advice: draw $X(0)$ uniformly from inside the region
    + Rationale 1: Maximum entropy distribution compatible with constraints
    + Rationale 2: If the process is rapidly mixing, memory of the distribution of $X(0)$ decays quickly anyway


## Backup: Improving on the Crude Confidence Interval

- Crude CI use distribution of $\tilde{\theta}$ under $\hat{\theta}$
- But really we want the distribution of $\hat{\theta}$ under $\theta$
- Mathematical observation: Generally speaking, (distributions of)
$\tilde{\theta} - \hat{\theta}$ is closer to $\hat{\theta}-\theta_0$ than the distribution of 
$\tilde{\theta}$ is to $\hat{\theta}$
   + ("centering" or "pivoting")


## Backup: The Basic, Pivotal CI

quantiles of $\tilde{\theta}$ = $q_{\alpha/2}, q_{1-\alpha/2}$

\[
\begin{eqnarray*}
1-\alpha & = & \Probwrt{q_{\alpha/2} \leq \tilde{\theta} \leq q_{1-\alpha/2}}{\hat{\theta}} \\
& = & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \tilde{\theta} - \hat{\theta} \leq q_{1-\alpha/2} - \hat{\theta}}{\hat{\theta}} \\
 & \approx & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \hat{\theta} - \theta_0 \leq q_{1-\alpha/2} - \hat{\theta}}{\theta_0}\\
& = & \Probwrt{q_{\alpha/2} - 2\hat{\theta} \leq -\theta_0 \leq q_{1-\alpha/2} - 2\hat{\theta}}{\theta_0}\\
& = & \Probwrt{2\hat{\theta} - q_{1-\alpha/2} \leq \theta_0 \leq 2\hat{\theta}-q_{\alpha/2}}{\theta_0}
\end{eqnarray*}
\]

Basically: re-center the simulations around the empirical data


```{r}
2*ar2.estimator(lynx)-apply(lynx.coef.sims.resamp, 1, quantile, prob=0.975)
2*ar2.estimator(lynx)-apply(lynx.coef.sims.resamp, 1, quantile, prob=0.025)
```

## Backup: Further reading

- Bootstrap begins with @Efron-on-boostrap-in-annals (drawing on an older idea from the 1940s called the "jackknife")
- Block bootstrap for time series comes from @Kunsch-bootstrap-for-stationary-obs
- @Buhlmann-on-bootstraps-for-time-series is a good review on time series bootstraps
- @Davison-Hinkley-bootstrap is the best over-all textbook on the bootstrap
- @Lahiri-resampling-for-dependent is _very_ theoretical



## References
