--- 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") ```