pymc_bayesian_linear_regression_mcmc_sampling_arviz_visualization.py

1import arviz as az
2import matplotlib.pyplot as plt
3import numpy as np
4import pymc as pm
5
6# Generate synthetic data
7size = 200
8true_intercept = 1
9true_slope = 2
10
11x = np.linspace(0, 1, size)
12# y = a + b*x + e
13true_regression_line = true_intercept + true_slope * x
14# add noise
15y = true_regression_line + np.random.normal(scale=0.5, size=size)
16
17# Define the model
18with pm.Model() as model:
19    # Priors for unknown model parameters
20    sigma = pm.HalfCauchy("sigma", beta=10)
21    intercept = pm.Normal("Intercept", 0, sigma=20)
22    slope = pm.Normal("slope", 0, sigma=20)
23
24    # Expected value of outcome
25    likelihood = pm.Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y)
26
27    # Inference!
28    # draw 1000 posterior samples using NUTS sampling
29    trace = pm.sample(1000, chains=4, cores=4)
30
31# Plot the results
32az.plot_trace(trace)
33plt.show()
34
35# Summarize the posterior distribution
36print(az.summary(trace, round_to=2))