---
title: "Markov Chain Monte Carlo"
author: "Prof. Richey and Prof. Wright"
date: "April 27, 2018"
output: pdf_document
header-includes:
- \usepackage{tikz}
- \usetikzlibrary{arrows,automata}
- \usepackage{amsmath}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# Introduction
We will build a Markov chain on a set of $n=4$ states with a given steady-state distribution.
Imagine the states lined up horizontally in a row, as in the diagram below.
Our goal is achieve a steady-state distribution defined by the frequencies $f(i) = i$, for $i=1,2,3,4$.
To start, define a simple proposal transition function $Q$ as given by the arrows in the diagram:
\begin{center}
\begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=3cm,thick]
\tikzstyle{every state}=[fill=cyan,draw=none,text=white,font=\bfseries,minimum size=1.2cm]
\tikzset{every edge/.append style={font=\small}}
\node[state] (A) {1};
\node[state] (B) [right of=A] {2};
\node[state] (C) [right of=B] {3};
\node[state] (D) [right of=C] {4};
\path
(A) edge [out=45+90,in=135+90,distance=60] node [pos=.5,left] {$\frac{1}{2}$} (A)
(A) edge [bend left] node [pos=.5,above] {$\frac{1}{2}$} (B)
(B) edge [bend left] node [pos=.5,below] {$\frac{1}{2}$} (A)
(B) edge [bend left] node [pos=.5,above] {$\frac{1}{2}$} (C)
(C) edge [bend left] node [pos=.5,below] {$\frac{1}{2}$}(B)
(C) edge [bend left] node [pos=.5,above] {$\frac{1}{2}$} (D)
(D) edge [bend left] node [pos=.5,below] {$\frac{1}{2}$} (C)
(D) edge [out=45-90,in=135-90,distance=60] node [pos=.5,right] {$\frac{1}{2}$} (D);
\end{tikzpicture}
\end{center}
## We can do this
Use the code from last time to build the Markov chain transition Matrix $T$ with steady-state vector $p_{ss}$ proportional to $\left[ f(1),f(2),f(3),f(4)\right]$.
First, make a proposal transition matrix $Q$, and then modify it to obtain $T$.
## Agents again
Imagine an agent $A$ located at any state $j$.
Our method for constructing $T$ provides a way to put $A$ into motion.
At state $j$, use $Q$ to select a site $i$ to *consider* moving to.
In other words, select $i$ with probability $Q_{i,j}$.
* If $p_{ss}[j] \le p_{ss}[i]$, then transition $j \rightarrow i$.
* If $p_{ss}[j] > p_{ss}[i]$, then transition $j \rightarrow i$ **with probability $\rho$**, where $\rho$ is defined by
$$\rho = \frac{p_{ss}[i]}{p_{ss}[j]}.$$
### Computational note
Here a couple thoughts on how to do some of this computationally.
**Selecting the Propoal Transition:**
You could, of course, use the proposal transition matrix $Q$ to do this. However, this approach doesn't scale well. Imagine building the $n\times n$ matrix $Q$ when the number of states is $n = 10000$ or so.
Instead, let's take advantage of the simple transition rule. Suppose the agent is in state $j$.
* If the state $j$ is "interior", i.e., $1 < j < n$, then, with equal probability the agent transitions to a new state $i = j + \Delta$ where $\Delta = \pm 1$.
* If the agent is in state $j$ at the border, i.e., $j = 1$ or $j = n$, then the agent either stays put or moves to the only adjacent state with equal probability. Thus, $i = j + \Delta$ where:
+ If $j=1$, then $\Delta \in \{0,1\}$.
+ If $j=n$, then $\Delta \in \{-1,0\}$.
Hence, to select the proposal transition, we just need to select a value of $\Delta$ for a set of two values defined by the state. Specifically, given a current state, `currState`, you can create the proposed next state, `propState`, via:
```{r}
n <- 4
currState <- 1
if(1 < currState && currState < n){
Delta <- c(-1,1)
}else if(currState == 1){
Delta <- c(0,1)
}else{
Delta <- c(-1,0)
}
propState <- currState + sample(Delta, 1)
```
**Do something with probability $\rho$:**
This is simple. If you have a probability $\rho \in [0,1]$, then you can decide (True or False) to do something with this probability pretty easily using the R `sample` function
```{r}
rho <- .55
sample(c(TRUE,FALSE), 1, prob=c(rho,1-rho))
```
For example, to do something repeatedly with probability $\rho$:
```{r}
for(i in 1:10){
print(sample(c(TRUE,FALSE), 1, prob=c(rho,1-rho)))
}
```
#### Simple example.
Suppose $a=1$ and $b=1$. With probability $\rho=a/(a+b)$, let $a=a+1$. Otherwise, let $b=b+1$. What happens to the ratio $a/b$ if you do this, say, 10000 times?
```{r}
a <- 1
b <- 1
M <- 10000
for(m in 1:M){
rho <- a/(a+b)
doIt <- sample(c(TRUE,FALSE), 1, prob=c(rho,1-rho))
if(doIt){
a <- a + 1
}else{
b <- b + 1
}
}
c(a, b, a/(a+b))
```
Interesting. What would happen if you repeat this experiment a large number of times a keep track of the various values of the final ratio $a/(a+b)$?
There are other ways to make a decision based on a probability $\rho$ using the built-in random number generators in R. For example, the following code uses the `runif` (random uniform) function.
```{r}
rho <- .55
runif(1) < rho
```
This is simpler than what we had previously.
## Back to Markov chains
**Exercise:** With our simple 4-state Markov chain, start an agent $A$ in any state (state 1 is fine). Use the transition rules defined above to send $A$ off on a long journey. Keep track of where the agent goes, say in a vector called `statesVisited` . Does the agent visit the states according to the steady state frequencies?
Starter code:
Setup the scenario.
```{r}
n <- 4
f <- function(i){ i }
currState <- 1
```
Now we implement a move for the current state to a new state. First the proposed state.
```{r}
currState <- 1
if(1 < currState && currState < n){
Delta <- c(-1,1)
}else if(currState == 1){
Delta <- c(0,1)
}else{
Delta <- c(-1,0)
}
propState <- currState + sample(Delta, 1)
```
Now the decision based on the probabilities.
```{r}
currProb <- f(currState)
propProb <- f(propState)
if(currProb <= propProb){
currState <- propState
}else{
rho <- propProb/currProb
doMove <- sample(c(TRUE,FALSE), 1, prob=c(rho,1-rho))
#doMove <- runif(1)