Gibbs sampler for a posterior distribution of hierarchical model

Question

I'm trying to get a posterior distribution of hierarchical model by using a Gibbs sampler. I haven't dealt with hierarchical models so far and I'm not really familiar with the Gibbs sampler so I'm a bit confused.

---

(1) $y_{i}|\theta_{i},\mu \stackrel{ind}\sim N(\theta_{i},\sigma^{2})$ for i =1,2,...,n ( $\sigma^{2}$ is known)

(2) $\theta_{i}|\mu \stackrel{iid}\sim N(\mu,\tau^{2})$ for i =1,2,...,n

(3) $\pi(\mu,\tau^{2})\propto \Large\frac{\tau^{2}}{\sigma^{2}+\tau^{2}}$

---

First Question

The thing I want to know is $P(\theta_{i}|y_{i})$ so I'm thinking of getting $(\theta_{(1)},\mu_{(1)},\tau^{2}_{(1)})$, $(\theta_{(2)},\mu_{(2)},\tau^{2}_{(2)})$, $(\theta_{(3)},\mu_{(3)},\tau^{2}_{(3)})$..... by using a gibbs sampler and then just pulling out only $\theta_{1},\theta_{2},\theta_{3},...$ Then does distribution of $\theta_{i}$ is P$(\theta_{i}|y_{i})$ that I'd like to know?

---

If that's the case... I actually calculated each of full conditional distributions to use a Gibbs sampler and I'm having trouble getting this one because of the leftmost term.

$P(\tau^{2}|\theta_{1},..,\theta_{n},y_{1},..,y_{n},\mu) \large\propto(\tau^{2})^{-\frac{n}{2}}exp[-\frac{\sum(\theta_{i}-\mu)^2}{2r^2}]\large\frac{\tau^{2}}{\sigma^{2}+\tau^{2}}$

I guess I could use Inverse Gamma distribution if I could just get rid of the leftmost term $\large\frac{\tau^{2}}{\sigma^{2}+\tau^{2}}$ and change a little bit.

$P(\tau^{2}|\theta_{1},..,\theta_{n},y_{1},..,y_{n},\mu) \large\propto(\tau^{2})^{-\frac{n-2}{2}-1}exp[-\frac{\sum(\theta_{i}-\mu)^2}{2r^2}] \sim IG(\frac{n-2}{2},\frac{\sum(\theta_{i}-\mu)^2}{2})$

Second Question

Since we are calculating in proportional, can I consider that term as just 1? ($\Large\frac{\tau^{2}}{\sigma^{2}+\tau^{2}} \propto 1$) so I could use Inverse Gamma distribution.

---

Can anyone help me with this? Anything is welcome. I'm kinda eager to know this:) Thank you.

Answer 1

The first question is answered in my comment, as you've seen.

For the second question, you cannot assume that $\tau^2 / (\sigma^2 + \tau^2) \propto 1$. You can use rejection sampling here.

Rejection sampling works as follows: you want to sample from a distribution $f(x)$ (in your case $x = \tau^2$) so you approximate it with another distribution $g(x)$, that you can sample from. You need to find the constant $M$ such that $f(x)/g(x) \leq M$ for all $x \geq 0$. Then,

1. Draw a sample $x^\prime$ from $g$ and a uniform $U \sim \textrm{Uniform}(0, 1)$.
2. Check if $f(x^\prime)/(Mg(x^\prime)) > U$.
3. a) If $f(x^\prime)/(Mg(x^\prime)) > U$ use $x^\prime$ as a sample from $f$ and you are finished. b) Otherwise, go back to step 1 until a sample $x^\prime$ is used as sample from $f$.

An explanation of why this process works is in the wikipedia article on rejection sampling.

In your case we can use the inverse gamma as $g$ to approximate $f(\tau^2) = P(\tau^2 \mid \theta, y, \mu)$. It's a bit more complicated here because I don't know the normalising constant for $f$. However in this case we can still do rejection sampling.

I'm going to let $a = (n - 2)/2$ and $b = \sum (\theta_i - \mu)^2 / 2$. I will let $g$ be $IG(a, b)$ to approximate $f$. Then $$f(\tau^2) = c g(\tau^2) \frac{\tau^2}{\sigma^2 + \tau^2},$$ where $c$ is an unknown constant. Then $$f(\tau^2)/g(\tau^2) = c\frac{\tau^2}{\sigma^2 + \tau^2} \leq c = M.$$

Then $$U < \frac{f(\tau^2)}{Mg(\tau^2)} \iff U < \frac{\tau^2}{\sigma^2 + \tau^2}.$$

So in conclusion, in order to sample from $P(\tau^2 \mid \theta, y, \mu)$, sample $x^\prime$ from $IG(a, b)$, and sample a uniform $U$. Accept $x^\prime$ as a draw from $P(\tau^2 \mid \theta, y, \mu)$ if $U < \frac{x^\prime}{\sigma^2 + x^\prime}.$ If you fail to accept, keep drawing pairs $(x^\prime, U)$ until you do accept.