Sampling from posterior predictive distribution

Question

First post. I'm working on this problem using Bayesian methods. In desperation I'm considering using p-values (shock horror), specifically posterior predictive p-values. So I need to simulate from the posterior predictive. I don't see a lot of statistics stuff here - it's mostly applied math at the moment. I could take this to stats.sx I suppose, but this is definitely a computational problem, so I thought I'd give it a shot here. If there are any Bayesians lurking here, now would be a good time to delurk. :-)

In what follows, the hats denote predictive values. Consider

Let $p(\mathbf{p})$ be a Dirichlet prior.

\begin{align*} p(\mathbf{p}) = (k-1)! I(\sum p_i = 1) \end{align*}

Let the likelihood be

\begin{align*} p(y| \mathbf{p}) = \prod_{i=1}^k p_i^{f_i} \end{align*}

This is the likelihood for a sequence of categorical variables. Assume the prior predictive is (calculation omitted)

\begin{align*} p(y) = (k-1)! \frac{\prod_i f_i!}{(N + k - 1)!} \end{align*}

We obtain

\begin{align*} p(\mathbf{p}|y) = p(\mathbf{p}) p(y| \mathbf{p}) / p(y) = \frac{ (k-1)! I(\sum p_i = 1) \prod_{i=1}^k p_i^{f_i} } { (k-1)! \frac{\prod_i f_i!}{(N + k - 1)!} } \end{align*}

Dividing we obtain the posterior

\begin{align*} p(\mathbf{p}|y) = \frac{ I(\sum p_i = 1) \prod_{i=1}^k p_i^{f_i} (N + k - 1)! }{ \prod_i f_i! } \end{align*}

Define $f^{\prime}_i = f_i + \hat{f}_i$, and $N^{\prime} = N+\hat{N}$. Then the posterior predictive is (calculation omitted but can be added if anyone cares)

\begin{align*} p(\hat{y}|y) = \int p(\hat{y}|\mathbf{p}) p(\mathbf{p}|y) d\mathbf{p} = \frac{ \prod_i f^{\prime}_i! }{ \prod_i f_i! } \frac{ (N + k - 1)! }{(N^{\prime}+k-1)!} \end{align*}

So, can anyone suggest methods, preferably as efficient and simple as possible, to simulate from $p(\tilde{y}|y)$? This is a finite distribution, but too big to enumerate. At the moment, I'm not even sure what distribution it is. Thanks in advance.

PS. Can some kind person create "bayesian-computation" and "statistics" tags, and tag this message with them? "random" is not very useful. Thanks.

Answer 1

I believe the best (only?) way to accomplish this will be through Markov Chain Monte Carlo sampling:

http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo

This will allow you to sample from the unknown posterior without knowing it's functional form.

---

These are really more like questions than a better answer, but this seemed the best place to put them:

First, you can reformulate the problem in terms of only $k-1$ unknowns since the last one is functionally dependent (the $p$'s must sum to one).

Assuming $f_i$ represents the observed frequency of each case (please correct that assumption if it's wrong), why include it as part of the prior? I'd suggest that your prior should not depend on the data. You might use an exponent of 0 to represent a vague prior instead.

Unless $k$ is huge, this problem is easily solved by MCMC. You could set it up in Winbugs or PyMC if you don't want to write your own code. Also, with MCMC, the above derivations aren't needed. You only need to specify the prior (which may just be constant) and the likelihood (which is pretty simple).

I'd also suggest not trying to simulate the posterior predictive directly, but instead simulate the posterior first and then the predictive is easily obtained from that by simulation.