Text: Bayesian Data Analysis 3E by Gelman
Section: 5.3 Bayesian analysis of conjugate hierarchical models
Pages: 109-110
I am trying to follow an example about rat tumor experiments.
Let $y = (y_1, \ldots, y_J)$ and $\theta = (\theta_1, \ldots, \theta_J)$ be such that $$y_j \overset{\text{ind}}{\sim} \mathrm{Bin}(n_j, \theta_j),$$ where $n_j$ is known for all $j = 1,2,\ldots, J,$
$$\theta_j \overset{\text{ind}}{\sim} \mathrm{Beta}(\alpha, \beta),$$ and $$(\alpha, \beta) \sim \text{noninformative hyperprior}$$
The text says the following:
Given $(\alpha, \beta)$, the components of $\theta$ have independent posterior densities that are of the form $\theta_j ^A (1 - \theta_j)^B$ -- that is, beta densities -- and the joint density is $$p(\theta | \alpha, \beta, y) = \prod\limits_{j = 1}^J \frac{\Gamma(\alpha + \beta + n_j)}{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)} \theta_j^{\alpha + y_j -1} (1 - \theta_j)^{\beta + n_j - y_j - 1}$$ How is the equality obtained?
My work:
I know that $$\begin{align} p(\theta | \alpha, \beta, y) &\propto p(\theta | \alpha, \beta) p(y | \theta, \alpha, \beta) \\ &\propto \prod \limits_{j = 1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \theta_j^{\alpha - 1} (1 - \theta_j)^{\beta - 1} \prod \limits_{j = 1}^J \theta_j^{y_j} (1 - \theta)^{n_j - y_j} \\ &= \prod \limits_{j = 1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \theta_j^{\alpha + y_j - 1} (1 - \theta_j)^{\beta +n_j - y_j - 1} \end{align}$$
I almost have the product of independent $\mathrm{Beta}(\alpha + y_j, \beta + n_j - y_j)$ random variables. Now I find the proper normalizing "constant".
$$\begin{align} &... \prod \limits_{j = 1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \theta_j^{\alpha + y_j - 1} (1 - \theta_j)^{\beta +n_j - y_j - 1} \\ &= \prod \limits_{j=1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \frac{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)}{\Gamma(\alpha + \beta + n_j)} \frac{\Gamma(\alpha + \beta + n_j)}{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)} \theta_j^{\alpha + y_j - 1} (1 - \theta_j)^{\beta +n_j - y_j - 1} \\ &= \prod \limits_{j=1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \frac{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)}{\Gamma(\alpha + \beta + n_j)} \prod \limits_{j=1}^J \frac{\Gamma(\alpha + \beta + n_j)}{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)} \theta_j^{\alpha + y_j - 1} (1 - \theta_j)^{\beta +n_j - y_j - 1} \end{align}$$
So now I need to show that $\prod \limits_{j=1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} \frac{\Gamma(\alpha + y_j) \Gamma(\beta + n_j - y_j)}{\Gamma(\alpha + \beta + n_j)} = 1$, but I don't believe this is true. Actually, I've tried a few cases where I fix $J = 1$ and use some easy values for $\alpha, \beta, y_j, n_j$ and I don't get that it equals one.
Can someone offer another approach?
You are there as soon as you have $$ p(\theta|\alpha,\beta,y)\propto\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1}, $$ in fact can just throw away the constant that seems to be confusing matters and write $$ p(\theta|\alpha,\beta,y)\propto\prod_{j=1}^{J}\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1}, $$ to find the normalising constant you have $$ \begin{align} \int\cdots\int\prod_{j=1}^J\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1}\operatorname{d}\theta_1\cdots\operatorname{d}\theta_J&=\prod_{j=1}^{J}\int\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1}\operatorname{d}\theta_j \\ &=\prod_{j=1}^{J} B(\alpha+y_j, \beta +n_j-y_j) \end{align} $$ and then just put it all together $$ \begin{align} p(\theta|\alpha,\beta,y)&=\frac{1}{\prod_{j=1}^{J}B(\alpha+y_j, \beta + n_j - y_j)}\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1} \\ &=\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta+n_j)}{\Gamma(\alpha+y_j)\Gamma(\beta + n_j - y_j)}\theta_j^{\alpha+y_j-1}(1-\theta_j)^{\beta+n_j-y_j-1}. \end{align} $$