MLE in R bivariate normal

Question

I'm experiencing a problem, possibly due to my coding mistake. I want to perform an MLE for a bivariate normal sample by an algorithm:

require(MASS)
require(tmvtnorm)
require(BB)
require(matrixcalc)

mu = c(0,0)
covmat = matrix(c(9,-2,-2,4),2,2)

set.seed(357)

sample = list(rmvnorm(10,mu,covmat),
          rmvnorm(100,mu,covmat),
          rmvnorm(500,mu,covmat))

I set the samples like above. I defined the negative log-likelihood as below;

neg_ll = function(mean_vec,cov_mat)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=mean_vec, sigma=cov_mat, log=TRUE))
  return(-log_ll)
}

when I run the following MLE code it returns an error;

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean_vec=xbar_vec, cov_mat=scov_mat))

Error in optim(start, f, method = method, hessian = TRUE, ...) : (list) object cannot be coerced to type 'double'

NLM function works (albeit only for the mean estimates, not for the covariance matrix). NLM returns:

nlm(f = neg_ll,xbar_vec, var(sample[[1]]),hessian = T)
$minimum
[1] 35.01874

$estimate
[1] -0.4036168  0.4703263

$gradient
[1] 1.463718e-06 3.886669e-06

$hessian
          [,1]      [,2]
[1,]  2.934318 -1.049366
[2,] -1.049366  7.769508

$code
[1] 1

$iterations
[1] 0

How do I obtain the estimates for all parameters? And what should I do to work with the MLE function?

EDIT: There was a type error on neg_ll function mean=mean is replaced by mean = mean_vec. Nonetheless, the problem is still on, nlm has estimated outputs only of mean vector.

Answer 1

If you look at the note of the optimize function, which the mle uses, it says, that the pars should be one-dimensional. That's why you are getting the error. If you rewrite your code like this:

neg_ll = function(mean1, mean2, cov11, cov12, cov21, cov22)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=c(mean1, mean2), 
                   sigma=matrix(c(cov11, cov12, cov21, cov22), 2, 2), log=TRUE))
  return(-log_ll)
}

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
             cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], 
             cov21 = scov_mat[2, 1], cov22 = scov_mat[2, 2]))

you can get past the error, but it will throw a new one, due to non-diagonalism in your covariance matrix. So, since the dmvnorm expects a diagonal matrix, you only need 1 of the upper or lower triangle, which in this case is 1 element, either cov12 or cov21.

So the code should look like this:

neg_ll = function(mean1, mean2, cov11, cov12, cov22)
{
  log_ll = sum(dmvnorm(x=sample[[1]], mean=c(mean1, mean2), 
                   sigma=matrix(c(cov11, cov12, cov12, cov22), 2, 2), log=TRUE))
  return(-log_ll)
}

xbar_vec = c(mean(sample[[1]][,1]),mean(sample[[1]][,2]))
scov_mat = var(sample[[1]])
mle(minuslogl = neg_ll, 
    start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
             cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], cov22 = scov_mat[2, 2]))

It gave me the output:

Call:
mle(minuslogl = neg_ll, start = list(mean1 = xbar_vec[1], mean2 = xbar_vec[2], 
    cov11 = scov_mat[1, 1], cov12 = scov_mat[1, 2], cov22 = scov_mat[2, 
        2]))

Coefficients:
     mean1      mean2      cov11      cov12      cov22 
-0.4036168  0.4703262  3.2228188  0.4352799  1.2171644