# COMP 136 Day 21 - Discrete Markov Chains and Stationary Distributions


## Outline

* [Part 1: Discrete Markov Models](#part1)
* [Part 2: Limiting behavior of the markov model as t -> infinity](#part2)
* [Part 3: Limiting behavior of many multiplies of transition matrix A](#part3)
* [Part 4: Eigenvectors of transition matrix A](#part4)

## Key takeaways

* New concept: 'Stationary distribution' or 'Equilibrium distribution', the limiting distribution of marginal $p(z_t)$ as $t \rightarrow \infty$

* When does the stationary distribution exist? As long it is "ergodic", which (intuitively) means each state has *some* greater-than-zero probability of reaching every other state after some number of steps.

* Below, we'll see 3 ways to compute a stationary distribution:
* * Manually compute the marginal $p(z_t)$ at each timestep $t \in 1, 2, \ldots T$. Observe it become stationary
* * Look at limits of multiplying the transition matrix: $A \cdot A \cdot A \cdot \ldots$. Eventually, will converge to a matrix where rows are the stationary distribution.
* * Look at eigenvectors of $A^T$ (A transposed). Find the eigenvector corresponding to eigenvalue 1 and renormalize it so it sums to one.



# Part 1: Discrete Markov model

#### Random variables

* $z_1, z_2, \ldots z_T$, where $t \in \{1, 2, \ldots T$ indexes the discrete timestep
* Each is a discrete random variable: $z_t \in \{1, 2, \ldots K\}$

#### Parameters

* $\pi$ is the $K$-length initial probability vector. $\pi \in \Delta^K$ (non-negative, sums to one)
* $A$ is a $K \times K$ transition matrix. For each row $j$, we have $A_j \in \Delta^K$ (all entries are non-negative and each *row* sums to one).

#### Probability mass function

* $p(z_1) = \text{CatPMF}(z_1 | \pi_1, \pi_2, \ldots \pi_K)$, which implies $p(z_1 = k) = \pi_k$
* $p(z_t | z_{t-1}=j) = \text{CatPMF}(z_t | A_{j1}, A_{j2}, \ldots A_{jK})$, which implies $p(z_t = k | z_{t-1}=j) = A_{jk}$


#### Things to remember

In the code below, we need to use zero-based indexing (like python always does). 

So, the "first" timestep is t=0.



```python
import numpy as np
```


```python
## Number of states
K = 2
```


```python
## Initial state probabilities
pi_K = np.asarray([0.5, 0.5])
```


```python
pi_K
```




    array([0.5, 0.5])




```python
## Transition probabilities
A_KK = np.asarray([[0.9, 0.1], [0.2, 0.8]])
```


```python
print(A_KK)
```

    [[0.9 0.1]
     [0.2 0.8]]


## Math Discussion 1a: Understanding marginal probabilities

Assume that for some given value of index $t$ and for all states $\ell \in 1, 2, \ldots K$, the following is known:

$$
p( z_{t-1} = \ell ) = \rho_{t-1,\ell}
$$

where $\rho_{t-1}$ is a given vector of size $K$ that sums to one.

What principle of probability can we use to show the following is true:

\begin{align}
p( z_t = k ) &= \sum_{\ell = 1}^K p(z_{t-1} = \ell) p( z_t = k | z_{t-1} = \ell)
\\
&=  \sum_{\ell = 1}^K \rho_\ell A_{\ell k}
\end{align}


```python
# TODO discuss
```

## Math to code exercise 1b: Understanding matrix math

Given a vector previous_rho_K and a matrix A_KK, write a compact one-line numpy expression that will compute the vector $\rho_t = [\rho_{t1}, \ldots \rho_{tK}]$ for the current timestep, such that:

$$
\rho_{t,k} = p(z_t = k)
$$


```python
current_rho_K = None # TODO write in terms of previous_rho_K and A_KK
# Hint: use a matrix multiply via np.dot(..., ....)
```

# Part 2: What happens to the marginal probability of z in the limit of many timesteps?

Let's look at what happens to the distribution over $z_t$ after dozens and dozens of timesteps....

We'll start with a uniform distribution over the initial state: $p(z_1 = 1) = p(z_1 = 2) = 0.5$


```python
pi_K = np.asarray([0.5, 0.5])
for t in range(100):
    if t == 0:
        proba_K = 1.0 * pi_K
    elif t > 0:
        proba_K = np.dot(proba_K, A_KK) # Matrix multiply computation for marginal of timestep t given previous timestep
    print("after t=%3d steps: p(z[t]) = Cat(%.4f, %.4f)" % (t, proba_K[0], proba_K[1]))
```

    after t=  0 steps: p(z[t]) = Cat(0.5000, 0.5000)
    after t=  1 steps: p(z[t]) = Cat(0.5500, 0.4500)
    after t=  2 steps: p(z[t]) = Cat(0.5850, 0.4150)
    after t=  3 steps: p(z[t]) = Cat(0.6095, 0.3905)
    after t=  4 steps: p(z[t]) = Cat(0.6267, 0.3734)
    after t=  5 steps: p(z[t]) = Cat(0.6387, 0.3613)
    after t=  6 steps: p(z[t]) = Cat(0.6471, 0.3529)
    after t=  7 steps: p(z[t]) = Cat(0.6529, 0.3471)
    after t=  8 steps: p(z[t]) = Cat(0.6571, 0.3429)
    after t=  9 steps: p(z[t]) = Cat(0.6599, 0.3401)
    after t= 10 steps: p(z[t]) = Cat(0.6620, 0.3380)
    after t= 11 steps: p(z[t]) = Cat(0.6634, 0.3366)
    after t= 12 steps: p(z[t]) = Cat(0.6644, 0.3356)
    after t= 13 steps: p(z[t]) = Cat(0.6651, 0.3349)
    after t= 14 steps: p(z[t]) = Cat(0.6655, 0.3345)
    after t= 15 steps: p(z[t]) = Cat(0.6659, 0.3341)
    after t= 16 steps: p(z[t]) = Cat(0.6661, 0.3339)
    after t= 17 steps: p(z[t]) = Cat(0.6663, 0.3337)
    after t= 18 steps: p(z[t]) = Cat(0.6664, 0.3336)
    after t= 19 steps: p(z[t]) = Cat(0.6665, 0.3335)
    after t= 20 steps: p(z[t]) = Cat(0.6665, 0.3335)
    after t= 21 steps: p(z[t]) = Cat(0.6666, 0.3334)
    after t= 22 steps: p(z[t]) = Cat(0.6666, 0.3334)
    after t= 23 steps: p(z[t]) = Cat(0.6666, 0.3334)
    after t= 24 steps: p(z[t]) = Cat(0.6666, 0.3334)
    after t= 25 steps: p(z[t]) = Cat(0.6666, 0.3334)
    after t= 26 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 27 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 28 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 29 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 30 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 31 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 32 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 33 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 34 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 35 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 36 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 37 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 38 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 39 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 40 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 41 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 42 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 43 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 44 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 45 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 46 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 47 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 48 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 49 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 50 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 51 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 52 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 53 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 54 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 55 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 56 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 57 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 58 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 59 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 60 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 61 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 62 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 63 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 64 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 65 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 66 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 67 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 68 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 69 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 70 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 71 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 72 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 73 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 74 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 75 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 76 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 77 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 78 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 79 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 80 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 81 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 82 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 83 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 84 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 85 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 86 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 87 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 88 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 89 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 90 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 91 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 92 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 93 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 94 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 95 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 96 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 97 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 98 steps: p(z[t]) = Cat(0.6667, 0.3333)
    after t= 99 steps: p(z[t]) = Cat(0.6667, 0.3333)


### Exercise 2b: Try starting instead from another initial condition: pi_K = [0.01, 0.99]? What happens?


```python
pi_K = np.asarray([0.01, 0.99])
for t in range(100):
    if t == 0:
        proba_K = 1.0 * pi_K
    elif t > 0:
        proba_K = np.zeros(K) # <<-- TODO determine the probability p(z_t = k); Hint: use a matrix product.
    print("after t=%3d steps: p(z[t]) = Cat(%.4f, %.4f)" % (t, proba_K[0], proba_K[1]))
```

    after t=  0 steps: p(z[t]) = Cat(0.0100, 0.9900)
    after t=  1 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  2 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  3 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  4 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  5 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  6 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  7 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  8 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  9 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 10 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 11 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 12 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 13 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 14 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 15 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 16 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 17 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 18 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 19 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 20 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 21 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 22 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 23 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 24 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 25 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 26 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 27 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 28 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 29 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 30 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 31 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 32 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 33 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 34 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 35 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 36 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 37 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 38 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 39 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 40 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 41 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 42 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 43 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 44 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 45 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 46 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 47 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 48 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 49 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 50 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 51 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 52 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 53 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 54 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 55 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 56 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 57 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 58 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 59 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 60 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 61 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 62 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 63 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 64 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 65 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 66 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 67 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 68 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 69 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 70 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 71 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 72 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 73 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 74 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 75 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 76 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 77 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 78 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 79 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 80 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 81 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 82 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 83 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 84 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 85 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 86 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 87 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 88 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 89 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 90 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 91 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 92 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 93 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 94 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 95 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 96 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 97 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 98 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 99 steps: p(z[t]) = Cat(0.0000, 0.0000)


### Exercise 2c: Try starting instead from pi_K = [0.99, 0.01]? What happens?


```python
for t in range(100):
    if t == 0:
        proba_K = np.zeros(K) # <<-- TODO set the initial distribution (see above)
    elif t > 0:
        proba_K = np.zeros(K) # <<-- TODO determine the probability p(z_t = k); Hint: use a matrix product.
    print("after t=%3d steps: p(z[t]) = Cat(%.4f, %.4f)" % (t, proba_K[0], proba_K[1]))
```

    after t=  0 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  1 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  2 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  3 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  4 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  5 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  6 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  7 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  8 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t=  9 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 10 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 11 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 12 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 13 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 14 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 15 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 16 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 17 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 18 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 19 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 20 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 21 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 22 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 23 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 24 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 25 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 26 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 27 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 28 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 29 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 30 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 31 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 32 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 33 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 34 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 35 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 36 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 37 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 38 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 39 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 40 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 41 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 42 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 43 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 44 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 45 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 46 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 47 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 48 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 49 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 50 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 51 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 52 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 53 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 54 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 55 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 56 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 57 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 58 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 59 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 60 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 61 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 62 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 63 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 64 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 65 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 66 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 67 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 68 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 69 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 70 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 71 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 72 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 73 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 74 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 75 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 76 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 77 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 78 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 79 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 80 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 81 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 82 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 83 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 84 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 85 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 86 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 87 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 88 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 89 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 90 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 91 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 92 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 93 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 94 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 95 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 96 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 97 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 98 steps: p(z[t]) = Cat(0.0000, 0.0000)
    after t= 99 steps: p(z[t]) = Cat(0.0000, 0.0000)


## Discussion 2d: The parts 2a-2c above all use the same transition matrix A.... do they all converge to the same limiting distribution?


```python
# TODO discuss
```

# Part 3: What is the limit of many products of the transition matrix?

Run the code below to see what happens if we define the following sequence of matrix multiplies:

\begin{align}
A_1 &= A
\\
A_2 &= A A
\\
\vdots
\\
A_t &= \underbrace{A A \ldots A}_{t ~\text{matrix multiplies}}
\end{align}


```python
np.set_printoptions(precision=2)
val_KK = A_KK
for t in range(100):
    val_KK = np.dot(val_KK, A_KK)
    msg = "--- after %3d steps: A_KK=\n%s" % (t, str(val_KK))
    print(msg)
```

    --- after   0 steps: A_KK=
    [[0.83 0.17]
     [0.34 0.66]]
    --- after   1 steps: A_KK=
    [[0.78 0.22]
     [0.44 0.56]]
    --- after   2 steps: A_KK=
    [[0.75 0.25]
     [0.51 0.49]]
    --- after   3 steps: A_KK=
    [[0.72 0.28]
     [0.55 0.45]]
    --- after   4 steps: A_KK=
    [[0.71 0.29]
     [0.59 0.41]]
    --- after   5 steps: A_KK=
    [[0.69 0.31]
     [0.61 0.39]]
    --- after   6 steps: A_KK=
    [[0.69 0.31]
     [0.63 0.37]]
    --- after   7 steps: A_KK=
    [[0.68 0.32]
     [0.64 0.36]]
    --- after   8 steps: A_KK=
    [[0.68 0.32]
     [0.65 0.35]]
    --- after   9 steps: A_KK=
    [[0.67 0.33]
     [0.65 0.35]]
    --- after  10 steps: A_KK=
    [[0.67 0.33]
     [0.66 0.34]]
    --- after  11 steps: A_KK=
    [[0.67 0.33]
     [0.66 0.34]]
    --- after  12 steps: A_KK=
    [[0.67 0.33]
     [0.66 0.34]]
    --- after  13 steps: A_KK=
    [[0.67 0.33]
     [0.66 0.34]]
    --- after  14 steps: A_KK=
    [[0.67 0.33]
     [0.66 0.34]]
    --- after  15 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  16 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  17 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  18 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  19 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  20 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  21 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  22 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  23 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  24 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  25 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  26 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  27 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  28 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  29 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  30 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  31 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  32 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  33 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  34 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  35 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  36 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  37 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  38 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  39 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  40 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  41 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  42 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  43 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  44 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  45 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  46 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  47 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  48 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  49 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  50 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  51 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  52 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  53 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  54 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  55 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  56 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  57 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  58 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  59 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  60 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  61 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  62 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  63 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  64 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  65 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  66 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  67 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  68 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  69 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  70 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  71 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  72 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  73 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  74 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  75 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  76 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  77 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  78 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  79 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  80 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  81 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  82 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  83 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  84 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  85 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  86 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  87 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  88 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  89 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  90 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  91 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  92 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  93 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  94 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  95 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  96 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  97 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  98 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]
    --- after  99 steps: A_KK=
    [[0.67 0.33]
     [0.67 0.33]]


## Discussion 3a: What do you notice about raising A to a power (many multiplies)? What is the connection to the limiting distribution you found in Part 2?


```python
# TODO discuss
```

# Part 4: Look at the eigenvalues of the transition matrix

Let's look at the eigenvalues of the (transposed) transition matrix.

Remember that any eigenvector is a a K-dimensional vector $v$ that satisfies:

$$
A^T v = \lambda v, \qquad v \in \mathbb{R}^K, \lambda \in \mathbb{R}
$$

We use transpose so that instead of left-multiplication $v A$ (as above), we can do right-multiplication: $A^T v$, and get the same answer (transposed).


```python
lam_K, V_KK = np.linalg.eig(A_KK.T)
```

Print out the raw eigenvalues and eigenvectors 


```python
lam_K
```




    array([1. , 0.7])




```python
V_KK
```




    array([[ 0.89, -0.71],
           [ 0.45,  0.71]])



Remember that by default, eigenvectors are vectors on the unit circle (L2 norm is 1).

But any scalar multiple of an eigenvector is also an eigenvector, for any scalar $s > 0$:

$$
A (s v) = s (A v) = s ( \lambda v)
$$

So let's find the eigenvector that sums to one (and thus is a valid Categorical distribution over K=2 states


```python
proba_K = V_KK[:,0] / np.sum(V_KK[:,0])
```


```python
proba_K
```




    array([0.67, 0.33])



Look, our stationary distribution is a (left)-eigenvector of the transition matrix! With eigenvalue 1!


```python
np.dot(proba_K, A_KK)
```




    array([0.67, 0.33])




```python

```


```python

```

8ff1b73f492521549a5ce73895e216e4ac34d33f
ipynb
Jupyter Notebook
notebooks/day21-MarkovChainsAndStationaryDistributions.ipynb
ErliCai/comp136-21s-assignments
da690dc23a09624a16ba96e1b34ef991f4001a6e
MIT
2021-02-17T08:09:17.000Z
2021-02-17T08:09:17.000Z
notebooks/day21-MarkovChainsAndStationaryDistributions.ipynb
ErliCai/comp136-21s-assignments
da690dc23a09624a16ba96e1b34ef991f4001a6e
MIT
2021-02-01T21:19:13.000Z
2021-04-17T23:25:22.000Z
notebooks/day21-MarkovChainsAndStationaryDistributions.ipynb
ErliCai/comp136-21s-assignments
da690dc23a09624a16ba96e1b34ef991f4001a6e
MIT
2021-02-08T21:54:43.000Z
2022-02-15T22:11:23.000Z
Qwen/Qwen-72B
1. YES
2. YES