I have a following homework in a subject called "Monte-Carlo Methods". I would be very thankful, if you could help me with this one, because I'm a bit stuck with this one ..

The task is as follows:

Use Gibbs method to generate uniformly distributed random vector in an ellipse $$ \frac{(x-5)^2}{25} + \frac{(y+1)^2}{10} \leq 1 $$

I understand that we should first fix $x$ and then examine how $y$ behaves (distribution-wise) etc. But I don't really understand, how this works and where should i start from. The code-part should be kind of short, but the theory beyond this is what confuses me at the moment.

I have a code example for an analog exercise from class.

n <- 100000  
x <- matrix(NA, nrow <- 3, ncol=n)  
x[,1] <- c(1, 1, 1)

for(i in 2:n) {  
x[1,i]<- rexp(1, rate = 1+x[2,i-1] * x[3,i-1])  
x[2,i]<- rexp(1, rate = 2+x[1,i] * x[3,i-1])  
x[3,i]<- rexp(1, rate = 3+x[1,i] * x[2,i])  
  • 1
    $\begingroup$ Gibbs' method requires you to know the conditional distributions. Here, conditional on some value of $x$, you can solve for the range of $y$. What is the distribution of $y$ in that range? $\endgroup$ – Christoph Hanck Apr 26 '18 at 11:49

Here is a possible (but lazy, and hence inefficient, as it lets R find the bounds of the marginal uniforms) implementation:


M <- 1e5
ff <- function(x,y) (x-5)^2/25 + (y+1)^2/10 - 1

draws <- matrix(NA,M,2)
draws[1,] <- c(5,0)

for (m in 2:M){
  interval.y <- uniroot.all(ff, c(-5,5), x=draws[m-1,1])
  y <- runif(1,interval.y[1],interval.y[2])

  interval.x <- uniroot.all(ff, c(-1,12), y=y)
  x <- runif(1,interval.x[1],interval.x[2])

  draws[m,] <- c(x,y)

plot(draws[,1], draws[,2], xlab="x", ylab="y", col="salmon", cex=0.5, pch=".")

enter image description here


The target is a density $\pi(x,y)$ proportional to the indicator $$\mathbb{I}_{\frac{(x-5)^2}{25} + \frac{(y+1)^2}{10} \leq 1}$$ hence the conditionals are also proportional to this indicator $$\pi(x|y)\propto \mathbb{I}_{\frac{(x-5)^2}{25} + \frac{(y+1)^2}{10} \leq 1}$$ and $$\pi(y|x)\propto \mathbb{I}_{\frac{(x-5)^2}{25} + \frac{(y+1)^2}{10} \leq 1}$$ which means that both conditionals are uniform: $$X|Y=y \sim {\cal U}\big(5-5\sqrt{1-.1(y+1)^2},5+5\sqrt{1-.1(y+1)^2}\big)$$ and $$Y|X=x \sim {\cal U}\big(-1-\sqrt{10(1-.04(x-5)^2)},-1+\sqrt{10(1-.04(x-5)^2)}\big)$$

  • $\begingroup$ Should the code for this should be something like this, then? $\endgroup$ – Martiiin Apr 26 '18 at 13:39
  • $\begingroup$ n <- 100000 x <- matrix(NA, nrow = 2, ncol=n) x[,1] <- c(1,1) for(i in 2:n){ x[1,i] <- runif(1, min = -1-sqrt(10*(1-0.04*(x[1,i-1] - 5)^2)), max = -1+sqrt(10*(1-0.04*(x[1,i-1] - 5)^2))) x[2,i] <- runif(1, min= 5-5*sqrt(1-0.1*(x[2,i-1]+1)^2), max = 5+5*sqrt(1-0.1*(x[2,i-1]+1)^2 )) print(x[2,i]) } $\endgroup$ – Martiiin Apr 26 '18 at 13:39
  • 1
    $\begingroup$ $(1,1)$ is not a point in the support, hence a poor choice of a starting value. $\endgroup$ – Christoph Hanck Apr 26 '18 at 13:44

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.