# Basic Bayesian MCMC to estimate two parameters from binomial distributions given unknown number of trials

This is a very basic question about Bayesian inference. I'm not grasping one or more fundamental concepts.

Let's say I have two observed outcomes, X and Y. I want to infer the probabilities (px and py, respectively) of each occurring given X and Y. I do not know N, the total number of trials. I'm assuming that X and Y are binomially distributed. How do I calculate the likelihood without N?

What I ultimately want is to show the bivariate posterior distribution of px and py via MCMC. I do not care about estimating N--I just want to show the chain in the plane of px and py. No convergence necessary.

Clarification: X and Y are drawn from the same N: X ~Binom(N, px) and Y ~Binom(N, py). We have no other information about px or py, although I'll use a beta prior to start. Also assume that X and Y are independent.

• Are X and Y both drawn from the same distibution $Bin(N, p)$? And do we have any information about $p$? – Cam.Davidson.Pilon Jan 7 '13 at 21:30
• Good questions! Revised post to clarify. – neophyte Jan 7 '13 at 21:32
• And are you familiar with Python? – Cam.Davidson.Pilon Jan 7 '13 at 21:37
• Not Python, but C++, Java, Javascript, etc. I am going to use Javascript for this. It might be best to use pseudocode here. – neophyte Jan 7 '13 at 21:37
• It seems to me like the problem is under-determined. You're trying to guess 3 unknowns (px, py, and N) from 2 data points, with no additional information. I don't think there's any way around that. – David J. Harris Jan 7 '13 at 21:38

## Model and Pseudocode

So I did some analysis in Python, though I used the pyMC library which hides all the MCMC mathy stuff. I'll show you how I modeled it in semi-pseudocode, and the results.

I set my observed data as $X=5, Y=10$.

X = 5
Y = 10


I assumed that $N$ has a Poisson prior, with the Poisson's rate a $EXP(1)$. This is a pretty fair prior. Though I could have chosen some uniform distribution on some interval:

rate = Exponential( mu = 1 )
N = Poisson( rate = rate)


You mention beta priors on $pX$ and $pY$, so I coded:

pX = Beta(1,1) #equivalent to a uniform
pY = Beta(1,1)


And I combine it all:

observed = Binomial(n = N, p = [pX, pY], value = [X, Y] )


Then I perform the MCMC over 50000 samples, burned-in about half of that. Below are the plots I generated after MCMC.

## Interpretation:

Let's examine the first graph for $N$. The N Trace graph are the samples, in order, I generated from the posterior distribution. The N acorr graph is the auto-correlation between samples. Perhaps there is still too much auto-correlation, and I should burn-in more. Finally, N-hist is the histogram of posterior samples. It looks like the mean is 13. Notice too that no samples were drawn from below 10. This is a good sign, as that would be impossible given the data observed was 5 and 10.

Similar observations can be made for the $pX$ and $pY$ graphs.

## Different Prior on $N$

If we restrict $N$ to be a Poisson( 20 ) random variable (and remove the Exponential heirarchy), we get different results. This is an important consideration, and reveals that the prior can make a large difference. See the plots below. Note the time to convergence was much larger here too.

On the other hand, using a Poisson( 10 ) prior produced similar results to the Exp. rate prior.

• I should note that I tried the same analysis with an uninformative prior on $N$ and got terrible results (chain did not converge to anything meaningful, and autocorrelations spiked). – Cam.Davidson.Pilon Jan 7 '13 at 22:11
• This looks really good--thanks. I hadn't been able to figure out what to do with N, but using a Poisson prior with some realism (i.e., >max(X,Y)) and fitting it too makes sense... and makes it easy to see how I'll calculate the likelihood. :P I'm going to try coding this up in js tomorrow and will almost certainly accept your answer. Thanks again. – neophyte Jan 7 '13 at 22:15
• +1 @Cam.Davidson.Pilon. This seems to have worked better than I expected. What happens if you give it a more reasonable prior like Poisson(10) or Poisson(20) that has more probability mass in the feasible region above 10? – David J. Harris Jan 7 '13 at 23:35
• With a prior of Poisson(10), the results look the same (and convergence is quick), not so with Poisson(20). See edit above. – Cam.Davidson.Pilon Jan 8 '13 at 13:14
• Bayesian methods are quite difficult in the details. The overall picture is easy: Return a distribution of the parameter in question, given the data and prior, but the small details like: how to sample, what to do with samples, what do I need to calculate and so on can be intimidating. It took me three straight days staring at a computer screen to finally break through ahha – Cam.Davidson.Pilon Jan 9 '13 at 13:16