error: sigma must be positive definite

Question

I am trying to create truncated multivariate normal r.vector with sigma that depends on some random vector z. Since initially sigma (in my code called nn) is not positive definite, i used function make.positive.definite() and then i got nn to be positive definite (and symmetric).

But in calling rtmvnorm i get the following error:

"Error in checkSymmetricPositiveDefinite(sigma) -sigma must be positive definite".

Any idea what might be wrong?

library(tmvtnorm)
library(matrixcalc)
library(corpcor)

zmean <- rep(0, 100)
zSigma <- diag(100)
z <- rmvnorm(n=1, mean=zmean, sigma=zSigma)

umean <- rep(0, 100)
usigma <- exp(-0.5 * rep(1, 100) + z)
nn <- t(usigma) %*% usigma

is.positive.definite(nn)

nn <- make.positive.definite(nn)
is.positive.definite(nn)
isSymmetric(nn)

a <- rep(0,100)
b <- rep(+Inf, 100)

U <- rtmvnorm(n=1, mean=umean, sigma=nn, lower=a, upper=b, algorithm="gibbs")

don't know but maybe you could try Matrix::nearPD instead, with ensureSymmetry=TRUE ? — Ben Bolker
it doesnt seem to help, assuming i am using nearPD in correct way... — nov

Karolis Koncevičius Karolis Koncevičius · Accepted Answer · 2014-11-27T22:43:05

Your problem is related to machine precision. Your matrix is positive definite after transformation, however you used different precision tolerance level in make.positive.definite than the one used in internal rtmvnorm checks.

With your example:

is.positive.definite(nn)
[1] TRUE

However:

det(nn)
[1] 0

So any function using det() to get determinant of your matrix will see it as non positive definite. You can overcome this by making sure your det(nn) returns a positive value.

One way is to add some variance in all directions:

nn <- nn + diag(ncol(nn))*0.01
det(nn)
[1] 9.98501e-253

Now rtmvnorm works.

If you want to use make.positive.definite() you have to change the tolerance:

nn <- make.positive.definite(nn, tol=1e-3)

error: sigma must be positive definite

1 Answers