---
title: "Markov Chain Model of Diffusion"
author: "Math 365 Tanya Leise"
date: "January 1, 2018"
output: html_document
---

The objective of these exercises is to explore large-time behavior and equilibria (invariant probability distributions) of finite-state Markov chains. Feel free to discuss problems with each other during lab in addition to asking me questions. Please email your completed work as a pdf knit from your updated R markdown file, to tleise@amherst.edu.

To take powers of matrices in R, remember to load the matrix exponentiation package expm. If you haven't installed it yet, do _install.packages('expm')_.

```{r results='hide', message=FALSE, warning=FALSE}
library(expm)
```

### Warm-up problem

Suppose an urn contains 2 balls, where balls can be either blue or red. Initially, both balls are red. At each stage, one ball is removed at random and replaced by a new ball, which with probability 0.8 is the same color as the removed ball and with probability 0.2 is the opposite color. 

1. Define $X_n$ to be the number of red balls in the urn after the $n$th selection. What are the possible states? What is the transition matrix $\mathbf{P}$ for the resulting Markov chain?

(Type possible states here and enter correct matrix values in R code below.)

```{r}
P <- matrix(0,3,3)
P[1,] <- c(.3,.3,0)
P[2,] <- c(.3,.3,.3)
P[3,] <- c(0,.3,.3)
P
```

2. What is the probability that $X_1=2$, given that $X_0=2$?

3. What is the probability that $X_5=2$, given that $X_0=2$? (Use a power of $\mathbf{P}$ to compute this.)

4. What is the probability in the long run that the chain is in state 2? Solve in two ways: (a) raise  $\mathbf{P}$ to a large power; (b) compute the left eigenvectors of  $\mathbf{P}$ and find the one corresponding to eigenvalue 1.

### Bernoulli-Laplace model of diffusion

Consider two urns $A$ and $B$, each of which contains $m$ balls; $b$ of the $2m$ total balls are black and the remaining $2m-b$ are white. At each stage, one ball is randomly chosen from each urn and put in the other urn (done simultaneously, so balls are being swapped between urns). Define $X_n$ to be the number of black balls in urn $A$ after the $n$th swap. Observe that $X_n=i$ implies $A$ contains $i$ black balls and $m-i$ white balls, while $B$ contains $b-i$ black balls and $m-(b-i)$ white balls. The resulting Markov chain is a probabilistic model of diffusion of two fluids.

Let $m=4$ and $b=4$ for purposes of this lab (to keep calculations manageable). Assume urn $A$ initially contains 4 white balls and $B$ contains 4 black balls.

5. What are the possible states? What is the transition matrix $\mathbf{P}$ for the resulting Markov chain? What kind of boundaries does this chain have?

```{r}
P <- matrix(0,5,5) # 5x5 matrix of zeros 
P[1,2]<-1
P[5,4]<-1 # enter the correct values for P
P[2,1:3]<-c(0/16,0/16,0/16)
P[3,2:4]<-c(0/16,0/16,0/16)
P[4,3:5]<-c(0/16,0/16,0/16)
```

6. What is the probability that $X_2=0$, given that $X_0=0$?

7. What is the probability that $X_3=0$, given that $X_0=0$?

8. What is the probability that $X_6=0$, given that $X_0=0$?

9. Let $\alpha$ be the initial probability distribution corresponding to $X_0=0$, and calculate $\alpha\mathbf{P}^n$, the probability distribution for $X_n$. As $n$ increases, what is $\alpha\mathbf{P}^n$ converging to?

```{r}
alpha<-c(1,0,0,0,0)
n<-10
alpha %*% (P %^% n)
```

10. Find the eigenvalues of $\mathbf{P}$.

```{r}
r <- eigen(t(P))
r$values
```

11. Find the left eigenvector of $\mathbf{P}$ corresponding to eigenvalue 1, and compare to the vector you found in Exercise 9.


