---
title: "Lecture 25"
author: "36-350"
date: "11/18/2014"
output: ioslides_presentation
---

## Previously...

We introduced: Markov Chains, Monte Carlo simulation, Bootstrap method

## Today: Applications to Probability

Where the Bayesian revolution takes over: drawing from complex probability distributions

- rejection sampling
- The Metropolis algorithm
- General introduction to "Markov Chain Monte Carlo" 

**Optional reading**: Geyer, ``Practical Markov Chain Monte Carlo'', _Statistical Science_ **7** (1992): 473--483 

## One-dimension probability distributions

We have [dpqr]{norm, binom, pois, ...} for known distributions. For many good reasons, we'd like to be able to: 

- draw from complex distributions
- find the solution to relevant integrals
- eventually, go to many higher dimensions

What's available?

## Integration by Simulation

Law of large numbers: if $X_1, X_2, \ldots X_n$ all IID with p.d.f. $p$,

  \[
  \frac{1}{n}\sum_{i=1}^{n}{f(X_i)} \rightarrow \mathbb{E}_{p}[f(X)] =
  \int{f(x) p(x) dx}
  \]

The **Monte Carlo principle**: to find $\int{g(x) dx}$, draw from $p$ and take the sample mean of $f(x) = g(x)/p(x)$.


## And it works, too

Central limit theorem:

  \[
  \frac{1}{n}\sum_{i=1}^{n}{\frac{g(x_i)}{p(x_i)}} \rightsquigarrow
  \mathcal{N}\left(\int{g(x)dx}, \frac{\sigma^2_{g/p}}{n}\right)
  \]

The Monte Carlo approximation to the integral is unbiased, with root mean square error $\propto n^{-1/2}$ -- just keep taking Monte Carlo draws, and the error gets as small as you like, even if $g$ or $x$ are very complicated.


## Taking unknown samples

Generating from $p$ is easy if it's a standard distribution or we have a nice, invertible CDF (quantile method). What can we do if all we've got is the probability density function $p$?

Suppose the pdf $f$ is zero outside an interval $[c,d]$, and $\leq M$ on the interval.

```{r echo=FALSE}
plot(c(0,1), c(0,3), ty="n", main="A Sample Distribution", ylab="Density f(x)", xlab="x")
curve (dbeta(x, 3, 6), add=TRUE)
lines(c(0,0,1,1), c(0,3,3,0))
```

## Rejection sampling

We know how to draw from uniform distributions in any dimension. Do it in two:

```{r}
x1 <- runif(300, 0, 1); y1 <- runif(300, 0, 2.6); 
selected <- y1 < dbeta(x1, 3, 6)
```

```{r echo=FALSE}
plot(c(0,1), c(0,3), ty="n", main="A Sample Distribution", 
     ylab="Density f(x)", xlab="x")
curve (dbeta(x, 3, 6), add=TRUE)
lines(c(0,0,1,1), c(0,3,3,0))
points (x1, y1, col=1+1*selected, cex=0.1)
```

## Rejection sampling

```{r}
mean(selected)
accepted.points <- x1[selected]
mean(accepted.points < 0.5)
pbeta(0.5, 3, 6)
```

## Rejection sampling

For this to work efficiently: we have to cover the target distribution with one that sits close to it. 

```{r}
x2 <- runif(100000, 0, 1); y2 <- runif(100000, 0, 10); 
selected <- y2 < dbeta(x2, 3, 6)
mean(selected)
```

## Rejection sampling

```{r echo=FALSE}
plot(c(0,1), c(0,6), ty="n", main="A Sample Distribution", 
     ylab="Density f(x)", xlab="x")
curve (dbeta(x, 3, 6), add=TRUE)
lines(c(0,0,1,1), c(0,3,3,0))
points (x2, y2, col=1+1*selected, cex=0.1)
```

## Metropolis Algorithm

(Really: Metropolis, Rosenbluth x2, Teller x2, 1953)

Rejection sampling works if we know the complete distribution $p(x)$, with $\int p(x) = 1$. For many, many problems, we know the function up to the normalizing constant.

Example: Equations of the canonical ensemble of energy formations in the state of uncontrolled hydrogen fusion,

\[ p(E) = \frac{1}{|Z|} exp(f(E)/kT) \]

$|Z|$ is legendarily hard to calculate precisely.

## Basic Metropolis

On integer steps:

```{r}
metropolis.one <- function (xx, fun=dpois, ...) {
  prop <- 2*rbinom(1,1,0.5) - 1 + xx
  acc.ratio <- fun(prop, ...)/fun(xx, ...)
  output <- if (acc.ratio > runif(1)) prop else xx
  output
}
replicate (10, metropolis.one(10, lambda=5))
```


## Basic Metropolis

On integer steps:

```{r}
start <- 50
draws <- rep(NA, 10000)
for (iteration in 1:10000) {
  start <- metropolis.one (start, lambda=10)
  draws[iteration] <- start
}
```

## Basic Metropolis

```{r}
plot(draws)
```



## Basic Metropolis

Have to discard the initial state for "burn-in"

```{r}
plot(draws[-(1:100)])
```



## Basic Metropolis

Have to discard the initial state for "burn-in"

```{r}
hist(draws[-(1:100)])
hist(rpois(10000, 10), border=2, add=TRUE)
```



## Basic Metropolis

On integer steps:

```{r}
metropolis.cts <- function (xx, fun=dgamma, ...) {
  prop <- rnorm(1, xx, 1)
  acc.ratio <- fun(prop, ...)/fun(xx, ...)
  output <- if (acc.ratio > runif(1)) prop else xx
  output
}
replicate (20, metropolis.cts(20, shape=100, rate=10))
```


## Basic Metropolis

```{r}
start <- 50
draws <- rep(NA, 100000)
for (iteration in 1:100000) {
  start <- metropolis.one (start, fun=dgamma, shape=100, rate=10)
  draws[iteration] <- start
}
```

## Basic Metropolis

```{r}
plot(draws)
```



## Basic Metropolis

Have to discard the initial state for "burn-in"

```{r}
plot(draws[-(1:200)])
```



## Basic Metropolis

Have to discard the initial state for "burn-in"

```{r}
hist(draws[-(1:200)], prob=TRUE)
hist(rgamma(10000, 100, 10), border=2, add=TRUE, prob=TRUE)
```




