<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <a href="https
Step1:  Data
Step2:  Centered model
Step3:  Non-centered
Step4:  Funnel of hell
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import sklearn
import scipy.stats as stats
import scipy.optimize
import matplotlib.pyplot as plt
import seaborn as sns
import time
import numpy as np
import os
import pandas as pd
!pip install -U pymc3>=3.8
import pymc3 as pm
print(pm.__version__)
import theano.tensor as tt
import theano
#!pip install arviz
import arviz as az
!mkdir ../figures
Explanation: <a href="https://colab.research.google.com/github/probml/pyprobml/blob/master/notebooks/schools8_pymc3.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
In this notebook, we fit a hierarchical Bayesian model to the "8 schools" dataset.
See also https://github.com/probml/pyprobml/blob/master/scripts/schools8_pymc3.py
End of explanation
# https://github.com/probml/pyprobml/blob/master/scripts/schools8_pymc3.py
# Data of the Eight Schools Model
J = 8
y = np.array([28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0])
sigma = np.array([15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0])
print(np.mean(y))
print(np.median(y))
names = []
for t in range(8):
    names.append("{}".format(t))
# Plot raw data
fig, ax = plt.subplots()
y_pos = np.arange(8)
ax.errorbar(y, y_pos, xerr=sigma, fmt="o")
ax.set_yticks(y_pos)
ax.set_yticklabels(names)
ax.invert_yaxis()  # labels read top-to-bottom
plt.title("8 schools")
plt.savefig("../figures/schools8_data.png")
plt.show()
Explanation: Data
End of explanation
# Centered model
with pm.Model() as Centered_eight:
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=5)
    sigma_alpha = pm.HalfCauchy("sigma_alpha", beta=5)
    alpha = pm.Normal("alpha", mu=mu_alpha, sigma=sigma_alpha, shape=J)
    obs = pm.Normal("obs", mu=alpha, sigma=sigma, observed=y)
    log_sigma_alpha = pm.Deterministic("log_sigma_alpha", tt.log(sigma_alpha))
np.random.seed(0)
with Centered_eight:
    trace_centered = pm.sample(1000, chains=4, return_inferencedata=False)
pm.summary(trace_centered).round(2)
# PyMC3 gives multiple warnings about  divergences
# Also, see r_hat ~ 1.01, ESS << nchains*1000, especially for sigma_alpha
# We can solve these problems below by using a non-centered parameterization.
# In practice, for this model, the results are very similar.
# Display the total number and percentage of divergent chains
diverging = trace_centered["diverging"]
print("Number of Divergent Chains: {}".format(diverging.nonzero()[0].size))
diverging_pct = diverging.nonzero()[0].size / len(trace_centered) * 100
print("Percentage of Divergent Chains: {:.1f}".format(diverging_pct))
dir(trace_centered)
trace_centered.varnames
with Centered_eight:
    # fig, ax = plt.subplots()
    az.plot_autocorr(trace_centered, var_names=["mu_alpha", "sigma_alpha"], combined=True)
    plt.savefig("schools8_centered_acf_combined.png", dpi=300)
with Centered_eight:
    # fig, ax = plt.subplots()
    az.plot_autocorr(trace_centered, var_names=["mu_alpha", "sigma_alpha"])
    plt.savefig("schools8_centered_acf.png", dpi=300)
with Centered_eight:
    az.plot_forest(trace_centered, var_names="alpha", hdi_prob=0.95, combined=True)
    plt.savefig("schools8_centered_forest_combined.png", dpi=300)
with Centered_eight:
    az.plot_forest(trace_centered, var_names="alpha", hdi_prob=0.95, combined=False)
    plt.savefig("schools8_centered_forest.png", dpi=300)
Explanation: Centered model
End of explanation
# Non-centered parameterization
with pm.Model() as NonCentered_eight:
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=5)
    sigma_alpha = pm.HalfCauchy("sigma_alpha", beta=5)
    alpha_offset = pm.Normal("alpha_offset", mu=0, sigma=1, shape=J)
    alpha = pm.Deterministic("alpha", mu_alpha + sigma_alpha * alpha_offset)
    # alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=J)
    obs = pm.Normal("obs", mu=alpha, sigma=sigma, observed=y)
    log_sigma_alpha = pm.Deterministic("log_sigma_alpha", tt.log(sigma_alpha))
np.random.seed(0)
with NonCentered_eight:
    trace_noncentered = pm.sample(1000, chains=4)
pm.summary(trace_noncentered).round(2)
# Samples look good: r_hat = 1, ESS ~= nchains*1000
with NonCentered_eight:
    az.plot_autocorr(trace_noncentered, var_names=["mu_alpha", "sigma_alpha"], combined=True)
    plt.savefig("schools8_noncentered_acf_combined.png", dpi=300)
with NonCentered_eight:
    az.plot_forest(trace_noncentered, var_names="alpha", combined=True, hdi_prob=0.95)
    plt.savefig("schools8_noncentered_forest_combined.png", dpi=300)
az.plot_forest(
    [trace_centered, trace_noncentered],
    model_names=["centered", "noncentered"],
    var_names="alpha",
    combined=True,
    hdi_prob=0.95,
)
plt.axvline(np.mean(y), color="k", linestyle="--")
az.plot_forest(
    [trace_centered, trace_noncentered],
    model_names=["centered", "noncentered"],
    var_names="alpha",
    kind="ridgeplot",
    combined=True,
    hdi_prob=0.95,
);
Explanation: Non-centered
End of explanation
# Plot the "funnel of hell"
# Based on
# https://github.com/twiecki/WhileMyMCMCGentlySamples/blob/master/content/downloads/notebooks/GLM_hierarchical_non_centered.ipynb
fig, axs = plt.subplots(ncols=2, sharex=True, sharey=True)
x = pd.Series(trace_centered["mu_alpha"], name="mu_alpha")
y = pd.Series(trace_centered["log_sigma_alpha"], name="log_sigma_alpha")
axs[0].plot(x, y, ".")
axs[0].set(title="Centered", xlabel="µ", ylabel="log(sigma)")
# axs[0].axhline(0.01)
x = pd.Series(trace_noncentered["mu_alpha"], name="mu")
y = pd.Series(trace_noncentered["log_sigma_alpha"], name="log_sigma_alpha")
axs[1].plot(x, y, ".")
axs[1].set(title="NonCentered", xlabel="µ", ylabel="log(sigma)")
# axs[1].axhline(0.01)
plt.savefig("schools8_funnel.png", dpi=300)
xlim = axs[0].get_xlim()
ylim = axs[0].get_ylim()
x = pd.Series(trace_centered["mu_alpha"], name="mu")
y = pd.Series(trace_centered["log_sigma_alpha"], name="log sigma_alpha")
sns.jointplot(x, y, xlim=xlim, ylim=ylim)
plt.suptitle("centered")
plt.savefig("schools8_centered_joint.png", dpi=300)
x = pd.Series(trace_noncentered["mu_alpha"], name="mu")
y = pd.Series(trace_noncentered["log_sigma_alpha"], name="log sigma_alpha")
sns.jointplot(x, y, xlim=xlim, ylim=ylim)
plt.suptitle("noncentered")
plt.savefig("schools8_noncentered_joint.png", dpi=300)
group = 0
fig, axs = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=(10, 5))
x = pd.Series(trace_centered["alpha"][:, group], name=f"alpha {group}")
y = pd.Series(trace_centered["log_sigma_alpha"], name="log_sigma_alpha")
axs[0].plot(x, y, ".")
axs[0].set(title="Centered", xlabel=r"$\alpha_0$", ylabel=r"$\log(\sigma_\alpha)$")
x = pd.Series(trace_noncentered["alpha"][:, group], name=f"alpha {group}")
y = pd.Series(trace_noncentered["log_sigma_alpha"], name="log_sigma_alpha")
axs[1].plot(x, y, ".")
axs[1].set(title="NonCentered", xlabel=r"$\alpha_0$", ylabel=r"$\log(\sigma_\alpha)$")
xlim = axs[0].get_xlim()
ylim = axs[0].get_ylim()
plt.savefig("schools8_funnel_group0.png", dpi=300)
Explanation: Funnel of hell
End of explanation
<END_TASK>