# Need help on Gibbs sampling with truncated normal and gamma

I am trying to use Gibbs Sampling to simulate a random sample from a joint distribution $f(\beta ,{{Z}_{1}},...,{{Z}_{75}},{{\lambda }_{1}},...,{{\lambda }_{75}})$, where the fully conditioned distribution function are:

$$\beta |{{Z}_{1}},...,{{Z}_{75}},{{\lambda }_{1}},...,{{\lambda }_{75}}\sim N\left( \frac{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}{{Z}_{i}}}}{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}}},\frac{1}{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}}} \right)$$

For i=1,...,24,

${{Z}_{i}}|\beta ,{{\lambda }_{1}},...,{{\lambda }_{75}}$ ~ left truncated normal at 0

$${{f}_{L}}\left( t;\beta ,0,\frac{1}{\lambda } \right)=\left\{ \begin{array}{*{35}{l}} 0 & \text{if t}\le 0 \\ \frac{{{e}^{-\frac{{{\lambda }_{i}}}{2}{{(t-\beta )}^{2}}}}}{\sqrt{2\pi /{{\lambda }_{i}}}[1-\Phi (-\beta \sqrt{{{\lambda }_{i}}})]} & \text{if t}>0 \\ \end{array} \right.$$



For i=25,...,75, ${{Z}_{i}}|\beta ,{{\lambda }_{1}},...,{{\lambda }_{75}} \sim$ right truncated normal at 0

$${{f}_{R}}\left( t;\beta ,0,\frac{1}{\lambda } \right)=\left\{ \begin{array}{*{35}{l}} \frac{{{e}^{-\frac{{{\lambda }_{i}}}{2}{{(t-\beta )}^{2}}}}}{\sqrt{2\pi /{{\lambda }_{i}}}[\Phi (-\beta \sqrt{{{\lambda }_{i}}})]} & \text{if t}\le 0 \\ 0 & \text{if t}>0 \\ \end{array} \right.$$ 

For i = 1,...,75, $${{\lambda }_{i}}|\beta ,{{Z}_{1}},...,{{Z}_{75}}\sim {\rm Gamma}\left( \frac{5}{2},\frac{2}{4+{{({{Z}_{i}}+\beta )}^{2}}} \right).$$

I am having trouble implementing this.

My algorithm: choose some fixed initial $Z_{i}$'s and $\beta$, and generate $\lambda_{i}$ from Gamma distribution.
Let $Z_{i}=0$ and $\beta=0$, I get $\lambda_{i} \sim {\rm Gamma}(5/2, 2/4)$. I then use these $\lambda_{i}$'s to generate $Z_{i}$.

Next step: how do I generate truncated normal, how do I know if my t is greater than zero or not? Any experts would like to do what I am doing right or wrong? Maybe suggestions? All welcome!

Note: I am using Java.

• You've done most of the hard work. For the truncated Normals, you can sample from them using e.g. the truncnorm package in R (among many other ways). Once you've got ways to sample all the conditionals, just write a loop to go through them one by one, storing the values and (if necessary) using them when you update the next variable. – guest Mar 4 '12 at 1:10
• he is using Java ... so for the truncated Normals, use the inverse CDF method. This is the corresponding R-Code, should be easy to adapt in Java: FA = pnorm(((a - mu)/sigma)) FB = pnorm(((b - mu)/sigma)) return(mu + sigma * qnorm(runif(length(mu)) * (FB - FA) + FA)) – joint_p Mar 4 '12 at 1:43
• yes, that works. Or more simply, generate from the untruncated Normal, and if your sample lies beyond the truncation point, try again... but your solution is better – guest Mar 4 '12 at 1:47
• @joint_p your FA and FB generates standard normal cdf. then, it return a t from f(t). what is (FB-FA)+FA doing? I understand you are trying to use it in qnorm() to get a quantile, then apply it to simga. – user1061210 Mar 4 '12 at 2:38
• It's in some paper and Rossi's book. Too lazy to look for it now, but this code is from R's bayesm package and I assure you it works. – joint_p Mar 4 '12 at 3:30

This should probably be a comment, but getting nicely formatted TeX in the comment field drives me nuts.

If your problem is generating a truncated normal distribution in java you have a number of choices.

1. Easiest, if your interval has high enough probability is just to generate normal random variables until one falls into the range. This is the essence of guest's selection.

2. Use accept-reject with some easy to simulate distribution, say a uniform over the interval in question, or an exponential if you are dealing with a tail.

3. Use the fact that the distribution function for the truncated normal is just a linear function of the distribution function for the normal. If $\Phi(x)$ is the distribution function of the normal distribution and $\phi(x)$ is the density, and you are truncating so that $a\le x \le b$ then the density function is $$\phi_{[a,b]}(x) = \frac{\phi(x)}{\Phi(b)-\Phi(a)}$$ and $$\Phi_{[a,b]}(x) = \frac{\Phi(x)-\Phi(a)}{\Phi(b)-\Phi(a)}.$$ both for $a\le x\le b$.

So if we want to simulate a truncated normal random variable, we choose a random $[0,1]$ uniform random variable $u$ and compute $$\Phi^{-1}(\Phi(A)+u(\Phi(b)-\Phi(a))$$ This is what @joint_p was doing. In R $\Phi(x)$ is pnorm(x) and $\Phi^{-1}(x)$ is qnorm(x). What these are in the java library you are using, I have no idea.

• I thought it was a very good answer, especially demonstrating how to construct phi^-1. – DWin Mar 4 '12 at 19:40
• deinst: thanks for taking the time to actually jot this down – joint_p Mar 4 '12 at 19:57
• deinst, thanks for taking the time to explain. In order for all of this $${{f}_{L}}\left( t;\beta ,0,\frac{1}{\lambda } \right)=\left\{ \begin{array}{*{35}{l}} 0 & \text{if t}\le 0 \\ \frac{{{e}^{-\frac{{{\lambda }_{i}}}{2}{{(t-\beta )}^{2}}}}}{\sqrt{2\pi /{{\lambda }_{i}}}[1-\Phi (-\beta \sqrt{{{\lambda }_{i}}})]} & \text{if t}>0 \\ \end{array} \right.$$ to become just this $$\Phi^{-1}(\Phi(a)+u(\Phi(b)-\Phi(a))$$ the a and b have to be standardized, correct? So, in my case, $a=(t-\beta)\lambda_{i}$ I am a little worried about truncated normal because I never worked with it before... – user1061210 Mar 4 '12 at 21:01
• @deinst: I posted below what I understood from your explanation. – user1061210 Mar 4 '12 at 21:03
• Yes, sort of. I am working with standard normals, it is just a linear transformation to get to any other normal. – deinst Mar 4 '12 at 21:29