Bivariate normal

Full conditionals

1 Bivariate normal

We say that the vector \((x,y) \in \Re^2\) follows a bivariate normal distribution with zero means, unit variances and correlation \(\rho \in (-1,1)\) if its probability density function is given by: \[ p(x,y) = \frac{1}{2\pi\sqrt{1-\rho^2}}\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2-2\rho x y + y^2)\right\}. \]

d2norm = function(x,y,rho){
  1/(2*pi*sqrt(1-rho^2))*exp(-0.5/(1-rho^2)*(x^2+y^2-2*rho*x*y))
}


rho = -0.95
N = 100
x = seq(-5,5,length=N)
y = seq(-5,5,length=N)
joint = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    joint[i,j] = d2norm(x[i],y[j],rho)
contour(x,y,joint,xlab="X",ylab="Y",drawlabels=FALSE)

2 Marginal distributions

We know that both marginal distributions are standard normal distributions: \[ x \sim N(0,1) \ \ \ \mbox{and} \ \ \ y \sim N(0,1). \]

3 Full conditional distributions

We also know that the full conditional distributions are normal as well: \[ x|y \sim N(\rho y , (1-\rho^2)) \ \ \ \mbox{and} \ \ \ y|x \sim N(\rho x , (1-\rho^2)). \]

3.1 A few steps of the Gibbs sampler

set.seed(1245)
y = 2

x = rnorm(1,rho*y,sqrt(1-rho^2))
x

## [1] -1.769124

y = rnorm(1,rho*x,sqrt(1-rho^2))
y

## [1] 1.226083

x = rnorm(1,rho*y,sqrt(1-rho^2))
x

## [1] -1.619378

y = rnorm(1,rho*x,sqrt(1-rho^2))
y

## [1] 1.221273

x = rnorm(1,rho*y,sqrt(1-rho^2))
x

## [1] -1.32667

y = rnorm(1,rho*x,sqrt(1-rho^2))
y

## [1] 1.04688

4 Gibbs sampler in action!

M     = 1000
burnin = 1000
skip  = 10
niter = burnin+M*skip
ys = rep(0,niter)
xs = rep(0,niter)
xs[1] = 0
ys[1] = 17
for (iter in 2:niter){
  xs[iter] = rnorm(1,rho*ys[iter-1],sqrt(1-rho^2))
  ys[iter] = rnorm(1,rho*xs[iter],sqrt(1-rho^2))
}

ind = seq(burnin+1,niter,by=skip)

par(mfrow=c(2,3))
ts.plot(ys[ind],xlab="Iterations",ylab="",main="Draws of X")
acf(ys[ind],main="")
hist(ys[ind],prob=TRUE,main="",xlab="X")
ts.plot(xs[ind],xlab="Iterations",ylab="",main="Draws of Y")
acf(xs[ind],main="")
hist(xs[ind],prob=TRUE,main="",xlab="Y")

4.1 Recovering the joint \(p(x,y)\)

par(mfrow=c(1,1))
plot(xs[ind],ys[ind],xlab="X",ylab="Y")
x = seq(-5,5,length=N)
y = seq(-5,5,length=N)
contour(x,y,joint,add=TRUE,col=2,drawlabels=FALSE)