<jupyter_start><jupyter_text>Gibbs sampling in 2DThis is BONUS content related to Day 22, where we introduce Gibbs sampling Random variables(We'll use 0-indexing so we have close alignment between math and python code)* 2D random variable $z = [z_0, z_1]$* each entry $z_d$ is a real scalar: $z_d \in \mathbb{R}$ Target distribution\begin{align}p^*(z_0, z_1) =  \mathcal{N}\left(    \left[ \begin{array}{c}    0 \\ 0    \end{array} \right],    \left[    \begin{array}{c c}    1 & 0.8 \\    0.8 & 2    \end{array} \right] \right)\end{align} Key takeaways* New concept: 'Gibbs sampling', which just iterates between two conditional sampling distributions:\begin{align}    z^{t+1}_0 &\sim p^* (z_0 | z_1 = z^t_1) \\    z^{t+1}_1 &\sim p^* (z_1 | z_0 = z^{t+1}_0)\end{align} Things to rememberThis is a simple example to illustrate the idea of how Gibbs sampling works.There are other "better" ways of sampling from a 2d normal.  Setup<jupyter_code>import numpy as np

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=2.0)<jupyter_output><empty_output><jupyter_text>Step 1: Prepare for Gibbs sampling Define functions to sample from target's conditionals<jupyter_code>def draw_z0_given_z1(z1, random_state):
    ## First, use Bishop textbook formulas to compute the conditional mean/var
    mean_01 = 0.4 * z1
    var_01 = 0.68
    
    ## Then, use simple transform to obtain a sample from this conditional
    ## Remember, if u ~ Normal(0, 1), a "standard" normal with mean 0 variance 1,
    ## then using transform: x <- T(u), with T(u) = \mu + \sigma * u
    ## we can say x ~ Normal(\mu, \sigma^2)
    u_samp = random_state.randn()
    z0_samp = mean_01 + np.sqrt(var_01) * u_samp
    return z0_samp
def draw_z1_given_z0(z0, random_state):
    ## First, use Bishop textbook formulas to compute conditional mean/var
    mean_10 = 0.8 * z0
    var_10 = 1.36
    
    ## Then, use simple transform to obtain a sample from this conditional
    ## Remember, if u ~ Normal(0, 1), a "standard" normal with mean 0 variance 1,
    ## then using transform: x <- T(u), with T(u) = \mu + \sigma * u
    ## we can say x ~ Normal(\mu, \sigma^2)
    u_samp = random_state.randn()
    z1_samp = mean_10 + np.sqrt(var_10) * u_samp
    return z1_samp<jupyter_output><empty_output><jupyter_text>Step 2: Execute the Gibbs sampling algorithmPerform 6000 iterations.Discard the first 1000 as "not yet burned in".<jupyter_code>S = 6000
sample_list = list()
z_D = np.zeros(2)

random_state = np.random.RandomState(0) # reproducible random seeds

for t in range(S):
    z_D[0] = draw_z0_given_z1(z_D[1], random_state)
    z_D[1] = draw_z1_given_z0(z_D[0], random_state)
    
    if t > 1000:
        sample_list.append(z_D.copy()) # save copies so we get different vectors
z_samples_SD = np.vstack(sample_list)<jupyter_output><empty_output><jupyter_text>Step 3: Compare to samples from built-in routines for 2D MVNormal sampling<jupyter_code>Cov_22 = np.asarray([[1.0, 0.8], [0.8, 2.0]])
true_samples_SD = random_state.multivariate_normal(np.zeros(2), Cov_22, size=S-1000)
fig, ax_grid = plt.subplots(nrows=1, ncols=2, sharex=True, sharey=True, figsize=(10,4))

ax_grid[0].plot(z_samples_SD[:,0], z_samples_SD[:,1], 'k.')
ax_grid[0].set_title('Gibbs sampler')
ax_grid[0].set_aspect('equal', 'box');

ax_grid[1].plot(true_samples_SD[:,0], true_samples_SD[:,1], 'k.')
ax_grid[1].set_title('np.random.multivariate_normal')
ax_grid[1].set_aspect('equal', 'box');
ax_grid[1].set_xlim([-6, 6]);
ax_grid[1].set_ylim([-6, 6]);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1349095966620306, 'fraction_numerical': 0.04089012517385257, 'mean_word_length': 4.407518796992481, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 18, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '487432', 'license': 'MIT', 'n_tokens_mistral': 1419, 'n_tokens_neox': 1326, 'n_words': 367, 'path': 'notebooks/GibbsSampling.ipynb', 'repo_name': 'tufts-ml-courses/comp136-spr-20s-assignments-', 'size': '3595'}