Conducting MLE for multivariate case (bivariate normal) in R

Question

The case is that I am trying to construct an MLE algortihm for a bivariate normal case. Yet, I stuck somewhere that seems there is no error, but when I run the script it ends up with a warning.

I have a sample of size n (a fixed constant, trained with 100, but can be anything else) from a bivariate normal distribution with mean vector = (0,0) and covariance matrix = matrix(c(2.2,1.8,1.8,3),2,2)

I've tried several optimization functions (including nlm(), mle(), spg() and optim()) to maximize the likelihood function (,or minimize neg-likelihood), but there are warnings or errors.

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

I've defined the first likelihood function as follows;

bvrt_ll = function(mu,sigma,rho,sample)
{
  n = nrow(sample)
  mu_hat = c(mu[1],mu[2])
  p = length(mu)
  if(sigma[1]>0 && sigma[2]>0)
  {
    if(rho<=1 && rho>=-1)
    {
    sigma_hat = matrix(c(sigma[1]^2
                     ,sigma[1]*sigma[2]*rho
                     ,sigma[1]*sigma[2]*rho
                     ,sigma[2]^2),2,2)
    stopifnot(is.positive.definite(sigma_hat))


    neg_likelihood = (n*p/2)*log(2*pi) + (n/2)*log(det(sigma_hat)) + 0.5*sum(((sample-mu_hat)%*%solve(sigma_hat)%*%t(sample-mu_hat)))

    return(neg_likelihood)
    }
  }
  else NA

}

I prefered this one since I could set the constraints for sigmas and rho, but when I use mle()

> mle(minuslogl = bvrt_ll  ,start = list(mu = mu_est,sigma=sigma_est,rho = 
rho_est)
+     ,method = "BFGS")
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  (list) object cannot be coerced to type 'double'

I also tried nlm and spg in package BB, but they did not help as well. I tried the same function without defining constraints (inside the likelihood, not in optimization function), I could have some results but with warnings, like in nlm and spg both said the process was failed due to covariance matrix being not positive definite while it was, I think that was due to iteration, when iterating covariance matrix might not have been positive definite, and the fact that I did not define the constraints.

Thus, as a result I need to construct an mle algorithm for bivariate normal, where do I do the mistake?

NOTE: I also tried the optimization functions with the following, (I am not sure I did it correct);

neg_likelihood = function(mu,sigma,rho)
{
    if(rho>=-1 && rho<=1)
        {
          -sum(mvtnorm::dmvnorm(x=sample_10,mean=mu
                    ,sigma = matrix(c(sigma[1]^2
                    ,sigma[1]*sigma[2]*rho,sigma[1]*sigma[2]*rho
                    ,sigma[2]^2),2,2),log = T))
        }
  else NA
}

Any help is appreciated.

Thanks.

EDIT : mu is a vector of length 2 specifying the population means, sigma is a vector of length 2 (specifying population standard deviations of the random variables), and rho is a scalar as correlation coefficient between the bivariate r.v s.

Benjamin Christoffersen Benjamin Christoffersen · Accepted Answer · 2017-11-05T12:35:02

You can do it in closed form so there is no need for numeric optimization. See wiki. Just use colMeans and cov and take note of the method argument in help("cov") and this comment

The denominator n - 1 is used which gives an unbiased estimator of the (co)variance for i.i.d. observations. These functions return NA when there is only one observation (whereas S-PLUS has been returning NaN), and fail if x has length zero.

Conducting MLE for multivariate case (bivariate normal) in R

1 Answers