Why is posterior same as prior despite observing data?

Questions
1

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).

2

Hi, it’s best to not use python if statements in PyMC

There are numerous ways to achieve what you want, here is one

with pm.Model() as model:
    weekend = pm.Bernoulli('weekend', 2/7)
    weekday_rate = 10
    weekend_rate = 3
    rate = weekend*weekend_rate + (1-weekend)*weekday_rate
    bus = pm.Poisson('buses', mu=rate, observed=4)
    trace = pm.sample(1000)

print('Is Weekday? {:.3f}'.format(np.mean(trace['weekend']==0)))
print('Is Weekend? {:.3f}'.format(np.mean(trace['weekend']==1)))
2 Likes
3

Thanks, this works!

I was also able to use theano's ifelse:

from theano.ifelse import ifelse

with pm.Model() as model:
    weekend = pm.Bernoulli('weekend', 2/7)
    rate = ifelse(weekend, 3, 10)
    bus = pm.Poisson('buses', mu=rate, observed=4)
    trace = pm.sample(10000)

Indeed, the docs for Deterministic say that the var parameter should be a “theano variable”:

pymc3.model.Deterministic(name, var, model=None)¶
Create a named deterministic variable

Parameters
name str
var theano variables

But perhaps it should feature more prominently in the docs (even if this is self-evident for someone who understands how PyMC works under the hood).

1 Like