Robust Modeling

!sudo /bin/bash -c "(source /venv/bin/activate; pip install --quiet jupyterlab-vim)"
!jupyter labextension enable

!sudo /bin/bash -c "(source /venv/bin/activate; pip install --quiet graphviz)"

!sudo /bin/bash -c "(source /venv/bin/activate; pip install --quiet dataframe_image)"

%load_ext autoreload
%autoreload 2


import arviz as az
import pymc as pm
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

dir_name = "./L07_data"

!ls $dir_name

import msml610.tutorials.msml610_utils as ut

ut.config_notebook()

data = np.loadtxt(f"{dir_name}/chemical_shifts.csv")
print(len(data), data)
# print(sorted(data))

# It looks Gaussian with a couple of outliers.
az.plot_kde(data, rug=True)

title = "Chemical shift"
ut.process_figure(title)

with pm.Model() as model_g:
    # The mean is Uniform in [40, 70] (which is larger than the data).
    mu = pm.Uniform("mu", lower=40, upper=70)
    # The std dev is half normal with a large value (which is a large value based on the data).
    sigma = pm.HalfNormal("sigma", sigma=10)
    # The model is N(mu, sigma).
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=data)
    # Sample.
    idata_g = pm.sample(1000)

pm.model_to_graphviz(model_g)

# There are 4 traces for 2 variables.
az.plot_trace(idata_g)

# The posterior distribution of the params is bi-dimensional, since it has mu and sigma.
az.plot_pair(idata_g, kind="kde", marginals=True)

# Report a summary of the inference.
az.summary(idata_g, kind="stats").round(2)

# Compute 100 posterior predictive samples.
y_pred_g = pm.sample_posterior_predictive(idata_g, model=model_g)

# Black: KDE of the data (observed)
# Blue: KDEs of the posterior predictive samples
# Orange: KDE of the posterior predictive mean
az.plot_ppc(y_pred_g, mean=True, num_pp_samples=100)

# Show the PDF for various values of \nu.

# Points to be used to sample the PDF.
x_values = np.linspace(-10, 10, 500)

# Plot t-student sweeping \nu.
for df in [0.1, 0.5, 1, 2, 5, 10, 30]:
    # Student-t with df.
    distr = stats.t(df)
    # Compute PDF.
    x_pdf = distr.pdf(x_values)
    plt.plot(x_values, x_pdf, label=f"nu={df}")

# Plot gaussian.
x_pdf = stats.norm.pdf(x_values)
plt.plot(x_values, x_pdf, "k--", label="Gauss / nu=infty")

plt.xlim(-5, 5)
plt.legend()

title = "Chap7: Student-t"
ut.process_figure(title)

# Use a Student-T model.
with pm.Model() as model_t:
    mu = pm.Uniform("mu", 40, 75)
    sigma = pm.HalfNormal("sigma", sigma=10)
    # A student with nu = 30 is close to a Gaussian.
    nu = pm.Exponential("nu", 1 / 30)
    #
    y = pm.StudentT("y", mu=mu, sigma=sigma, nu=nu, observed=data)
    idata_t = pm.sample(1_000)

az.plot_trace(idata_t)

az.summary(idata_t, kind="stats").round(2)

# Compute 100 posterior predictive samples.
y_ppc_t = pm.sample_posterior_predictive(idata_t, model_t)

ax = az.plot_ppc(y_ppc_t, num_pp_samples=100, mean=True)
ax.set_xlim(40, 70)