I have two beta-distributions:

$H_1 = Beta(\alpha_1, \beta_1) $


$H_2 = Beta(\alpha_2, \beta_2) $

(parameters are known), and I'd like to estimate whether a new sample $D$ rather comes from $H_1$ or $H_2$.

It seems that Bayes factor is the solution

$K = \frac{\int \Pr(\theta_1|H_1)\Pr(D|\theta_1,H_1)\,d\theta_1} {\int \Pr(\theta_2|H_2)\Pr(D|\theta_2,H_2)\,d\theta_2}$

but I don't really see how to compute it. How do I get and integrate $Pr(\theta)$ and $Pr(D|\theta)$?

Also, how does it generalize to Dirichlets?

Thanks, JP.

  • $\begingroup$ What is the likelihood of your data? I'm assuming Binomial? $\endgroup$
    – user44764
    May 15 '14 at 14:22
  • $\begingroup$ @Matthew Well, actually my real data is multinomial, that's why I ask for the generalization, but I'd like first to understand the binomial/beta-distribution case. $\endgroup$
    – user45607
    May 15 '14 at 14:27

You said in the comments that this is for a Binomial likelihood, which you want to extend to the Multinomial-Dirichlet situation.

To evaluate $K$ you're looking to evaluate the likelihood of your data given each model. In the Beta-Binomial case of $N$ trials where $$ p \sim \text{Beta}(\alpha, \beta) \\ y \sim \text{Binomial}(N, p), $$ then the data likelihood is $$ \begin{eqnarray*} \pi(y \mid \alpha, \beta) &=& \int_0^1 \text{Binomial}(y \mid N, p) \text{Beta}(p \mid \alpha, \beta) \,dp. \\ &=& \int_0^1 {N \choose y} p^y (1-p)^{N-y} \times \frac{1}{\text{B}(\alpha,\beta)} p^{\alpha-1} (1-p)^{\beta-1} \,dp \\ &=& {N \choose y} \frac{1}{\text{B}(\alpha,\beta)} \int_0^1 p^{\alpha - 1 + y} (1-p)^{\beta - 1 + N -y} \,dp \\ &=& {N \choose y} \frac{\text{B}(\alpha+y,\beta + N - y)}{\text{B}(\alpha,\beta)} \end{eqnarray*} $$

where $\text{B}$ is the beta function; you should be able to compute it with some software package or even some online tool. This probability distribution is called the Beta-Binomial distribution. In this case, $$ K = \frac{\text{B}(\alpha_1 +y,\beta_1 + N - y) \text{B}(\alpha_2,\beta_2)}{ \text{B}(\alpha_2 +y,\beta_2 + N - y) \text{B}(\alpha_1,\beta_1)} $$

In the multiple category case there is the Dirichlet-Multinomial distribution which gives $K$ in a similar form with Gamma functions.

  • 1
    $\begingroup$ Thanks a lot! This answer exactly covers what I missed from my readings so far, that's why I mark it as best answer. But new readers should definitely read catchmeifyoutry's answer as well, which is perfectly complementary. $\endgroup$
    – user45607
    May 15 '14 at 15:58
  • $\begingroup$ @Matthew: Shouldn't the data likelihood as you computed it be the numerator of your result? The normalizing constant of the posterior beta distribution is B(a+y,b+N-y). In my opinion, it is not a ratio of beta functions $\endgroup$
    – beginneR
    Oct 30 '14 at 17:24
  • $\begingroup$ @beginneR The p.m.f. for the Beta-Binomial distribution is as I derived. If you instead mean that I have incorrectly determined the model likelihood, then observe from the law of total probability that $P(D\mid M) = \int_\mathbb{R} P(D \mid \theta) p(\theta \mid M) d\theta$, where $p(\theta \mid M)$ is just the prior probability of the parameter $\theta$ under the model $M$. $\endgroup$
    – user44764
    Oct 30 '14 at 20:51
  • $\begingroup$ Maybe I am missing something here. Is it correct that you derived the denominator of Bayes Law for binomial data, also called the prior predictive distribution? I thought that the BF is a ratio of the denominator terms of two models (p(D|M1)/p(D|M2)). Or, put differently, the BF is the ratio of two normalizing constants for the posterior distributions? $\endgroup$
    – beginneR
    Nov 1 '14 at 12:34
  • $\begingroup$ @beginneR The models being discussed here are probability distributions on the parameters; that is, $M_1 \equiv p \sim \text{Beta}(\alpha_1, \beta_1)$ and $M_2 \equiv p \sim \text{Beta}(\alpha_2, \beta_2)$. So in order to obtain $P(M_i \mid D)$ it is necessary to integrate over the parameter $p$ as shown $\endgroup$
    – user44764
    Nov 3 '14 at 3:44

I'll answer for the more general case of a Dirichlet prior and Multinomial likelihood, which is what you are interested in anyway.

The integral has a closed form solution, which can be found in this Wikipedia article on the Drichlet-Multinomial.

$$p(X|\alpha_0) = DirMult(X| \alpha_0) = \int_\theta Cat(X|\theta) Dir(\theta|\alpha_0) d \theta$$ $$= \frac{\Gamma(\sum_k \alpha_{0k})}{\Gamma(\sum_k n_k + \sum_k \alpha_{0k})} \prod_k \frac{\Gamma(n_k + \alpha_{0k})}{\Gamma(\alpha_{0k})}$$

where X are the observations, which we can summarize into occurrence counts $n_k$ ($n_i$ is the number of times event $i$ occurred in observations $X$), $a0$ is the vector of parameters of the Dirichlet prior. $\Gamma$ is the Gamma function.

An important but subtitle thing to note is that this solution actually assumes a Categorical likelihood instead of Multinomial, as the name suggests. That means that for a proper multinomial, you should include the multinomial coefficient.


It can be good to study how to obtain this result when confronted with other conjugate pairs

  • From the Dirichlet distribution we know that

$$\int_\theta \left [ \frac{\Gamma(\sum_k \alpha_{0k})}{\prod_k \Gamma(\alpha_{0k})} \prod_k \theta_k^{(\alpha_{0k} - 1)} \right ] d\theta = 1$$ thus $$\int_{\theta} \left [ \prod_k \theta_k^{(\alpha_k - 1)} \right ] d\theta = \frac{\prod_k \Gamma(\alpha_k)} {\Gamma(\sum_k \alpha_k)}$$

  • Therefore,

$$p(X|\alpha_0) = \int_\theta Cat(X|\theta) Dir(\theta|\alpha_0) d \theta$$

$$= \int_\theta \left [ \prod_k \theta_k^{n_k} \times \frac{\Gamma(\sum_k \alpha_{0k})}{\prod_k \Gamma(\alpha_{0k})} \prod_k \theta_k^{(\alpha_{0k} - 1)} \right ] d\theta$$

$$= \frac{\Gamma(\sum_k \alpha_{0k})}{\prod_k \Gamma(\alpha_{0k})} \int_\theta \left [ \prod_k \theta_k^{(n_k + \alpha_{0k} - 1)} \right ] d\theta$$

$$= \frac{\Gamma(\sum_k \alpha_{0k})}{\prod_k \Gamma(\alpha_{0k})} \times \frac{\prod_k \Gamma(\alpha_{0k} + n_k)} {\Gamma(\sum_k \alpha_{0k} + \sum_k n_k)}$$

$$= \frac{\Gamma(\sum_k \alpha_{0k})}{\Gamma(\sum_k n_k + \sum_k \alpha_{0k})} \prod_k \frac{\Gamma(n_k + \alpha_{0k})}{\Gamma(\alpha_{0k})}$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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