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

# Data
y = np.array([28, 8, -3, 7, -1, 1, 18, 12])
sigma = np.array([15, 10, 16, 11, 9, 11, 10, 18])

with pm.Model() as m:
	# Priors
	mu = pm.Normal('mu', mu=0, sigma=10)
	tau = pm.HalfCauchy('tau', beta=10)

	# Likelihood
	theta = pm.Normal('theta', mu=mu, sigma=tau, shape=len(y))
	y_obs = pm.Normal('y_obs', mu=theta, sigma=sigma, observed=y)

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

	# Posterior diagnostics
	summary = az.summary(trace)