MLE for bivariate normal data with known error variances

Question

I am working with a set of bivariate data arranged into columns labelled 'x' and 'y'. I also have measurements for the error variances corresponding to each observation, labelled 'sx' and 'sy'. An example of the layout is given below:

x   sx   y   sy
1   0.1  2   0.05

Assuming that the measurements for (x,y) correspond to the bivariate normal distribution

$$(\mathbf{x},\mathbf{y})^T \sim N((\mu_1,\mu_2)^T,\Sigma)$$

$$\Sigma=\begin{bmatrix} \tau_{11}+s_{x,i} & \tau_{12} \\ \tau_{12} & \tau_{22}+s_{y,i} \end{bmatrix}$$

where $s_{x,i}$ corresponds to the known error variance of the $i$th $x$ observation etc, how does one go about finding maximum likelihood estimates for the model parameters $(\mu_1, \mu_2, \tau_{11}, \tau_{12}, \tau_{22})$ in R in this case?

Can someone provide code in R which achieves this? The inclusion of the known errors in the covariance matrix is certainly complicating things.

Answer 1

Here is some code to find the maximum likelihood estimates of $\mu_1$, $\mu_2$, $\tau_{11}$, $\tau_{22}$, and $\tau_{12}=\rho \tau_{11}^{1/2} \tau_{22}^{1/2}$.

First generate some data from a bivariate normal distribution.

# Generate some data
set.seed(12345)
n <- 100  # Sample size
mu10 <- 4
mu20 <- 0
tau110 <- 1
tau220 <- 2
rho0 <- 0.7
# "Known" variances
sx <- 0.25 + 0.25*runif(n)
sy <- 0.25 + 0.25*runif(n)

library(mvtnorm)
x <- NULL
y <- NULL
for (i in 1:n) {
    # Covariance
    tau120 <- rho0 * tau110^0.5 * tau220^0.5
    # Random sample from multivariate normal
    xy <- rmvnorm(1, c(mu10, mu20), 
      matrix(c(tau110+sx[i], tau120, tau120, tau220+sy[i]), nrow=2))
    x[i] <- xy[1]
    y[i] <- xy[2]            
}

Now two functions used to calculate the maximum likelihood estimates.

logL = function(z, x, y, sx, sy) {
  mu1 <- z[1]
  mu2 <- z[2]
  tau11 <- z[3]
  tau22 <- z[4]
  rho <- z[5]

  sum((((mu2 - y)*(-(mu2*(sx + tau11)) + rho*sqrt(tau11)*sqrt(tau22)*(mu1 - x) + 
    (sx + tau11)*y))/(sx*(sy + tau22) + tau11*(sy + tau22 - rho^2*tau22)) - 
    ((mu1 - x)*((sy + tau22)*(mu1 - x) + rho*sqrt(tau11)*sqrt(tau22)*(-mu2 + y)))/
    (sx*(sy + tau22) + tau11*(sy + tau22 - rho^2*tau22)) - 
    2*log(2*pi) - log(sx*(sy + tau22) + tau11*(sy + tau22 - rho^2*tau22)))/2)
}

mle = function(x, y) {

# Get starting values
  mu1.initial <- mean(x)
  mu2.initial <- mean(y)
  tau11.initial <- var(x) - mean(sx)
  tau22.initial <- var(y) - mean(sy)
  rho.initial <- cor(x, y)
  inits <- c(mu1.initial, mu2.initial, tau11.initial, tau22.initial, rho.initial)

# Find maximum likelihood estimates
  results <- optim(inits, logL, x=x, y=y, sx=sx, sy=sy, 
    method = "L-BFGS-B", hessian=TRUE,
    lower=c(-Inf, -Inf, 0, 0, -1),
    upper=c(Inf, Inf, Inf, Inf, 1),
    control=list(fnscale=-1))

# Standard errors
  se <- diag(solve(-results$hessian))^0.5 

  list(results=results, se=se)
}

Getting the estimates:

mle(x, y)

with output