# Gibbs sampling in R

Let X a continous variable and Y a binary variable with joint distribution : $$p(x,y;\beta,\rho_1,\rho_2,\phi_1,\phi_2)=\frac{1}{Z(\beta,\rho_1,\rho_2,\phi_1,\phi_2)}\exp(-0.5 \beta x^2+1_{y=0}\rho_1 x+1_{y=1} \rho_2 x+1_{y=0} \phi_1+1_{y=1} \phi_2)$$ with Z a normalising constant.

The conditional distribution are :

$$f(x|y;\beta,\rho_1,\rho_2,\phi_1,\phi_2)\sim N(\frac{\alpha+1_{y=0}\rho_1 +1_{y=1} \rho_2}{\beta},\frac{1}{\sqrt(\beta)})$$ $$P(y=0|x;\beta,\rho_1,\rho_2,\phi_1,\phi_2)=\frac{\exp(\rho_1 x+\phi_1)}{\exp(\rho_1 x+\phi_1)+\exp(\rho_2 x+\phi_2)}$$

I tried to use Gibbs sampling to simulate from the joint distribution in R. The size of sample is 100'000, the burn in period is 1000 and every 100th is taken. The initial values are x=0 and y=0. Simulations have been performed with other initial values and it gives the same results as for x=0 and y=0. The problem is when I change the order of variables by putting y before x, it gives me different results. How is it possible? Here is the function in R :

b=2
alpha=0
rho=c(0,4)
phi=c(6,2)

sampler1=function(n,x,y){
mat=matrix(,ncol=2,nrow=n)
mat[1,]=c(x,y)
for(i in 1:n){
x=rnorm(1,(alpha+rho[y+1])/b,1/sqrt(b))
y=rbinom(1,1,exp(rho[y+1]*x+phi[y+1])/(exp(rho*x+phi)+exp(rho*x+phi)    ))
mat[i,]=c(x,y)
}
mat
}

sampler2=function(n,x,y){
mat=matrix(,ncol=2,nrow=n)
mat[1,]=c(x,y)
for(i in 1:n){
y=rbinom(1,1,exp(rho[y+1]*x+phi[y+1])/(exp(rho*x+phi)+exp(rho*x+phi)    ))
x=rnorm(1,(alpha+rho[y+1])/b,1/sqrt(b))
mat[i,]=c(x,y)
}
mat
}

par(mfrow=c(2,1))

sa_1=sampler1(100000,0,0)
sa_2=sampler2(100000,0,0)
sa_1=sa_1[seq(1000,100000,100),]
sa_2=sa_2[seq(1000,100000,100),]
plot(sa_1,main="order x y",xlab="x",ylab="y")
plot(sa_2,main="order y x",xlab="x",ylab="y")

• This isn't sufficient -- what did you derive the conditionals to be? what were the different results? What were your settings (how long was your warm-up period, for example and how long did you sample for?)? What were your starting points? Did varying the starting point make a difference? – Glen_b May 4 '16 at 14:41
• First, I would check carefully that I got the conditional distributions right. Additionally, what does different results mean? I.e. was this simply slightly different random numbers (perhaps depending on the starting point), but would have turned out to be the same distribution, if you only sampled for long enough? – Björn May 4 '16 at 14:48
• Your variables do not share the same name in the R code and in the formula. E.g., what is the connection between $\beta$, alpha and b? – Xi'an May 4 '16 at 15:23

## 1 Answer

There is a mistake in the simulation code

y=rbinom(1,1,exp(rho[y+1]*x+phi[y+1])/(exp(rho*x+phi)+exp(rho*x+phi)))

in that y cannot appear in the probability of the binomial. It should be

y=rbinom(1,1,exp(rho*x+phi)/(exp(rho*x+phi)+exp(rho*x+phi)))

And the other simulation code

x=rnorm(1,(alpha+rho[y+1])/b,1/sqrt(b))

should be

x=rnorm(1,rho[y+1]/b,1/sqrt(b))

if b means $\beta$.