Hi,

I am working on writing a Gibbs sampler to do hypothesis testing in a Bayesian framework. I have a set of data which give the number of occurrences $\displaystyle x$ in a sample of $\displaystyle n$ trials. The particular model of interest is a three-layered hierarchical model with the following distributions:

1) The number of occurrences $\displaystyle x$ given $\displaystyle p$ is drawn from a Binomial distribution, $\displaystyle f_{X|P}(x|p) = {n \choose x}p^{x}(1-p)^{n-x}$.

2) The distribution of $\displaystyle p$ given the hyperparameter $\displaystyle m$ is a Beta distribution, $\displaystyle f_{P|M}(p|m) = \frac{\Gamma(m)}{\Gamma(mq)\Gamma(m(1-q))} \; p^{mq-1}(1-p)^{m(1-q)-1}$.

Here $\displaystyle q$ is a known constant between 0 and 1. Note that the shape paramters $\displaystyle \alpha = mq$ and $\displaystyle \beta=m(1-q)$ are chosen so that the mean is $\displaystyle q$.

3) $\displaystyle m$ is drawn from a uniform distribution $\displaystyle f_{M}(m)$ on $\displaystyle (\ell,\infty)$ where $\displaystyle \ell=\max(1/q,1/(1-q))$ which ensures that the distribution $\displaystyle f_{P|M}(p|m)$ is concave.

If I am understanding the Gibbs sampling procedure correctly, it goes something like this.

Start by generating a value for $\displaystyle m$ from the uniform distribution. Then:

1) Given this $\displaystyle m$ generate a value for $\displaystyle p$ from $\displaystyle f_{P|M}(p|m)$.

2) Given $\displaystyle p$ generate a value for $\displaystyle x$ from $\displaystyle f_{X|P}(x|p)$.

3) Given $\displaystyle x$ generate a new value for $\displaystyle p$ from $\displaystyle f_{P|X}(p|x)$.

4) Given this new $\displaystyle p$ generate a new value for $\displaystyle m$ from $\displaystyle f_{M|P}(m|p)$.

5) Repeat steps 1-4.

I am having trouble with step 4, since I need the conditional distribution $\displaystyle f_{M|P}(m|p)$. I should be able to get this from Bayes theorem, since $\displaystyle f_{M|P}(m|p) = \frac{f_{P|M}(p|m)f_{M}(m)}{\int f_{P|M}(p|m)f_{M}(m)\mathrm{d}m}$.

I have been unable to calculate the normalization integral that appears in the denominator above:

$\displaystyle \int_{\ell}^{\infty} \; \frac{\Gamma(m)}{\Gamma(mq)

\Gamma(m(1-q))} \;p^{mq-1}(1-p)^{m(1-q)-1}\; \mathrm{d}m$.

Question 1: Does anyone know of a method to perform this integration of the Beta distribution with respect to $\displaystyle m$ analytically? I have been unable to do so or find much on integration of the Beta distribution with respect to the shape parameter.

Question 2: For those readers who are familiar with Gibbs sampling and the Bayesian framework, please feel free to comment on the method of approach I outlined here. I am rather new at this and not entirely confident in the way I am attempting to do the sampling.

Thanks in advance,

basmati