The joint distribution for a hierarchical Bayes model Consider the following hierarchical Bayes model. We observe random variables $X_1, X_2, ..., X_n$ conditionally independent and having Poisson distributions $X_i \sim Poiss(m_i\theta_i)$, where $m_i$ are known multipliers and $\theta_i$ are unknown parameters. Assume that $\theta_1, ..., \theta_n$ are independet random variables with the same distribution $Gamma(\alpha, \lambda)$. Parameter $\lambda$ is a random variable with $Gamma(\beta, \sigma)$ distribution. Hyperparameters $\alpha, \beta, \sigma$ are known.
I would like to compute the the joint distribution.
Is the following computation correct?
$P(X_1, X_2, ..., X_n|\theta_1, \theta_2, ..., \theta_n) P(\theta_1, \theta_2, ..., \theta_n| \lambda) P(\lambda) =$
$= \prod \limits_{i=1}^n P(X_i|\theta_1, \theta_2, ..., \theta_n) \prod \limits_{i=1}^n P(\theta_i|\lambda) P(\lambda) = $
$
\prod \limits_{i=1}^n \frac{(m_i\theta_i)^{x_i}}{x_i!} e^{-m_i\theta_i} \prod \limits_{i=1}^n \frac{\lambda^{\alpha}}{\Gamma(\alpha)} \theta_i^{\alpha-1}e^{-\lambda\theta_i} \frac{\sigma^{\beta}}{\Gamma(\beta)} \lambda^{\beta-1} e^{-\sigma \lambda}
$
(Next step is to propose Gibbs sampler.)
 A: You are on the right way. Observe that
$$p(\theta_{1:n}, \lambda \vert x_{1:n}) = \prod \limits_{i=1}^n \frac{(m_i\theta_i)^{x_i}}{x_i!} e^{-m_i\theta_i} \prod \limits_{i=1}^n \frac{\lambda^{\alpha}}{\Gamma(\alpha)} \theta_i^{\alpha-1}e^{-\lambda\theta_i} \frac{\sigma^{\beta}}{\Gamma(\beta)} \lambda^{\beta-1} e^{-\sigma \lambda} \\
\propto \prod \limits_{i=1}^n (m_i\theta_i)^{x_i} e^{-m_i\theta_i} \prod \limits_{i=1}^n \lambda^{\alpha} \theta_i^{\alpha-1}e^{-\lambda\theta_i}  \lambda^{\beta-1} e^{-\sigma \lambda}.$$
To sample $\lambda | \theta_{1:n}$, it is easy to see that
$$p(\lambda \vert x_{1:n}, \theta_{1:n}) \propto \prod \limits_{i=1}^n \lambda^{\alpha} e^{-\lambda\theta_i}  \lambda^{\beta-1} e^{-\sigma \lambda} = \lambda^{\alpha+\beta-1} e^{-\lambda(\sigma + \sum_{i=1}^n \theta_i)},$$
which is a Gamma distribution with parameters $(\alpha + \beta, \sigma + \sum_{i=1}^n \theta_i)$. Note that the data does not inform us about $\lambda$. 
Conversely, to sample $\theta_{1:n} | \lambda, x_{1:n}$ we have
$$p(\theta_{1:n} | \lambda, x_{1:n}) \propto \prod \limits_{i=1}^n (m_i\theta_i)^{x_i} e^{-m_i\theta_i}  \theta_i^{\alpha-1}e^{-\lambda\theta_i} = \frac{1}{\prod_{i=1}^n m_i^{\alpha-1}} \prod \limits_{i=1}^n (m_i\theta_i)^{x_i} e^{-m_i\theta_i}  (m_i\theta_i)^{\alpha-1}e^{-\lambda\theta_i} \\
\propto  \prod \limits_{i=1}^n (m_i\theta_i)^{x_i + \alpha - 1} e^{-(m_i + \lambda)\theta_i} \\
\propto \prod \limits_{i=1}^n \theta_i^{x_i + \alpha - 1} e^{-(m_i + \lambda)\theta_i}, $$
such that each $\theta_i | \lambda, x_i$ is proportional to a Gamma distribution with parameters $(x_i + \alpha, m_i + \lambda)$. These can be sampled independently across multiple processors.
