<jupyter_start><jupyter_text>Markov ChainsA Markov chain is a sequence of events in which the probability of the next event depends only on the state of the current event. For example, we have previously encountered Markov chains in the random walk and Google Page Rank algorithm.   Example: Random walk and diffusion<jupyter_code>n = 1000
k = 100
xs = np.c_[np.zeros(n), np.random.choice([-1,1], (n, k))]
xs = np.cumsum(xs, axis=1)
for x in xs:
    plt.plot(np.arange(k+1), x)<jupyter_output><empty_output><jupyter_text>Snapshots of displacement after $k$ steps follow normal distribution with $\sigma = \sqrt{k}$<jupyter_code>fig, axes = plt.subplots(1, 5, figsize=(15, 3), sharex=True)
for i, k in enumerate(range(20, 101, 20)):
    axes[i].hist(xs[:, k], alpha=0.3, normed=True)
    axes[i].set_title('t = %d' % k)
    xp = np.linspace(-25, 25, 100)
    axes[i].plot(xp, stats.norm(loc=0, scale=np.sqrt(k)).pdf(xp))
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Properties of states and Markov chains A Markov chain is irreducible if it is possible to get from any state to any state. Otherwise it is reducible. A state has period $k$ if it must return to that state in multiples of $k$ moves. If $k=1$, the state is aperiodic. If all states are aperiodic, then the Markov chain is aperiodic. If any state in an irreducible Markov chain is aperiodic, then all states are aperiodic.   Reducible and aperiodicA is a transient state. B and C are recurrent states.<jupyter_code>%%dot

digraph g {
    node [shape=circle]
    A -> B
    B -> C
    C -> B
    C -> C
}<jupyter_output><empty_output><jupyter_text>Irreducible and periodicA, B, and C are recurrent states.<jupyter_code>%%dot

digraph g {
    layout="circo";
    node [shape=circle]
    A -> B
    B -> C
    C -> A
}<jupyter_output><empty_output><jupyter_text>Absorbing Markov chainD is an absorbing state. A, B and C are transient states. Since every state can reach an absorbing state, this is an absorbing Markov chain<jupyter_code>%%dot

digraph g {
    node [shape=circle]
    A -> B
    B -> C
    C -> B
    C -> C
    B -> D
}<jupyter_output><empty_output><jupyter_text>ErgodicAn ergodic state is aperiodic and positive recurrent. If all states are ergodic, then we have an ergodic Markov chain. A finite chain that is aperiodic and irreducible is ergodic.<jupyter_code>%%dot

digraph g {
    layout="circo";
    node [shape=circle]
    A -> B
    B -> C
    C -> A
    C -> C
}<jupyter_output><empty_output><jupyter_text>Markov ChainSuppose we have the following transition table:|             | Raleigh | Chapel Hill | Durham ||:------------|:--------|:------------|:-------|| Raleigh     | 0.9     | 0.05        | 0.05   || Chapel Hill | 0.1     | 0.8         | 0.1    || Durham      | 0.04    | 0.01        | 0.95   |This says that a resident of Raleigh this year has a 90% chance of stying in Raleigh and a 5% chance of relocating to Chapel Hill or Durham the following year. Note that the probabilities only depend on the current year. The transitions between states (Raleigh, Chapel Hill and Durham) are then said to be modeled by a Markov chain. If Raleigh, Chapel Hill and Durham each started with 300,000 residents, what is the long run distribution assuming no immigration or emigration?<jupyter_code>%%dot

digraph g {
    node [shape=box]
    R [label="Raleigh"]
    C [label="Chapel Hill"]
    D [label="Duham"]
    R -> R [label=0.9]
    R -> C [label=0.05]
    R -> D [label=0.05]
    C -> R [label=0.1]
    C -> C [label=0.8]
    C -> D [label=0.1]
    D -> R [label=0.04]
    D -> C [label=0.01]
    D -> D [label=0.95]
}
A = np.array([
    [0.9, 0.05, 0.05],
    [0.1, 0.8, 0.1],
    [0.04, 0.01, 0.95]
])
x = np.array([300000, 300000, 300000]).reshape(-1,1)<jupyter_output><empty_output><jupyter_text>After one year<jupyter_code>x.T @ A<jupyter_output><empty_output><jupyter_text>After 2 years<jupyter_code>x.T @ A @ A<jupyter_output><empty_output><jupyter_text>Brute force solution<jupyter_code>n = 100
x.T @ np.linalg.matrix_power(A, n)<jupyter_output><empty_output><jupyter_text>Eigenvector solutionAt steady state, we have$$p.T A = p.T$$Taking the transpose on both sides, we see that$$A^Tp = p$$suggesting that the steady state probability vector is the normalized eigenvector of $A^T$ with $\lambda = 1$.<jupyter_code>E, V = np.linalg.eig(A.T)
E
p = V[:, 2:] / V[:, 2:].sum()
p
p * x.sum()<jupyter_output><empty_output><jupyter_text>Linear system solutionWe cannot solve $A^Tp - p= \mathbb{0}$ directly because $A^T - I$ does not have full rank (since the last column of $A$ is fixed given the other columns).We need to add the constraint $\sum{p} = 1$ and then we can solve the linear system to find $p$. $$\begin{bmatrix}(A^T - I) \\\mathbb{1}^T \end{bmatrix} p = \begin{bmatrix}\mathbb{0} \\1 \end{bmatrix}$$<jupyter_code>np.linalg.matrix_rank(A.T - np.eye(3))
M = np.r_[A.T - np.eye(3), [np.ones(3)]]
b = np.r_[np.zeros(3), 1].reshape(-1,1)
M
b
x, res, rank, s = np.linalg.lstsq(M, b, rcond=None)
x<jupyter_output><empty_output><jupyter_text>Markov chain averagesWe can also consider the perspective of a single individual in terms of the frequencies of places visited.The individual starts from one of the 3 places (Raleigh, Chapel Hill or Durham) and moves from place to place according to the probabilities in $A$ over a long time. In the long run, the average frequency of visits to a place is the steady state probability of visiting that place.This suggests that we can estimate the equilibrium or steady state distribution by using the long-run averages, which we'll revisit when we look at Markov chain Monte Carlo (MCMC). Note that the initial averages depend on the starting conditions, so it is common practice to discard some number of the initial iterations (burn-in).<jupyter_code>P = np.cumsum(A).reshape(A.shape) - np.array([0,1,2]).reshape(-1,1)
P
n = int(1e6)
xs = np.zeros(n+1).astype('int')
xs[0] = np.random.choice([0,1,2])

for t in range(n):
    r = np.random.rand()
    xs[t+1] = np.argmax(r < P[xs[t]])
np.mean(xs==0), np.mean(xs==1), np.mean(xs==2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13137254901960785, 'fraction_numerical': 0.027450980392156862, 'mean_word_length': 3.7671339563862927, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 78, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '167766', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2286, 'n_tokens_neox': 2096, 'n_words': 739, 'path': 'notebooks/S10B_MarkovChains.ipynb', 'repo_name': 'taotangtt/sta-663-2018', 'size': '6120'}