---
title: "Iteration and Simulation"
author: "Statistical Computing, 36-350"
date: "Friday October 7, 2016"
---

Drug effect example
===

- Let's start with a motivating example: suppose we had a model for the way a drug affected certain patients
- E.g., all patients will undergo chemotherapy. We believe those who aren't given the drug experience a reduction in tumor size of percentage
$$
X_{\mathrm{no\,drug}} \sim 100 \cdot \mathrm{Exp}(\mathrm{mean}=R), 
\;\;\; R \sim \mathrm{Unif}(0,1) 
$$
- And those who were given the drug experience a reduction in tumor size of percentage
$$
X_{\mathrm{drug}} \sim 100 \cdot \mathrm{Exp}(\mathrm{mean}=2)
$$
(Here $\mathrm{Exp}$ denotes the exponential distribution, and $\mathrm{Unif}$ the uniform distribution)

(Continued)
===

```{r}
# Simulate from model, supposing 60 subjects in each group 
set.seed(0)
n = 60 
mu.drug = 2
mu.nodrug = runif(n, min=0, max=1)
x.drug = 100*rexp(n, rate=1/mu.drug) 
x.nodrug = 100*rexp(n, rate=1/mu.nodrug)
```

(Continued)
===

```{r}
# Find the range of all the measurements together, and define breaks
x.range = range(c(x.nodrug,x.drug))
breaks = seq(min(x.range),max(x.range),length=20)
# Produce hist of the non drug measurements, then drug measurements on top
hist(x.nodrug, breaks=breaks, probability=TRUE, xlim=x.range, 
     col="lightgray", xlab="Percentage reduction in tumor size", 
     main="Comparison of tumor reduction")
# Plot a histogram of the drug measurements, on top
hist(x.drug, breaks=breaks, probability=TRUE, col=rgb(1,0,0,0.2), add=TRUE) 
# Draw estimated densities on top, for each dist
lines(density(x.nodrug), lwd=3, col=1)
lines(density(x.drug), lwd=3, col=2)
legend("topright", legend=c("No drug","Drug"), lty=1, lwd=3, col=1:2)
```

Repetition and reproducibility
===

- One single simulation is not always trustworthy (depends on the situation, of course)
- In general, simulations should be repeated and aggregate results reported---requires iteration!
- To make random number draws reproducible, we must set the seed with `set.seed()`
- More than this, to make simulation results reproducible, we need to follow **good programming practices**
- Gold standard: any time you show a simulation result (a figure, a table, etc.), you have code that can be run (by anyone) to produce exactly the same result

Iteration and simulation (and functions): good friends
===

- Writing a **function** to complete a single run of your simulation is often very helpful
- This allows the simulation itself to be intricate (e.g., intricate steps, several simulation parameters), but makes running the simulation simple
- Then you can use **iteration** to run your simulation over and over again
- Good design practice: write another function for this last part (running your simulation many times)

Code sketch
===

Consider the code below for a generic simulation. (Think about how you would frame this for the drug effect example)

```
# Function to do one simulation run
one.sim = function(param1, param2=value2, param3=value3) {
  # Possibly intricate simulation code goes here
}

# Function to do repeated simulation runs
rep.sim = function(nreps, param1, param2=value2, param3=value3, seed=NULL) {
  # Set the seed, if we need to
  (if !is.null(seed)) set.seed(seed)
  
  # Run the simulation over and over
  sim.objs = vector(length=nreps, mode="list")
  for (r in 1:nreps) {
    sim.objs[r] = one.sim(param1, param2, param3)
  }
  
  # Aggregate the results somehow, and then return something
}
```

Saving results
===

Sometimes simulations take a long time to run, and we want to save intermediate or final output, for quick reference later

There two different ways of saving things from R (there are more than two, but here are two useful ones):

- `saveRDS()`: allows us to save single R objects (like a vector, matrix, list, etc.), in (say) .rds format. E.g.,

```
saveRDS(my.mat, file="my.matrix.rds")
```

- `save()`: allows us to save any number of R objects in (say) .rdata format. E.g., 

```
save(mat.x, mat.y, list.z, file="my.objects.rdata")
```

Note: there is a big difference between how these two treat variable names

Loading results
===

Corresponding to the two different ways of saving, we have two ways of loading things into R:

- `readRDS()`: allows us to load an object that has been saved by `savedRDS()`, and **assign a new variable name**, e.g.,

```{r, error=FALSE}
ls() # Show all the variables in our workspace, none so far
my.new.mat = readRDS("my.matrix.rds")
dim(my.new.mat)
my.new.mat[1:5,1:5]
```

- `load()`: allows us to load all objects that have been saved through `save()`, **according to their original variables names**, e.g.,

```{r, error=FALSE}
ls() # Show all the variables in our workspace, only my.new.mat
load("my.objects.rdata")
ls() # Show all the variables in our workspace, some are new
```