<jupyter_start><jupyter_text>Theory and Practice of Markov Chain Monte Carlo<jupyter_code>%pylab inline<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Motivation MCMC is a computational tool to solve the following problem: suppose you can calculate a function $p(\theta)$ that is proportional to a probability density but do not know its normalization,$$\text{prob}(\theta) = \frac{1}{N} p(\theta) \quad, \quad N = \int p(\theta)\,d\theta \; ,$$how can you calculate integrals like$$\langle f(\theta)\rangle = \frac{1}{N} \int p(\theta) f(\theta)\, d\theta \; .$$The argument $\theta$ is typically a vector of real-valued model parameters. An alternative method for solving this problem is [variational inference](https://arxiv.org/abs/1601.00670).This problem is central to Bayesian inference, where the relevant probability density is the posterior$$P(\theta\,|\,D) = \frac{P(D\,|\,\theta) P(\theta)}{P(D)}$$but we generally only know how to calculate its numerator, the product of the likelihood and prior,$$p(\theta) = P(D\,|\,\theta) P(\theta)$$but not its denominator, the evidence$$N = P(D) = \int P(D\,|\,\theta) P(\theta)\, d\theta \; .$$Integrals like $\langle f(\theta)\rangle$ are required to estimate the mean values and covariances of parameters and to marginalize out nuisance parameters.  Theory of Markov Chains A Markov chain is a sequence of randomly generated parameter values:$$\theta_0\rightarrow \theta_1\rightarrow \theta_2\rightarrow \ldots\rightarrow \theta_{N-1}\rightarrow \theta_N\rightarrow \ldots$$where the probability of obtaining $\theta$ at step $N$ only depends on the previous value $\theta_{N-1}$ and not on any earlier values. In other words, the chain has no memory.In general, the probability also depends on $N$, unless the chain is **time homogeneous** (but most practical MCMC schemes do dependend on $N$).If $\theta$ is limited to a finite set of values, we can think of $\theta$ as a state label and the chain as a sequence of state transitions.  If the chain is also **time homogeneous**, the state transition probabilities can be represented as a matrix, e.g., with three states $A, B, C$:$$T = \begin{pmatrix}p_{AA} & p_{AB} & p_{AC} \\p_{BA} & p_{BB} & p_{BC} \\p_{CA} & p_{CB} & p_{CC}\end{pmatrix} \; .$$Note that probabilities sum to one along each row (but not down each column). If we represent the states with vectors$$A = \begin{pmatrix}1\\ 0\\ 0\end{pmatrix} \quad, \quadB = \begin{pmatrix}0\\ 1\\ 0\end{pmatrix} \quad, \quadC = \begin{pmatrix}0\\ 0\\ 1\end{pmatrix} \; ,$$then each step in the chain has probabilities$$\theta_{N}^t = \theta_{N-1}^t T  \; ,$$and multiple steps are given by matrix multiplication$$\theta_N^t = \theta_0^t T^N \; .$$<jupyter_code>def generate_random_transition_matrix(n):
    T = np.random.uniform(size=(n, n))
    return T / T.sum(axis=1)[:, np.newaxis]
def plot_markov_chain(T, N):
    n = len(T)
    fig, axes = plt.subplots(nrows=1, ncols=N, sharex=True, sharey=False, squeeze=True, figsize=(15., 15. / N))
    TN = T.copy()
    for ax in axes:
        ax.imshow(TN, interpolation='none', origin='upper', extent=(0., 1., 0., 1.), cmap='viridis')
        ax.set_axis_off()
        TN = TN.dot(T)
    plt.subplots_adjust(wspace=0.05)<jupyter_output><empty_output><jupyter_text>Display transition probabilities $T^N$ after $N=1-10$ steps of a finite (time homogeneous) Markov chain. Row $i$ shows the probability of ending in each state $j$ starting from state $i$.<jupyter_code>plot_markov_chain(generate_random_transition_matrix(10), 10)<jupyter_output><empty_output><jupyter_text>Vertical stripes appear because the chain has converged to a stationary distribution that does not depend on the initial state. This is a generic feature of a randomly chosen $T$.  The vector $\theta_s$ of stationary distribution probabilities is a left eigenvector of $T$:$$\theta_s^t = \theta_s^t T \; .$$<jupyter_code>plot_markov_chain(generate_random_transition_matrix(10), 10)
plot_markov_chain(generate_random_transition_matrix(100), 10)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13321926179418236, 'fraction_numerical': 0.013688584698117819, 'mean_word_length': 4.845714285714286, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 1, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '136934', 'license': 'MIT', 'n_tokens_mistral': 1300, 'n_tokens_neox': 1213, 'n_words': 456, 'path': 'Extra/MCMC.ipynb', 'repo_name': 'dkirkby/astroml-study', 'size': '4091'}