[R] dmvnorm returns NaN

[Steven LeBlanc]

Greets,

I'm using nlminb() to estimate the parameters of a multivariate normal random sample with missing values and ran into an unexpected result from my call to dmvnorm() within the likelihood function. Particular details are provided below. It appears that dmvnorm() makes a call to log(eigen(sigma)). Whereas eigen(sigma) is returning a negative number, I understand log()'s complaint. However, it is a mystery to me why this data set should produce such a result.

Any suggestions?

Best Regards,
Steven

> complete
             [,1]      [,2]
 [1,]  0.84761637  3.994261
 [2,]  0.91487059  4.952595
 [3,]  0.84527267  4.521837
 [4,]  2.53821358  8.374880
 [5,]  1.16646209  6.255022
 [6,]  0.94706527  4.169510
 [7,]  0.48813564  3.349230
 [8,]  3.71828469  9.441518
 [9,]  0.08953357  1.651497
[10,]  0.68530515  5.498403
[11,]  1.52771645  8.484671
[12,]  1.55710697  5.231272
[13,]  1.89091603  4.152658
[14,]  1.08483541  5.401544
[15,]  0.58125385  5.340141
[16,]  0.24473250  2.965046
[17,]  1.59954401  8.095561
[18,]  1.57656436  5.335744
[19,]  2.73976992  8.572871
[20,]  0.87720252  6.067468
[21,]  1.18403087  3.526790
[22,] -1.03145244  1.776478
[23,]  2.88197343  7.720838
[24,]  0.60705218  4.406073
[25,]  0.58083464  3.374075
[26,]  0.87913427  5.247637
[27,]  1.10832692  3.534508
[28,]  2.92698371  8.682130
[29,]  4.04115277 11.827360
[30,] -0.57913297  1.476586
[31,]  0.84804365  7.009075
[32,]  0.79497940  3.671164
[33,]  1.58837762  5.535409
[34,]  0.63412821  3.932767
[35,]  3.14032433  9.271014
[36,] -0.18183869  1.666647
[37,]  0.57535770  6.881830
[38,]  3.21417723 10.901636
[39,]  0.29207932  4.120408
[40,]  0.65938218  5.209301
> u
[1] 1.267198 5.475045
> sigma
          [,1]     [,2]
[1,] 0.6461647 2.228951
[2,] 2.2289513 5.697834
> dmvnorm(x=complete,mean=u,sigma=sigma)
 [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
[30] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
Warning message:
In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced