# MCMC Modelling - can this even be solved?

I am very new to Bayesian modelling and MCMC - I would like to know if the problem I describe below can be solved. It seems to be there is too much missing information but I wanted to get your thoughts. Consider the following:

I have a road intersection with one entrance and two exits, A and B. My goal is to estimate the number of cars that pass through this intersection in a given day, which equals to the number of cars that pass through the entrance. I post two people, one at exit A and the other at exit B, to count the number of cars coming out of their exits. They are both not very good so they only capture C percent of the cars that actually pass through their respective exits. I do not know what C is. To make matters worse, the person at exit B lost his record so I only have the numbers from A. I think this situation should be described by the following graph: From historical data, I know on average p% of ppl go through A and (1-p)% people go through B, but on this given day, I have no information. In this example I only have two exits, but in general I may have more (e.g. 3-5).

Is is possible to estimate the distribution for "# cars thru entrance" with the data I have? If so, what distributions would you assign to each random variable?

• By definition, $100*p$ percent of people must use A and $100*(1-p)$ percent of people must use B, unless there is a third exit, so that change in notation will simplify your graph. The same goes for number of cars thru A and B: they must add up to the total number. – Sycorax says Reinstate Monica Jun 1 '14 at 15:02

I'll record my modelling thoughts as I create it:

import pymc as pm


How can I recreate this data? Well, N_T cars enter. Do I know N_T? No, so it's a random variable. For simplicity, I'll say its a discrete uniform with max 1000 (you can change this to, say, a Poisson)

N_T = pm.DiscreteUniform('N_T', 0, 1000)


Of these N_T cars that enter, I know a car will enter A with probability p. Thus the number of cars through A is a binomial rv. But do I know p for certain? Probably not, you mentioned it was a historical average, so we should model this stochastically (using the beta-binomial model):

p_A = pm.Beta('p_A', 50, 50) #centered at 0.5, but with some variance. This depends on your historical data.
N_A = pm.Binomial('N_A', N_T, p_A)


Next, only C percent of cars are recorded by A's observer. Let's just assume that C is completely accurate, i.e. there is no uncertainty about what C is:

C=0.95


Then the number of cars recorded is again a binomial, with N_A trials. Furthermore, suppose observer A record 20 cars.

obs_A = pm.Binomial('obs', N_A, C, observed=True, value = 20)


Let's fire the mcmc canon:

mcmc = pm.MCMC( [obs_A, N_A, p_A, N_T] )
mcmc.sample(10000, 5000)

from pymc.Matplot import plot as mcplot
mcplot(mcmc)


And I get this posterior for N_T: • Fun question, btw =) – Cam.Davidson.Pilon Jun 1 '14 at 16:34
• Thanks for the reply and glad you found it interesting. All I know about this topic is from your EXCELLENT Bayesian Methods for Hackers book. I definitely need to re-read it again soon as I made a very quick pass (mainly because it was so engaging and readable) the first time. Any plans on completing chapters X1 and X2? – mchangun Jun 1 '14 at 16:44
• Yes, but they will be in a different medium. Details soon ;) – Cam.Davidson.Pilon Jun 1 '14 at 16:47
• What if I don't know the value of C and need to model that as a random variable as well? It's actually quite crucial to this problem that I don't assume I know it's value. Would I still be able to get reasonable answers? – mchangun Jun 1 '14 at 16:50
• Of course: you can make C a pymc stochastic variable: C = pm.Beta('C', some_alpha, some_beta). The implication of doing this would be wider posteriors for N_T. – Cam.Davidson.Pilon Jun 1 '14 at 17:20