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)

 A: 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...
