---
title: "Bootstrap"
author: "JJB + Course"
date: "02/25/2019"
output:
   html_document:
     toc: true
     toc_float:
       collapsed: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Non-parametric Bootstrap

## Example: Bootstrapping subject-heights mean

```{r bootstrap-heights, cache = TRUE}
# Step 1: Obtain samples from population
sample_data = data.frame(id     = c(1, 2, 3, 4, 5),
                         sex    = c("M","F","F","M","M"),
                         height = c(6.1, 5.5, 5.2, 5.6, 5.9))                            
theta_hat = mean(sample_data$height)   # Compute the mean for the height data
n_obs = nrow(sample_data)              # Length of data
boot_iter = 500L                       # Number of bootstrap iterations
theta_star = rep(NA, boot_iter)        # Bootstrapped estimate of theta

for (i in seq_len(boot_iter)) {
  set.seed(11882 + i)                  # Set seed for reproducibility
  # Step 2: Randomly sample observations positions from 1 to n_obs
  indexes = sample(n_obs, n_obs, replace = TRUE)

  # Extract out the observation positions
  sample_data_star = sample_data[indexes,, drop = FALSE]

  # Step 3: Compute the desired statistic on the bootstrapped values
  theta_star[i] = mean(sample_data_star$height)
  
  
} # Step 4: Repeat until i matches boot_iter 

# See first portion of output
head(theta_star)

# See compute the mean of the bootstrapped distribution
mean(theta_star)

# Compare with the mean of the population sample
theta_hat
```

## Examples: Compute Percentiles

```{r example-quantiles}
# Sample Data
x = c(1, 2, 3, 4, 5, 6)

quantile(x, 
  probs = c(0.25, 0.5, 0.75, 1)
)

# Median is the 50% quantile
median(x)
```

### Example: Distribution Quantiles

```{r my-dist-quantile}
# Retrieve real quantiles
qnorm(c(0.25, 0.5, 0.75, 1))
```


### Example: Quantile CI

```{r ci-bootstrap-quantile, cache = TRUE, dependson="bootstrap-heights"}
alpha    = 0.05
alpha / 2
1 - alpha / 2
ci_range = quantile(theta_star, probs = c(alpha / 2, 1 - alpha / 2))
ci_range
```


## Example: Plotting Bootstrapped Samples

```{r view-bootstrap-samples}
# Graph results
library("ggplot2")

graph_bootstrap = data.frame(iter       = seq_along(theta_star),
                             theta_star = theta_star)
ggplot(graph_bootstrap) + 
  geom_histogram(aes(theta_star)) + labs(
    title = "Non-parametric Bootstrapped Data",
    sub   = "Example simulation",
    x     = "Values of Theta^*",
    y     = "Frequency of Theta^* Values"
  )
```


## Exercise: Standard Deviation of `iris`' `Sepal.Width`

```{r check-out-iris}
# Data set is included with base R
head(iris)

# Step 1: Obtain samples from population
sample_data = iris

# Extract all of the Sepal.Width observations
# grab observations
```

```{r view-bootstrap}
# Graph results
library("ggplot2")

graph_bootstrap = data.frame(iter       = seq_along(theta_star),
                             theta_star = theta_star)
ggplot(graph_bootstrap) + 
  geom_histogram(aes(theta_star)) + labs(
    title = "Non-parametric Bootstrapped Data",
    sub   = "Example simulation",
    x     = "Values of Theta^*",
    y     = "Frequency of Theta^* Values"
  )
```



## Example: Non-parametric bootstrap with `boot`

```{r r-boot-pkg, eval = FALSE}
# install.packages('boot')
library('boot')

# Create a data.frame w/ distance data
nsim_obs = 100

# Set seed to reproduce generated data
set.seed(981)

problem_data = 
  data.frame(distance = rchisq(nsim_obs, df = 5)
             # other data here 
)

# Create a sampling or subset function for the data
sampling_function = function(d, ind) {
  # Extract out the observation positions (this allows for multiple indices
  # to be selected)
  problem_data_star = problem_data[ind,, drop = FALSE]
  
  # Compute the desired statistic on the bootstrapped data
  mean(problem_data_star$distance)
}

# Run the bootstrapping procedure
booted_means = boot(
  data = problem_data,            # Pass the data
  statistic = sampling_function,  # Pass a function that computes a sample and statistic
  R = 200                         # Number of iterations
)

# Calculate different confidence intervals
boot.ci(booted_means)

# Show underlying statistic distribution (e.g. t is the mean statistic)
plot(booted_means)
```


# Parametric Bootstrap

## Example: Parametric bootstrap with `sd` and `mean` of a Normal Distribution

```{r r-parametric, cache = TRUE}
sample_values = rnorm(1000)           # Step 1: Obtain samples from known
                                      # population distribution.
                                      # Step 2: Obtain statistics
theta_mean_hat = mean(sample_values)  # Compute sample mean
theta_sd_hat = sd(sample_values)      # Compute sample standard deviation

n_obs = length(sample_values)         # Length of data
boot_iter = 250L                      # Number of bootstrap iterations
theta_mean_star = rep(NA, boot_iter)  # Bootstrapped estimate of mean
theta_sd_star = rep(NA, boot_iter)    # Bootstrapped estimate of standard dev

for (i in seq_len(boot_iter)) {
  set.seed(385 + i)                   # Set seed for reproducibility

  # Step 3: Randomly generate observations under distribution
  sample_values_star = rnorm(n_obs, mean = theta_mean_hat, sd = theta_sd_hat )

  # Step 4: Compute the desired statistic on the bootstrapped values
  theta_mean_star[i] = mean(sample_values_star)
  theta_sd_star[i] = sd(sample_values_star)

} # Step 5: Repeat until i matches boot_iter 

mean(theta_mean_star)
theta_mean_hat
mean(theta_sd_star)
theta_sd_hat
```

```{r my-data}
# install.packages("tidyr")
library("tidyr")

graph_bootstrap = data.frame(iter            = seq_len(boot_iter),
                             theta_mean_star = theta_mean_star,
                             theta_sd_star   = theta_sd_star)

# Tidy data
tidy_bootstrap = gather(graph_bootstrap,
                        key = theta_type_estimated, value = theta_star,
                        theta_mean_star:theta_sd_star)

theta_vlines = data.frame(theta_type_estimated = c("theta_mean_star", "theta_sd_star"), 
                          theta_val = c(0, 1))

library("ggplot2")
ggplot(tidy_bootstrap) + 
  geom_histogram(aes(theta_star)) +
  geom_vline(data = theta_vlines,
             aes(xintercept = theta_val, color = theta_type_estimated)) +
  facet_wrap(~theta_type_estimated) + 
  labs(
    title = "Non-parametric Bootstrapped Data",
    color = "Type of Theta Estimated",
    sub   = "Example simulation",
    x     = "Values of Theta^*",
    y     = "Frequency of Theta^* Values"
  ) + theme_bw()
```

## Exercise: Parametric bootstrap with `mean`, `sd`, and `median` of a Poisson Distribution

Estimate the `mean`, `sd`, and `median` of a Lambda distribution with an initial
parameter of `lambda = 3`.

```{r r-parametric-pois, cache = TRUE, eval = FALSE}

```


```{r overview-data, eval = FALSE}


```

