---
title: "Lecture 24: Random Simulation, For Fun And Profit"
output: ioslides_presentation
---

## Previously...

Topics we got to earlier in the course:

- How to make draws from, and use random distributions
- Writing functions with known run-times and verifiable results
- Why not to use `for` loops


## Today...

- Writing simulations using random number generation
- From the jackknife to the bootstrap
- Processes that have memory
- Finally, `for` loops!


## `[dpqr]unif`, etc.

Recall how we drew from various distributions: `runif()`: draws from the uniform distribution (others follow)

Home-baked discrete distributions: Use `sample()`

```{r}
population.values <- 1:3
probs <- c(5,2,3)
my.draw <- sample (population.values, 100000, probs, replace=TRUE)
table(my.draw)
```

## Permutations with `sample()`

`sample()` is powerful -- it works on any object that has a defined `length()`. Permutations:

```{r}
sample(5)
sample(1:6)
replicate(3,sample(c("Curly","Larry","Moe","Shemp")))
sample(list("A",3,sum))
```

## Resampling with `sample()`

Resampling from any existing distribution gets us the **bootstrap** family of estimators:

```{r}
bootstrap.resample <- function (object) sample (object, length(object), replace=TRUE)
replicate(5, bootstrap.resample (6:10))
```

Recall: the *jackknife* worked by removing one point from the sample and recalculating the statistic of interest. Here we resample the same length with replacement.

## Bootstrap test

The 2-sample `t`-test checks for differences in means according to a known null distribution. Let's resample and generate the sampling distribution under the bootstrap assumption:

```{r}
library(MASS)
diff.in.means <- function(df) {
  mean(df[df$Sex=="M","Hwt"]) - mean(df[df$Sex=="F","Hwt"])
}
resample.diffs <- replicate(1000, diff.in.means(cats[bootstrap.resample(1:nrow(cats)),]))
hist(resample.diffs); abline(v=diff.in.means(cats), col=2, lwd=3)

```

## Dependence

Most interesting examples in probability have a little dependence added in: "If it rained yesterday, what is the probability it rains today?"

Use this to generate weather patterns and probabilities for some time in the future. Almost always occurs with time series; can occur with other dependence (spatial -- if it's sunny in Dallas today, will it also be sunny in Fort Worth?)

## Markov Dependence

Suppose we have a sequence of observations that are dependent. In a time series, what happens next depends on what happened before:

\[ p(X_1, X_2, ..., X_n) = p(X_1)p(X_2|X_1)...p(X_n|X_{n-1},...,X_1) \]

(Note: you could, of course, condition on the future to predict the past, if you had a time machine.)

Markov dependence: each outcome only depends on the one that came before.

\[ p(X_1, X_2, ..., X_n) = p(X_1)\prod_{i=2}^np(X_i|X_{i-1})  \]

## Generating a Markov Chain

1. Set up the conditional distribution.

2. Draw the initial state of the chain.

3. For every additional draw, use the previous draw to inform the new one.


## And `for` loops are back in play!

Alluded to but not taught explicitly because we had better ways of doing things in R -- until now.

Simple weather model: 
- if it was sunny yesterday, today's chance of sun is 80%.
- if it wasn't sunny yesterday, today's chance of sun is 20%.

Simulate for one year (at the equator?)

```{r}
sunny.year <- rep(NA, 365)
sunny.year[1] <- 1
for (day in 2:365) sunny.year[day] <- rbinom(1,1,0.2 + 0.6*sunny.year[day-1])
```

## And `for` loops are back in play!

```{r}
plot(sunny.year, main="Sunny Days in An Equatorial City", xlab="Day", ylab="Sunshine?", ty="l")
```

## Different From Independence


```{r}
boring.year <- rbinom(365, 1, 0.5)
plot(boring.year, main="Sunny Days in A Coin Flip City", xlab="Day", ylab="Sunshine?", ty="l")
```

## The above chain

Transitions are represented as a matrix: $Q_{ij}$ is $P(X_t = j|X_{t-1} = i)$.

```{r echo=FALSE}
matrix (c(0.8, 0.2, 0.2, 0.8), nrow=2)
```

## Flippy chain

```{r}
weird.year <- rep(NA, 365)
weird.year[1] <- 1
transition.matrix <- matrix (c(0.2, 0.8, 0.8, 0.2), nrow=2)
for (day in 2:365) weird.year[day] <- sample(1:2, 1, prob=transition.matrix[weird.year[day-1],])
```

## Flippy chain

```{r}
plot(weird.year, main="Sunny Days in Al Gore's Nightmare", xlab="Day", ylab="Sunshine?", ty="l")
```

## General Markov Chain

```{r}
rmarkovchain <- function (nn, 
                          transition.matrix, 
                          start=sample(1:nrow(transition.matrix), 1)) {
  output <- rep (NA, nn)
  output[1] <- start
  for (day in 2:nn) output[day] <- 
    sample(ncol(transition.matrix), 1, 
           prob=transition.matrix[output[day-1],])
}
```

## Simple Unbounded Markov Chain

"Unbiased Random Walk": Independent events atop a dependent structure.

```{r}
randomwalk <- function (nn, upprob=0.5, start=50) {
  output <- rep (NA, nn)
  output[1] <- start
  for (iteration in 2:nn) 
    output[iteration] <- 
      output[iteration-1] + 2*rbinom(1, 1, upprob) - 1
  output  
}
```

## Simple Unbounded Markov Chain

```{r}
plot (randomwalk(10000, start=200), main="Simple Random Walk")
```

## Simple Unbounded Markov Chain

```{r}
plot (randomwalk(10000, start=200), main="Simple Random Walk")
```

## Simple Unbounded Markov Chain

```{r}
plot (randomwalk(10000, start=200), main="Simple Random Walk")
```

## Simple Unbounded Markov Chain

```{r}
plot (randomwalk(10000, start=200), main="Simple Random Walk")
```