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.

  • $\begingroup$ Try asking on stats.stackexchange.com - actually that seems to be my answer to a lot of qs on here.... problem? $\endgroup$
    – Spacedman
    Dec 7, 2011 at 10:47
  • $\begingroup$ @Spacedman I think the two overlap nicely - I've started suggesting for a small number of qs over there that they come here. $\endgroup$
    – Fomite
    Dec 7, 2011 at 13:40
  • $\begingroup$ I thought the idea with SE was that overlapping fields would live on one site with a tag. I'm wondering if this scicomp site is a good idea now... We'll see. $\endgroup$
    – Spacedman
    Dec 7, 2011 at 13:56
  • 1
    $\begingroup$ @Spacedman: Well, I thought the idea was for this site to be a sort of meta-computation site akin to SO for programming. I'd be disappointed if it turns into a purely applied maths affair. Maybe part of the reason is that currently applied math people have nowhere else to go on SE. Which makes me wonder if lots of bioinformatics people are going to show up too. I just posted on stats.sx chat suggesting people visit here. Maybe I should post on meta stats too. $\endgroup$ Dec 7, 2011 at 15:11
  • 2
    $\begingroup$ @Spacedman I think this fits squarely in the realm for this site since the question is "Can anyone suggest methods, preferably as efficient and simple as possible, to simulate from" a given probability distribution. That is a computational questions (though one that does come up in statistics). $\endgroup$ Dec 7, 2011 at 20:30

1 Answer 1


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


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.

  • 1
    $\begingroup$ Hmm, I thought someone might say MCMC. I was hoping for something else. MCMC is notoriously computation intensive. Any idea what distribution this is? Also, any tips about doing it specifically in this case? $\endgroup$ Dec 7, 2011 at 15:04
  • $\begingroup$ @FaheemMitha Ah, I did miss the "preferably as efficient and simple as possible" part of your question, so I see how MCMC might not be a great answer. As for your other questions, I'll think more about them, but unfortunately, as my "(only?)" remark suggests, I may not be much help. $\endgroup$
    – Barron
    Dec 7, 2011 at 15:46
  • $\begingroup$ Sure, no problem. Thanks for answering. I actually did my thesis in an MCMC related area (perfect sampling), so I know it is not that easy to use in practice, and to be honest, tend to run in the opposite direction when it is mentioned.:-) But maybe this question is easier than my thesis topic - I'll take a closer look. $\endgroup$ Dec 7, 2011 at 15:55
  • $\begingroup$ @FahheemMitha If you do use MCMC, depending on what kind of software you're using to perform the computation, you can sample directly from the posterior. See confounding.net/2011/08/19/…, though it's fairly inefficient as these things go. $\endgroup$
    – Fomite
    Dec 7, 2011 at 16:30
  • 1
    $\begingroup$ The posterior $p(\mathbf{p}|y)$ is Dirichlet, no need to use simulation for that. If you are suggesting generating from $p(\mathbf{p}|y)$, then plugging the simulated thetas into $p(\hat{y}|\mathbf{p})$ to get values from $p(\hat{y}|y)$, that's what I'm currently doing (sorry, I should have mentioned that) but it is very ugly and inefficient. Anyway, this approach doesn't actually need MCMC at all, as far as I can see. $\endgroup$ Dec 7, 2011 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.