[R] estimation of covariance matrix of a bivariate normal distribution using maximization of the log-likelihood

[Hertzog Gladys]

Dear all,
I’m new in R and I’m trying to estimate the covariance matrix of a bivariate normal distribution by maximizing the log-likelihood. The maximization really has to be performed with the non-linear minimization routine (nlm). The 2 means of the distribution are known and equal to 1.
I’ve already tried 2 different ways to compute this covariance but for each of them I obtained a lot of warnings and illogical values for the covariance matrix.
 
In the first one, I defined the bivariate normal distribution with the command dmvnorm:
 
x<-rnorm(6000, 2.4, 0.6)
x <- matrix(c(x), ncol=1)
y<-rlnorm(6000, 1.3,0.1)
y <- matrix(c(y), ncol=1)
XY <- cbind(x,y) 
 
L <- function(par,x,y) {
return (-sum(log(par[4]*dmvnorm(XY, mean=c(1,1), sigma= matrix(c(par[1], par[1]*par[2]*par[3],par[1]*par[2]*par[3], par[2] ),nrow=2, ncol=2))            )))
}
par.start<- c(0.5, 0.5 ,0.5 ,0.5)
result<-nlm(L,par.start,y=y,x=x, hessian=TRUE)
par.hat <- result$estimate
 
par.hatIl y a eu 32 avis (utilisez warnings() pour les visionner)
> par.hat <- result$estimate
> par.hat
[1] 5.149919e+01 2.520721e+02 8.734212e-03 3.996771e+02
> warnings()
Messages d'avis :
1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
2: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
3: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
4: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
5: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
6: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
7: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
8: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
9: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
10: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
11: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
12: In nlm(L, par.start, y = y, x = x, hessian = TRUE) :
  NA/Inf replaced by maximum positive value
13: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  production de NaN
 
In the second one, I wrote step by step the bivariate normal distribution in order to have each parameter separately (not in a matrix) but it didn’t work as well:
 
x<-rnorm(6000, 2.4, 0.6)
y<-rlnorm(6000, 1.3,0.1)
L <- function(par,x,y) {
return (-sum(log((1-par[4])*( (1/(2*pi*par[1]*par[2]*sqrt(1-par[3])))*exp(   (-1/2*(1-par[3]^2))* ((y-1)/par[2])^2 +((x-1)/par[1])^2 - 2*(y-1)*(x-1)/(par[2]*par[1])   )) )))
}
#par [1]= sigma_x , par [2]= sigma_y  par [3]= rho_xy  par[4] is a mixing parameter. The final step of my calculation will be to have a mixture of bivariate normal and log-normal distributions.
par.start<- c(0.5, 0.5 ,0.5 ,0.5)
result<-nlm(L,par.start,y=y,x=x, hessian=T)
par.hat <- result$estimate
par.ha
 
When I run this script, I get always 50 advices like those below:
Messages d'avis :
1: In sqrt(1 - par[3]) : production de NaN
2: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
3: In sqrt(1 - par[3]) : production de NaN
4: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
5: In sqrt(1 - par[3]) : production de NaN
6: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
7: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] *  ... : production de NaN
8: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
9: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] *  ... : production de NaN
10: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
11: In log((1 - par[4]) * ((1/(2 * pi * par[1] * par[2] *  ... : production de NaN
12: In nlm(L, par.start, y = y, x = x, hessian = T) :
  NA/Inf replaced by maximum positive value
……
Does one of you know how to use the nlm method to estimate the covariance matrix (and mixing parameter ) of a bivariate normal distribution?
 
Thank you in advance for your help and answers.
 
Best regards,
 
Gladys Hertzog
Master student in environmental engineering at ETH Zurich
-------------- next part --------------
A non-text attachment was scrubbed...
Name: message_to_mailing_list.pdf
Type: application/pdf
Size: 200565 bytes
Desc: message_to_mailing_list.pdf
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130607/3fcf8ef9/attachment.pdf>