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

# Given Data
N = 11  # Number of observations
x = np.array([-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10])
y = np.array([4.75906, 1.59423, 2.99548, 5.27501, 1.66472, 2.24347, 2.8914, 4.08681, 4.60588, 0.802364, 3.92136])
k = np.array([40, 37, 29, 12, 4, 3, 9, 19, 77, 82, 33])

with pm.Model() as m:
	# Define the prior distributions for the parameters
	alpha = pm.HalfNormal("alpha", sigma=1.0)
	beta = pm.Normal("beta", mu=0, sigma=1.0)
	sigma = pm.HalfNormal("sigma", sigma=1.0)

	# Define the likelihood function
	mu = alpha * np.exp(beta * x)
	likelihood = pm.Poisson("likelihood", mu=mu, observed=y)

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

	# Posterior diagnostics
	summary = az.summary(trace)