
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)
	σ_σ = pm.HalfNormal('σ_σ', sd=10)

	μ = pm.Normal('μ', mu=μ_σ, sd=σ_σ, shape=len(y))

	obs = pm.Normal('obs', mu=μ, sd=sigma, observed=y)

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

	# Posterior diagnostics
	summary = az.summary(idata)
