# 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? – guest Mar 7 '12 at 0:11
• I could be mistaken, let me rewrite! – user1061210 Mar 7 '12 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. – Xi'an Mar 7 '12 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! – user1061210 Mar 8 '12 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. – user1061210 Mar 8 '12 at 19:38
• Yes- the theta's were for your manual tweaking :) – hearse Mar 8 '12 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. – jbowman Mar 8 '12 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. – user1061210 Mar 8 '12 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. – jbowman Mar 8 '12 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}})\}$ – user1061210 Mar 25 '12 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.) – jbowman Mar 26 '12 at 1:35