# Gibbs sampler from conditional distribution

I am trying to propose Gibbs sampling with the density below,

$$p(y_1,y_2,y_3)\propto \exp [-({{y}_{1}}+{{y}_{2}}+{{y}_{3}}+{{\theta}_{12}{y_1}{y_2}}+{{\theta }_{13}{y_3}{y_1}}+{{\theta }_{23}{y_2}{y_3}})]$$

where, $({{y}_{1}},{{y}_{2}},{{y}_{3}})\in R_{+}^{3}$ and ${{\theta }_{ij}}=i+j$

How do I find the full conditional distribution?

And then, I'll generate sample $\{(y_{1}^{i},y_{2}^{i},y_{3}^{i})\}$ for $i=1,...n$.

I understand Gibbs sampling, sample one variable while keep others fixed.

• Are the double subscripts in the displaymath deliberate? If so, what do they mean? Mar 7, 2012 at 0:11
• I could be mistaken, let me rewrite! Mar 7, 2012 at 0:37
• I have edited your formula as I think you mean the double exponential distribution (used by Julian Besag in several of his papers). I do not understand the $\theta_{ij}=i+j$. Otherwise, the conditionals are determined by fixing all $y_i$'s but one. You then get a standard distribution, as detailed in our Monte Carlo Statistical Methods. Mar 7, 2012 at 10:54

There you go- Gibbs Sampler: The burning period is to reach some stationarity in the samples

burning_period=5000
iterations=1000
y=matrix(nrow=(burning_period+iterations),ncol=3)
a=matrix(nrow=iterations,ncol=3)

y[1,1]=0.5        #Initial Sample
y[1,2]=0.6
y[1,3]=0.2

for(i in 2:(burning_period+iterations)){

#I have put 3,3,4 as my theta's. You may make the code generic for any choice of theta's.
t= 1+3*y[i-1,2]+4*y[i-1,3]
# Use t or t*-1 based on your requirement.
t=rexp(1,t)

y[i,1]=t
t= 1+3*y[i,1]+5*y[i-1,3]
y[i,2]=rexp(1,t)

t=1+4*y[i,1]+5*y[i,2]
y[i,3]=rexp(1,t)
}
posteriorSample=y[(burning_period+1) : (burning_period+iterations), ]

• I am working my way to there, will reply soon once I get there. thanks! Mar 8, 2012 at 18:28
• I think rexp() really does the trick, without it, generating y1, y2, y3 would take more work, maybe something like inverse transformation? just a random comment. Mar 8, 2012 at 19:38
• Yes- the theta's were for your manual tweaking :) Mar 8, 2012 at 20:21

The full conditional distributions can be found by fixing the values of the two "other" variables, then combining all the terms that you can combine and seeing what you get:

$p(y_1|y_2,y_3) \propto \exp\{-(1+\theta_{12}y_2+\theta_{13}y_3)y_1\}$

which is evidently an exponential distribution with parameter $1+\theta_{12}y_2+\theta_{13}y_3$, and similarly for $p(y_2|y_1,y_3)$, which is an exponential distribution with parameter $1 + \theta_{12}y_1 + \theta_{23}y_3$, and $p(y_3|y_1,y_2)$, which is an exponential distribution with parameter $1 + \theta_{13}y_1 + \theta_{23}y_2$.

Your Gibbs sampler will start with some values for $y_1, y_2, y_3$ and just loop over the three conditional distributions, again and again, generating exponential variates from the appropriate distribution and substituting in the current values for the the appropriate $y_i$ as it goes. Naturally you'll have to discard some initial block of variates to get rid of burn-in effects.

• Yes, you should, given that you know what they are. Mar 8, 2012 at 19:37
• I think I can answer that myself with a yes, as pointed out Praneeth's answer. Correct me if I am wrong. Mar 8, 2012 at 19:39
• Our comments crossed, I see! Yes, Praneeth's code looks good to me, except of course your thetas are 3, 4, and 5 instead of 3, 3, and 4. Mar 8, 2012 at 20:06
• Why is it ok to ignore the last term inside the exponential?(see below) I know it is a constant, but it will affect the probability in the end when you raise e to the power of that constant. $p({{y}_{1}}|{{y}_{2}},{{y}_{3}})\propto \exp \{-(1+{{\theta }_{12}}{{y}_{2}}+{{\theta }_{13}}{{y}_{3}}){{y}_{1}}+({{y}_{2}}+{{y}_{3}}+{{\theta }_{23}}{{y}_{2}}{{y}_{3}})\}$ Mar 25, 2012 at 22:13
• We can ignore the last term inside the exponential because it disappears when you apply the constant of integration to get the exact distribution. (Also, because when using the proportionality symbol, it's ok to neglect multiplicative constants.) Mar 26, 2012 at 1:35