Inferring prior distribution

Question

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?

Answer 1

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.

Answer 2

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.