pymc3_bayesian_linear_regression_map_mcmc_sampling.py

1import pymc3 as pm
2import numpy as np
3import matplotlib.pyplot as plt
4
5# Generating simulated data
6np.random.seed(123)
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# Building the model
18with pm.Model() as model:
19    # Priors for unknown model parameters
20    sigma = pm.HalfNormal("sigma", sd=1)
21    intercept = pm.Normal("Intercept", mu=0, sd=20)
22    x_coeff = pm.Normal("x", mu=0, sd=20)
23
24    # Expected value of outcome
25    likelihood = intercept + x_coeff * x
26
27    # Likelihood (sampling distribution) of observations
28    y_obs = pm.Normal("y_obs", mu=likelihood, sd=sigma, observed=y)
29
30    # Obtain posterior samples using MCMC
31    trace = pm.sample(2000, return_inferencedata=True)
32
33# Plotting the results
34import arviz as az
35az.plot_trace(trace)
36plt.show()
37
38# Print summary statistics
39print(az.summary(trace, round_to=2))