2
$\begingroup$

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)
$\endgroup$
1
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.