# Inferring prior distribution

Suppose that we take a sample ($X_1, X_2, ... X_n$) from a distribution where we assume that $X_i$~$Bin(n_i, p_i)$ and $n_i$ is known for every $i$. We also assume that $p_i$'s are independent and identically distributed, $p_i$ ~ $D$, where $D$ is some unknown distribution. $n_i$ cannot be assumed to be large.

My goal is to get a Bayesian estimate (or a probability distribution) for $p_i$. But this requires coming up with a distribution for $D$.

One option is to make an empirical distribution that uses frequentist estimates for each $p_i$ (i.e. $p_i = X_i/n_i$). This is a rather intuitive and potentially reasonable idea. Unfortunately, the presence of small $n_i$'s would make the tails heavier than they should be (lots of extreme values close to 0 or 1).

I'm looking for another option that doesn't have the problems of the aforementioned solution.

One possibility I have in mind is to use the following algorithm:

1. Generate prior distribution as explained earlier.
2. Get MAP or EAP estimate for every $p_i$.
3. Generate new empirical prior from the probabilities obtained in 2.
4. Go back to 2 (continue for a set number of steps, or possibly until convergence?)

Is this method similar to any method out there? Is it reasonable?

• One thing to keep in mind is that a prior is supposed to be (at least in principle) information not contained in the data but that you have available through some other channel (previous study, intuition, etc). Obviously that raises a bigger issue about empirical priors in general, but right now it sounds like you just don't have much outside info. So why not use some kind of weak prior? You could use the beta distribution, which happens to be conjugate, but is also very flexible. Then you're also free to control the tails and apply hyperpriors (e.g. the half-Cauchy). – shadowtalker Sep 3 '14 at 4:42
• This sounds like a hierarchical model. If I wanted to recreate the dataset, here's what I'd do: Let $D$ be $Beta(\alpha, \beta)$ (reasonable since we are dealing with probabilities). We don't know $\alpha, \beta$, so we assign priors to them, say exponential for both with some $\lambda$ hyperparameter. Then we draw the $p_i$ for each $i$, and sample $X_i$ from the binomials. Let me write something up... – Cam.Davidson.Pilon Sep 3 '14 at 18:03

I hope you like Python! I'll recite my comment here:

This sounds like a hierarchical model. If I wanted to recreate the dataset, here's what I'd do: Let $D$ be a $Beta(\alpha, \beta)$ distribution (reasonable since we are dealing with probabilities). We don't know $\alpha, \beta$, we assign priors to them, say exponential for both with some $\lambda$ hyperparameter. Then we draw the $p_i$ for each $i$, and sample $X_i$ from the binomials.

That's how I would recreate the dataset. To make inference, we go backwards. Here's the model in PyMC:

import pymc as pm

#fake data
X = np.array([3,2,2,5,7,10,11])
n = np.array([5, 4, 4, 6, 10, 19, 12])

#here I make sure I fulfill fake-data constraints
assert X.shape == n.shape
assert (X <= n).all()

alpha = pm.Exponential("alpha", 1)
beta = pm.Exponential("beta", 1)

p = pm.Beta( "p", alpha, beta, size=X.shape[0])

obs = pm.Binomial("obs", n, p, value=X, observed=True)

mcmc = pm.MCMC([obs,p,beta,alpha])
mcmc.sample(10000, 5000)


And some output:

With samples from the posteriors of $\alpha$ and $\beta$, we can reconstruct possible distributions of $D$, the unknown distribution:

Edit: Apologies, the x-axis should be between 0-1, not 0-500, it's Python thing I forgot to change.

• Here's the ipython %hist: gist.github.com/CamDavidsonPilon/1fed2295083e660b776a – Cam.Davidson.Pilon Sep 3 '14 at 18:23
• Thank you very much about your answer. Unfortunately I don't really have the math knowledge to really understand it yet. Is MCMC (which I don't know about) typically used for that kind of stuff? What specifically should I read about? You say you can "reconstruct possible distributions of D". Would it make sense to take an average of a bunch of reconstructions or would it not be any more valid than any one reconstruction? – Uwat Sep 4 '14 at 0:14
• @Uwat, when you have a even somewhat non-trivial model, the only way to perform Bayesian inference is with MCMC. I would recommend Bayesian Methods for Hackers (I'm the author btw) for an intro to Bayesian methods (chapter 1 & 2) and MCMC (chapter 3). For "reconstruct", it's best to not use an individual reconstruction. There are a few choices you could make: 1) use the average values of alpha, beta. 2) Use the MAP of alpha, beta (better). 3) Use the average over the distributions (also better) – Cam.Davidson.Pilon Sep 4 '14 at 2:42
• Assuming that $p_i ~ \stackrel{iid}{\sim}\mbox{Beta}(\alpha, \beta)$ seems like a very strong assumption, and if the ultimate goal is inference about $D$ I wouldn't recommend it as a default inference. Either a Dirichlet mixture of Beta distributions (as in Liu, AOS, 1996) or some nonparametric Empirical Bayes technique seems much better. There isn't really any reason to expect $D$ to be a simple 2 parameter distribution. – guy Sep 4 '14 at 14:55
• @Cam.Davidson.Pilon I disagree that a Beta can capture "most interesting distributions." For example, no Beta distribution is bimodal, which is the scenario in the data considered by Liu in the paper I mentioned. – guy Sep 5 '14 at 0:06

The algorithm that you describe is treating $D$ like a variable in the problem, but using a method other than Bayesian inference to deal with $D$. A Bayesian solution requires handling $D$ in a Bayesian manner, i.e. giving it a prior and integrating it out to get the marginal for $p_i$. For example, you could use a Dirichlet process as the prior for $D$ and do inference for $p_i$ via Gibbs sampling.