<jupyter_start><jupyter_text>Metropolis-Hastings Example 2Adapted from [Codes for Unit 5: metro2.m](https://www2.isye.gatech.edu/isye6420/supporting.html).Associated lecture video: Unit 5 Lesson 9<jupyter_code>rng = np.random.default_rng(1)

obs = 1000000
burn = 500
data = np.array([0.2, 0.1, 0.25])  # observed
data_prod = np.prod(data)  # prod1 in MATLAB code
n = len(data)

# inits
alpha = 2
eta = 2
alphas = np.zeros(obs)
etas = np.zeros(obs)

beta = 2
xi = 2

# pre-generate 3 sets of uniform randoms
unif = rng.uniform(size=obs)
log_unif_alpha = np.log(rng.uniform(size=obs))
log_unif_eta = np.log(rng.uniform(size=obs))

for i in track(range(obs)):
    # proposals
    alpha_prop = -alpha * log_unif_alpha[i]
    eta_prop = -eta * log_unif_eta[i]

    # prod2 in MATLAB code
    prop_prod = np.prod(np.exp(eta * data**alpha - eta_prop * data**alpha_prop))

    # acceptance ratio
    rr = (
        (eta_prop / eta) ** (beta - 1)
        * np.exp(alpha - alpha_prop - xi * (eta_prop - eta))
        * np.exp(
            -alpha / alpha_prop - eta / eta_prop + alpha_prop / alpha + eta_prop / eta
        )
        * data_prod ** (alpha_prop - alpha)
        * prop_prod
        * ((alpha_prop * eta_prop) / (alpha * eta)) ** (n - 1)
    )

    rho = min(rr, 1)
    if unif[i] < rho:
        alpha = alpha_prop
        eta = eta_prop

    alphas[i] = alpha
    etas[i] = eta

alphas = alphas[burn:]
etas = etas[burn:]

print(f"{np.mean(alphas)=}")
print(f"{np.var(alphas)=}")
print(f"{np.mean(etas)=}")
print(f"{np.var(etas)=}")

plt.hist(alphas, 50)
plt.xlabel("alpha")
plt.show()

plt.hist(etas, 50)
plt.xlabel("eta")
plt.show()
%watermark --iversions -v<jupyter_output>Python implementation: CPython
Python version       : 3.10.4
IPython version      : 8.4.0

matplotlib: 3.5.2
numpy     : 1.22.4
{'fraction_non_alphanumeric': 0.15559246954595793, 'fraction_numerical': 0.031007751937984496, 'mean_word_length': 2.962719298245614, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '434821', 'license': 'Apache-2.0', 'n_tokens_mistral': 797, 'n_tokens_neox': 743, 'n_words': 154, 'path': 'unit5/Unit5-MetropolisHastings2.ipynb', 'repo_name': 'areding/6420-pymc', 'size': '1806'}