# Generating samples from Gibbs method

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

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])
}

• 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? – Christoph Hanck Apr 26 '18 at 11:49

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)$$

• Should the code for this should be something like this, then? – Martiiin Apr 26 '18 at 13:39
• 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]) } – Martiiin Apr 26 '18 at 13:39
• $(1,1)$ is not a point in the support, hence a poor choice of a starting value. – Christoph Hanck Apr 26 '18 at 13:44

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

library(rootSolve)

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=".")