import pymc as pm
import corner
import numpy as np
import arviz as az
from pymc.distributions.transforms import Interval
The last line is only needed for some of the codes
For this example, we fit the model from chapter 6.1.2. First of all, we need to read the data:
(obstot, obsbkg, C) = np.loadtxt('/Dataset/data6.1.2.dat', delimiter = ",", unpack = True)
The full model is specified in the library:
with pm.Model() as model:
pm.Uniform("s", lower = 0, upper = 1.0e7)
pm.Uniform("bkg", lower = 0, upper = 1.0e7)
pm.Poisson("obstot", mu = (model['s'] + model['bkg']/C), observed = obstot)
pm.Poisson("obsbkg", mu = model['bkg'], observed = obsbkg)
To sample the posterior distributions we use this instruction:
chain = pm.sample(return_inferencedata = True)
To get a summary of the posterior you can use:
print(az.summary(chain.posterior['s']))
print(az.summary(chain.posterior['bkg']))
To draw a plot of the posterior you can use:
figure = corner.corner(chain.posterior['s'][0])
figure = corner.corner(chain.posterior['bkg'][0])