
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:
	μ_σ = pm.Normal('μ_σ', mu=0, sd=10, shape=8)
	μ = pm.Normal('μ', mu=μ_σ, sd=σ, observed=y)
	σ = pm.HalfNormal('σ', sd=1)

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

	# Posterior diagnostics
	summary = az.summary(trace)