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

# Given Data
N = 12  # Number of observations
n = np.array([47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360])
r = np.array([0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24])

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
	mu = alpha + beta * r
	likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=n)

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

	# Posterior diagnostics
	summary = az.summary(trace)