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

# Data
X = np.array([1, 1.5, 1.5, 1.5, 2.5, 4, 5, 5, 7, 8, 8.5, 9, 9.5, 9.5, 10, 12, 12, 13, 13, 14.5, 15.5, 15.5, 16.5, 17, 22.5, 29, 31.5])
y = np.array([1.8, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47, 2.19, 2.26, 2.4, 2.39, 2.41, 2.5, 2.32, 2.32, 2.43, 2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.7, 2.72, 2.57])

with pm.Model() as m:
	# Define priors
	alpha = pm.Normal('alpha', mu=0, sigma=10)
	beta = pm.Normal('beta', mu=0, sigma=10)
	sigma = pm.HalfNormal('sigma', sigma=10)

	# Define likelihood
	y_obs = pm.Normal('y_obs', mu=alpha + beta * X, sigma=sigma, observed=y)

	# Sample posterior
	trace = pm.sample(1000, tune=1000, chains=4, return_inferencedata=True, idata_kwargs={"log_likelihood": True})

	# Posterior diagnostics
	summary = az.summary(trace)