Learning correlation in the bivariate Gaussian

Likelihood, MLE and the bootstrap

1 Bivariate Gaussian correlation

In this exercise, we will assume we have a random sample \(\mbox{data}=\{(x_1,y_1),\ldots,(x_n,y_n)\}\) from a bivariate normal with marginal means and variances equal to zero and one, respectively, and correlation equal to \(\rho\). For the pair \((x,y)\), the density of such bivariate distribution is given by \[ \frac{1}{2\pi\sqrt{1-\rho^2}}\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2+y^2-2\rho x y))\right\}, \] which, as a function of \(\rho\), will be denoted \(L(\rho;\mbox{data})\), for \(\rho \in (-1,1)\).

1.1 Some bivariate normal data

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

func = function(gamma){
  rho = (exp(gamma)-1)/(exp(gamma)+1)
  return(-sum(log(dnorm2(rho))))
}

set.seed(12345)
n   = 30
rho = 0.95
x   = rnorm(n)
y   = rnorm(n,rho*x,sqrt(1-rho^2))

M    = 1000
rhos = seq(0.7,1,length=M)
den  = rep(0,M)
for (i in 1:M)
  den[i] = prod(dnorm2(rhos[i]))

# MLE
run       = nlm(func,1)
gamma.mle = run$estimate
rho.mle   = (exp(gamma.mle)-1)/(exp(gamma.mle)+1)
rho.cor   = cor(x,y)

par(mfrow=c(1,2))
plot(x,y)
plot(rhos,den,type="l",xlab="rho",main="",ylab="Likelihood for rho")
abline(v=rho.mle,col=2,lwd=3)
abline(v=rho.cor,col=6,lwd=3)

2 Bootstrap

dnorm2 = function(rho){1/(2*pi*sqrt(1-rho^2))*exp(-0.5/(1-rho^2)*(x1^2+y1^2-2*rho*x1*y1))}
Rep = 10000
rho.b = rep(0,Rep)
rho.s = rep(0,Rep)
for (r in 1:Rep){
  ind = sample(1:n,size=n,replace=TRUE) 
  x1 = x[ind]
  y1 = y[ind]
  run = nlm(func,1)
  gamma.mle = run$estimate
  rho.b[r] = (exp(gamma.mle)-1)/(exp(gamma.mle)+1)
  rho.s[r] = cor(x1,y1)
}
qrho.b = quantile(rho.b,c(0.025,0.975))
qrho.s = quantile(rho.s,c(0.025,0.975))

par(mfrow=c(1,2))
hist(rho.b,breaks=seq(min(rho.b),max(rho.b),length=50),prob=TRUE,xlab="rho",main="Sampling distribution of the MLE\nvia bootstrap")
lines(density(rho.b),lwd=2,col=6)
points(qrho.b[1],0,pch=16,col=2,cex=2)
points(qrho.b[2],0,pch=16,col=2,cex=2)
points(rho.mle,0,pch=16,cex=2)
hist(rho.s,breaks=seq(min(rho.s),max(rho.s),length=50),prob=TRUE,xlab="rho",main="Sampling distribution of the sample correlation\nvia bootstrap")
lines(density(rho.s),lwd=2,col=6)
points(qrho.s[1],0,pch=16,col=2,cex=2)
points(qrho.s[2],0,pch=16,col=2,cex=2)
points(rho.cor,0,pch=16,cex=2)