# MLE for bivariate normal data with known error variances

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.

• While this is a programming question, it "requires statistical expertise to understand or answer" and so falls within the scope, it would be more helpful for a wider variety of later readers (and more likely to attract an answer) if it focused more on the approach than the specific code implementation of that. (It's still okay to ask for code) Oct 19 '19 at 2:28
• Now to clarify the question: it's the joint distribution of all of the observations, not just the bivariate margin for the $i$th observation that will be important (without that you can't even write down a likelihood). Is it your intention that the observations be independent? Oct 19 '19 at 2:31

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