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.

  • $\begingroup$ 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) $\endgroup$
    – Glen_b
    Oct 19, 2019 at 2:28
  • $\begingroup$ 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? $\endgroup$
    – Glen_b
    Oct 19, 2019 at 2:31

1 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
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)

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),

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

  list(results=results, se=se)

Getting the estimates:

mle(x, y)

with output

Maximum likelihood output


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.