I have the following model. (I’ll provide a short “story” which describes what is being modelled, not just the raw model)
Suppose we don’t know what day of the week it is. But we know that the number of buses per hour on the street is distributed according to a Poisson distribution, with rate 10 on weekdays and 3 on the weekends. We observe 4 buses in an hour. What is the probability that it is the weekend?
This is my PyMC3 model:
import pymc3 as pm import numpy as np with pm.Model() as model: weekend = pm.Bernoulli("weekend", 2/7) rate = 3 if weekend else 10 bus = pm.Poisson("buses", mu=rate, observed=4) trace = pm.sample(1000) np.mean(trace["weekend"])
When I run it, I get the output
Multiprocess sampling (4 chains in 4 jobs) BinaryGibbsMetropolis: [weekend] Sampling 4 chains, 0 divergences: 100%|██████████| 6000/6000 [00:01<00:00, 4433.79draws/s] 0.2785
where 0.2785 is about 2/7. In other words, the trace contains samples from the prior distribution of whether it is the weekend (Bernoulli with p=2/7).
How do I obtain samples from the posterior distribution (after observing 4 from either of the Poisson distributions)?
I am new to PyMC3 and probabilistic programming (I only followed the first few chapters of Bayesian Methods for Hackers).