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