<jupyter_start><jupyter_text>Monte Carlo Markov Chain  Sample from a bivariate normal distribution<jupyter_code>import numpy as np 
import matplotlib.pyplot as plt 
plt.rcParams.update({'font.size': 16}) 
import scipy.stats as stats
# parameters
muX = 78.8
sigX = 3.668
muY = 211
sigY = 26.904
rho = 0.81
n=10000
def logp(x):
    if x < 0:
        return 0
    else :
        return np.log(x)
# Gibbs sampler 
XgY = np.empty(n)
YgX = np.empty(n)
# initialize 
XgY[0] = muX

for i in range(1,n):
    # update 
    x = XgY[i-1]
    # sample 
    YgX[i] = np.random.normal(muY + (sigY/sigX) * rho * (x - muX), np.sqrt((1 - rho**2)*sigY**2))
    # update 
    y = YgX[i]
    # sample
    XgY[i] = np.random.normal(muX + (sigX/sigY) * rho * (y - muY), np.sqrt((1 - rho**2)*sigX**2))
acc_rate = (len(np.unique(XgY)) / len(XgY)) * 100
print(f'Acceptance rate = {acc_rate:.0f}%')
# Metropolis-Hastings 
XgYMH = np.empty(n)
YgXMH = np.empty(n)
# initialize 
XgYMH[0] = muX
YgXMH[0] = muY
# iterate 
for i in range(1,n): 
    # 1. Generate random candidate
    xgy = np.random.normal(XgYMH[i-1], 5)
    ygx = np.random.normal(YgXMH[i-1], 30)

    # 2. Calculate acceptance probability 
    num = stats.multivariate_normal.pdf(np.array([xgy, ygx]), np.array([muX, muY]), np.array([[sigX**2, rho*sigX*sigY],[rho*sigX*sigY, sigY**2]]))
    denom = stats.multivariate_normal.pdf(np.array([XgYMH[i-1], YgXMH[i-1]]), np.array([muX, muY]), np.array([[sigX**2, rho*sigX*sigY],[rho*sigX*sigY, sigY**2]]))
    A = np.min([np.log(1), logp(num/denom)])

    # 3. Accept or reject 
    u = np.log(np.random.uniform(low=0, high=1))
    if u <= A: # accept 
        YgXMH[i] = ygx
        XgYMH[i] = xgy    
    else : # reject 
        YgXMH[i] = YgXMH[i-1]
        XgYMH[i] = XgYMH[i-1]
acc_rate = (len(np.unique(YgXMH)) / len(YgXMH)) * 100
print(f'Acceptance rate = {acc_rate:.0f}%')<jupyter_output>Acceptance rate = 33%<jupyter_text>Marginal distributions<jupyter_code>plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.hist(XgYMH[1000:]) # 1000 burn in 
plt.xlabel('x')
plt.ylabel('Frequency')
plt.subplot(1,2,2)
plt.hist(YgXMH[1000:]) #
plt.xlabel('y')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Joint distribution<jupyter_code>plt.figure()
plt.scatter(XgYMH[1000:], YgXMH[1000:])
plt.xlabel('X')
plt.ylabel('Y')
plt.show()<jupyter_output><empty_output><jupyter_text>Numpy sampler<jupyter_code>test = np.random.multivariate_normal(np.array([muX, muY]), np.array([[sigX**2, rho*sigX*sigY],[rho*sigX*sigY, sigY**2]]), size=(n,))<jupyter_output><empty_output><jupyter_text>Numpy v Gibbs v Metropolis-Hastings<jupyter_code># numpy v gibbs v mh 
plt.figure(figsize=(9,7))
plt.scatter(test[:,0], test[:,1], facecolors='none', edgecolor='r', label='numpy')
plt.scatter(XgY[1000:], YgX[1000:], facecolors='none', edgecolor='b', label='MCMC-Gibbs')
plt.scatter(XgYMH[1000:], YgXMH[1000:], facecolors='none', edgecolor='g', label='MCMC-MH')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc=0)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19492084432717677, 'fraction_numerical': 0.0395778364116095, 'mean_word_length': 4.606284658040665, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 20, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '454034', 'license': 'MIT', 'n_tokens_mistral': 1502, 'n_tokens_neox': 1352, 'n_words': 222, 'path': 'Bivariate_normal.ipynb', 'repo_name': 'tulimid1/MCMC', 'size': '3032'}