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

[Rolf Turner]

See in line below.

On 08/06/13 00:00, Hertzog Gladys wrote:
> 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).
Why?
> The 2 means of the distribution are known and equal to 1.
But you simulate your data using means 2.4 and 1.3.
> 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)
Why do you use matrix() and c() here? Totally unnecessary
(and confusing).
> y<-rlnorm(6000, 1.3,0.1)
> y <- matrix(c(y), ncol=1)
> XY <- cbind(x,y)
To estimate the covariance matrix underlying the data XY you could
simply do:

M <- var(XY).

If you want to impose the constraint that the true mean vector is c(1,1)
you could do:

mu <- c(1,1)
xbar <- apply(XY,2,mean)
M <- (5999/6000)*var(XY) + (xbar-mu)%*%t(xbar-mu)

It would of course make more sense in terms of your simulated data
to use:

mu <- c(2.4,1.3)

I haven't looked at your code below any further. Too much of a mess.

cheers,

Rolf Turner
>   
> 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.