# Gibbs sampler for conditionals that are exponential: Example from Casella & George paper

I am trying to work out Example 2 from Casella and George's paper "Explaining the Gibbs Sampler" in R.

The example is:

f(x | y) = y*exp(-yx)
f(y | x) = x*exp(-xy)

where x, y are defined on [0,B]


The idea is to generate the distribution of f(x).

My code (below) works, but what I do not understand is the role of the interval on which x and y are defined. In the paper this interval is [0, B]. In their example, the authors say that their simulation uses B = 5, but I am not clear where to incorporate this into my code...

X = rep(0, 500)
Y = rep(0, 500)

k = 15

for (i in 1:500) {
x = rep(1, k)
y = rep(1, k)

for (j in 2:k) {
x[j] = y[j-1]*rexp(1, y[j-1])
y[j] = x[j]*rexp(1, x[j])
}

X[i] = x[k]
Y[i] = y[k]
}
print(max(X))
print(max(Y))

hist(X, breaks=40, freq=F)


I think this works:

X = rep(0, 500)
Y = rep(0, 500)

k = 15

for (i in 1:500) {
x = rep(1, k)
y = rep(1, k)

for (j in 2:k) {
temp_x = 6
while(temp_x>5) {
x[j] = rexp(1,y[j-1])
temp_x = x[j]
}

temp_y = 6
while(temp_y>5) {
y[j] = rexp(1,x[j])
temp_y = y[j]
}
}

X[i] = x[k]
Y[i] = y[k]
}
print(max(X))
print(max(Y))

hist(X, breaks=40, freq=F)


The bits I've changed are:

• The conditional samplers for $x$ and $y$, which weren't actually drawing from exponential distributions (they should be, the density for an $Exp(\lambda)$ random variable is $\lambda e^{-\lambda x}$ not just $e^{-\lambda x}$, so you don't need to multiply your samples by $\lambda$ again after drawing them from rexp)
• I've used a rejection sampler (very simple here) to draw samples from $x$ and $y$ conditional on them being less than 5. You should be able to see what I've done from the code. The resulting marginal histogram I get for $x$ looks like the one in the Casella paper you mention, so I think it works ok! They mention in the paper that you have to do this in order for $x$ to have a marginal distribution.

Hope that helps...

• Yes, I understood the first correction... But why while(temp_x>5)? Should that not be a <5?
– Ravi
Jan 3 '14 at 11:05
• No - you keep drawing samples until you get one less than 5, so if temp_x<5 you break the loop... Jan 3 '14 at 11:18