Bayesian Polynomial Regression#
%pip install pymc pytensor
Requirement already satisfied: pymc in /usr/local/lib/python3.10/dist-packages (5.7.2)
Requirement already satisfied: pytensor in /usr/local/lib/python3.10/dist-packages (2.14.2)
Requirement already satisfied: arviz>=0.13.0 in /usr/local/lib/python3.10/dist-packages (from pymc) (0.15.1)
Requirement already satisfied: cachetools>=4.2.1 in /usr/local/lib/python3.10/dist-packages (from pymc) (5.3.2)
Requirement already satisfied: cloudpickle in /usr/local/lib/python3.10/dist-packages (from pymc) (2.2.1)
Requirement already satisfied: fastprogress>=0.2.0 in /usr/local/lib/python3.10/dist-packages (from pymc) (1.0.3)
Requirement already satisfied: numpy>=1.15.0 in /usr/local/lib/python3.10/dist-packages (from pymc) (1.25.2)
Requirement already satisfied: pandas>=0.24.0 in /usr/local/lib/python3.10/dist-packages (from pymc) (1.5.3)
Requirement already satisfied: scipy>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from pymc) (1.11.4)
Requirement already satisfied: typing-extensions>=3.7.4 in /usr/local/lib/python3.10/dist-packages (from pymc) (4.9.0)
Requirement already satisfied: setuptools>=48.0.0 in /usr/local/lib/python3.10/dist-packages (from pytensor) (67.7.2)
Requirement already satisfied: filelock in /usr/local/lib/python3.10/dist-packages (from pytensor) (3.13.1)
Requirement already satisfied: etuples in /usr/local/lib/python3.10/dist-packages (from pytensor) (0.3.9)
Requirement already satisfied: logical-unification in /usr/local/lib/python3.10/dist-packages (from pytensor) (0.4.6)
Requirement already satisfied: miniKanren in /usr/local/lib/python3.10/dist-packages (from pytensor) (1.0.3)
Requirement already satisfied: cons in /usr/local/lib/python3.10/dist-packages (from pytensor) (0.4.6)
Requirement already satisfied: matplotlib>=3.2 in /usr/local/lib/python3.10/dist-packages (from arviz>=0.13.0->pymc) (3.7.1)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from arviz>=0.13.0->pymc) (23.2)
Requirement already satisfied: xarray>=0.21.0 in /usr/local/lib/python3.10/dist-packages (from arviz>=0.13.0->pymc) (2023.7.0)
Requirement already satisfied: h5netcdf>=1.0.2 in /usr/local/lib/python3.10/dist-packages (from arviz>=0.13.0->pymc) (1.3.0)
Requirement already satisfied: xarray-einstats>=0.3 in /usr/local/lib/python3.10/dist-packages (from arviz>=0.13.0->pymc) (0.7.0)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=0.24.0->pymc) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=0.24.0->pymc) (2023.4)
Requirement already satisfied: toolz in /usr/local/lib/python3.10/dist-packages (from logical-unification->pytensor) (0.12.1)
Requirement already satisfied: multipledispatch in /usr/local/lib/python3.10/dist-packages (from logical-unification->pytensor) (1.0.0)
Requirement already satisfied: h5py in /usr/local/lib/python3.10/dist-packages (from h5netcdf>=1.0.2->arviz>=0.13.0->pymc) (3.9.0)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (4.48.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (1.4.5)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.2->arviz>=0.13.0->pymc) (3.1.1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.1->pandas>=0.24.0->pymc) (1.16.0)
import numpy as np
import arviz as az
import matplotlib.pyplot as plt
import requests
import io
import csv
import pandas as pd
import pymc as pm
target_url = 'https://raw.githubusercontent.com/cfteach/brds/main/datasets/anscombe.csv'
download = requests.get(target_url).content
ans = pd.read_csv(io.StringIO(download.decode('utf-8')))
x = ans[ans.group == 'II']['x'].values
y = ans[ans.group == 'II']['y'].values
x = x - x.mean()
f=plt.figure()
f.set_figwidth(4)
f.set_figheight(2)
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x,y,'C0.',alpha=0.6, markersize=10)
[<matplotlib.lines.Line2D at 0x7e4a779e3100>]
with pm.Model() as model_poly:
alpha = pm.Normal('alpha',mu=y.mean(),sigma=1)
beta1 = pm.Normal('beta1',mu=0.,sigma=1)
beta2 = pm.Normal('beta2',mu=0.,sigma=1)
epsilon = pm.HalfCauchy('epsilon', 5)
mu = pm.Deterministic('mu',alpha + beta1*x + beta2*x*x)
y_pred = pm.Normal('y_pred', mu=mu, sigma = epsilon, observed = y)
idata_poly = pm.sample(2000, tune = 2000, return_inferencedata=True)
100.00% [4000/4000 00:09<00:00 Sampling chain 0, 0 divergences]
100.00% [4000/4000 00:09<00:00 Sampling chain 1, 0 divergences]
az.plot_trace(idata_poly, var_names = ['alpha','beta1','beta2'])
array([[<Axes: title={'center': 'alpha'}>,
<Axes: title={'center': 'alpha'}>],
[<Axes: title={'center': 'beta1'}>,
<Axes: title={'center': 'beta1'}>],
[<Axes: title={'center': 'beta2'}>,
<Axes: title={'center': 'beta2'}>]], dtype=object)
Perturbed Dataset#
yn = y.copy()
yn += np.random.normal(0,.5, len(yn))
f = plt.figure(figsize=(7,3))
plt.plot(x,yn,'C0.')
[<matplotlib.lines.Line2D at 0x7e4a79047850>]
with pm.Model() as model_poly2:
alpha = pm.Normal('alpha',mu=y.mean(),sigma=1)
beta1 = pm.Normal('beta1',mu=0.,sigma=1)
beta2 = pm.Normal('beta2',mu=0.,sigma=1)
epsilon = pm.HalfCauchy('epsilon', 5)
mu = pm.Deterministic('mu',alpha + beta1*x + beta2*x*x)
y_pred_per = pm.Normal('y_pred_per', mu=mu, sigma = epsilon, observed = yn)
idata_poly_per = pm.sample(2000, tune = 2000, return_inferencedata=True)
100.00% [4000/4000 00:08<00:00 Sampling chain 0, 0 divergences]
100.00% [4000/4000 00:07<00:00 Sampling chain 1, 0 divergences]
az.plot_trace(idata_poly_per)
array([[<Axes: title={'center': 'alpha'}>,
<Axes: title={'center': 'alpha'}>],
[<Axes: title={'center': 'beta1'}>,
<Axes: title={'center': 'beta1'}>],
[<Axes: title={'center': 'beta2'}>,
<Axes: title={'center': 'beta2'}>],
[<Axes: title={'center': 'epsilon'}>,
<Axes: title={'center': 'epsilon'}>],
[<Axes: title={'center': 'mu'}>, <Axes: title={'center': 'mu'}>]],
dtype=object)
res2 = az.summary(idata_poly_per)
print(res2)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \
alpha 8.926 0.284 8.394 9.458 0.007 0.005 1673.0
beta1 0.453 0.060 0.340 0.566 0.001 0.001 2647.0
beta2 -0.129 0.021 -0.168 -0.090 0.001 0.000 1688.0
epsilon 0.613 0.196 0.315 0.962 0.005 0.004 1716.0
mu[0] 9.250 0.277 8.716 9.742 0.007 0.005 1709.0
mu[1] 8.343 0.273 7.856 8.892 0.007 0.005 1785.0
mu[2] 8.672 0.336 8.004 9.252 0.006 0.005 2892.0
mu[3] 8.926 0.284 8.394 9.458 0.007 0.005 1673.0
mu[4] 9.315 0.262 8.835 9.794 0.006 0.004 1958.0
mu[5] 7.962 0.480 7.057 8.867 0.010 0.007 2518.0
mu[6] 6.404 0.256 5.923 6.892 0.005 0.003 3055.0
mu[7] 3.431 0.475 2.546 4.362 0.009 0.007 2557.0
mu[8] 9.123 0.267 8.633 9.628 0.005 0.004 2635.0
mu[9] 7.503 0.253 7.017 7.968 0.005 0.004 2176.0
mu[10] 5.047 0.327 4.398 5.643 0.006 0.004 3316.0
ess_tail r_hat
alpha 1792.0 1.0
beta1 1945.0 1.0
beta2 1907.0 1.0
epsilon 1940.0 1.0
mu[0] 1704.0 1.0
mu[1] 1831.0 1.0
mu[2] 2202.0 1.0
mu[3] 1792.0 1.0
mu[4] 1978.0 1.0
mu[5] 2528.0 1.0
mu[6] 2512.0 1.0
mu[7] 2535.0 1.0
mu[8] 1866.0 1.0
mu[9] 2145.0 1.0
mu[10] 2717.0 1.0
ppc = pm.sample_posterior_predictive(idata_poly_per, model=model_poly2)
100.00% [4000/4000 00:00<00:00]
plt.plot(x,yn,'C0.')
az.plot_hdi(x,idata_poly_per.posterior['mu'],color='C4',hdi_prob=.68)
az.plot_hdi(x,idata_poly_per.posterior['mu'],color='C2',hdi_prob=.95)
<Axes: >
plt.plot(x,yn,'C0.')
az.plot_hdi(x,ppc.posterior_predictive['y_pred_per'],color='C4',hdi_prob=.68)
az.plot_hdi(x,ppc.posterior_predictive['y_pred_per'],color='C2',hdi_prob=.95)
<Axes: >
Notice there is a difference between the “uncertainty band” built with the posterior distribution and with the posterior predictive distribution, the latter being larger.
Both the posterior and the posterior predictive distributions incorporate uncertainty due to noise in the data, but they do so in different ways. The posterior reflects this uncertainty in the context of parameter estimation, while the posterior predictive distribution goes a step further, and reflects it in the context of predicting new, unseen outcomes.
Posterior predictive distribution inherently includes both the uncertainty in the parameter estimates and the noise in the data.
import xarray as xr
tmp_y_pred = ppc.posterior_predictive['y_pred_per']
y_pred_merged = xr.concat([tmp_y_pred.sel(chain=i) for i in tmp_y_pred.chain], dim='draw')
# In the context of Bayesian statistics and probabilistic programming (such as when using libraries like PyMC or xarray for
# handling posterior distributions), the term "draw" typically refers to a single sample from the posterior distribution.
np.shape(y_pred_merged)
(4000, 11)
y_pred_means = y_pred_merged.mean(axis=0)
y_pred_std = y_pred_merged.std(axis=0)
print(np.shape(y_pred_means), np.shape(y_pred_std))
(11,) (11,)
y
array([9.14, 8.14, 8.74, 8.77, 9.26, 8.1 , 6.13, 3.1 , 9.13, 7.26, 4.74])
dx = np.linspace(0,11,100)
dy = dx
# Plotting
plt.figure(figsize=(6.5, 6))
plt.errorbar(y, y_pred_means, yerr=y_pred_std, fmt='o', label='Predicted (with uncertainty)', alpha=0.7)
plt.plot(dx,dy,'--',color='r', label = "")
plt.scatter(y, yn, color='r', label='Observed')
plt.xlabel('observed y')
plt.ylabel('predicted y')
plt.title('Predicted vs. Observed Y Values')
plt.xlim((0,11))
plt.ylim((0,11))
plt.legend()
plt.show()