[R] dmvnorm returns NaN

[David Winsemius]

On Oct 18, 2013, at 8:26 AM, Steven LeBlanc wrote:

> On Oct 17, 2013, at 11:37 PM, David Winsemius <dwinsemius at comcast.net> wrote:
> 
>> 
>> On Oct 17, 2013, at 9:11 PM, Steven LeBlanc wrote:
>> 
>>> 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()
>> 
>> There are at least 5 different version of dmvnorm. None of them are in the default packages.
>> 
>>> within the likelihood function. Particular details are provided below.
>> 
>> Complete? Except for the name of the package that has `dmvnorm`.
>> 
> 
> Package: ‘mvtnorm’ version 0.9-9992
> complete was the name of the data set.

I was clear that "complete" was the name of the dataset:

library(mvtnorm)
# First five rows of your complete:

complete <-
structure(c(0.84761637, 0.91487059, 0.84527267, 2.53821358, 1.16646209, 
3.994261, 4.952595, 4.521837, 8.37488, 6.255022), .Dim = c(5L, 
2L))

dmvnorm( complete, mean=c(1.267198, 5.475045), 
                   sigma= matrix( c( 0.6461647, 2.2289513, 2.228951,  5.697834), 2) )
Error in dmvnorm(complete, mean = c(1.267198, 5.475045), sigma = matrix(c(0.6461647,  : 
  sigma must be a symmetric matrix

So trimming the covariance elements to be exactly equal:
> matrix( c( 0.6461647, 2.2289513, 2.228951,  5.697834), 2)
          [,1]     [,2]
[1,] 0.6461647 2.228951
[2,] 2.2289513 5.697834
> dmvnorm( complete, mean=c(1.267198, 5.475045), 
                     sigma= matrix( c( 0.6461647, 2.228951, 2.228951,  5.697834), 2) )
[1] NaN NaN NaN NaN NaN
Warning message:
In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced

> dmvnorm(x=c(0,0))
[1] 0.1591549
> dmvnorm( complete, mean=c(1.267198, 5.475045) )
[1] 0.048690952 0.130494869 0.092440480 0.001059309 0.116818598
> eigen(sigma, symmetric = TRUE, only.values = TRUE)$values
[1]  6.5406882 -0.1966895

So the specified variance covariance matrix is not invertible. This can happen if you use sample statistics:

http://stats.stackexchange.com/questions/49826/what-to-do-when-sample-covariance-matrix-is-not-invertible

I'm not sure what mixtools::dvnorm is doing that avoids the problem that mvtnorm::dvnorm is identifying. Perhaps a pseudo-inverse if being constructed and use as a substitute for sigma.


> 
> Best Regards,
> Steven

David Winsemius
Alameda, CA, USA