In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) : NaNs produced

Question

I am estimating a log-likelihood using optim(). I am having some problems with the eigenvalues that don't let me find a valid hessian matrix and, therefore, the standard errors can't be calculated.

Here are the "warning" messages I get:

Warning messages:
1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
2: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
3: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
4: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced
5: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
  NaNs produced

Can somebody help me find the problem?

Here's my code:

tobit.ll <- function(theta,x,y1,y2,y3)
{
library("copula")
library(mvtnorm)
    p <- nrow(theta)
    n<-nrow(x)
    x<-as.matrix(cbind( 1,x))
    km<- as.numeric(ncol(x))
    beta1 <- theta[1:km]
      beta2 <- theta[(km+1):(2*km)]
    beta3 <- theta[(2*km+1):(km*3)]
        s1  <- theta[(3*km+1)]   ### If sigma^2 is given, take sqrt()
            if (s1 < 0)   return(1e-10)
        s2  <- theta[(3*km+2)]   ### If sigma^2 is given, take sqrt()
            if (s2 < 0)   return(1e-10)
        s3  <- theta[(3*km+3)]  
            if (s3 < 0)   return(1e-10) 
        rho12   <- theta[(3*km+4)] 
            if (rho12< -1 || rho12 > 1)   return(NA)  
        rho13   <- theta[(3*km+5)] 
            if (rho13< -1 || rho13> 1)   return(NA)  
    rho23   <- theta[(3*km+6)]  
            if (rho23< -1 || rho23> 1)   return(NA) 
    fy1 <- ifelse(y1 >0, dnorm(y1,x%*%beta1, s1), (1-pnorm((x%*%beta1)/s1)))
    fy2 <- ifelse(y2 >0, dnorm(y2,x%*%beta2, s2), (1-pnorm((x%*%beta2)/s2)))
    fy3 <- ifelse(y3 >0, dnorm(y3,x%*%beta3, s3), (1-pnorm((x%*%beta3)/s3)))
    Fy1 <- ifelse(y1>0, pnorm((y1-x%*%beta1)/s1),(1-pnorm((x%*%beta1)/s1)))
    Fy2 <- ifelse(y2>0, pnorm((y2-x%*%beta2)/s2),(1-pnorm((x%*%beta2)/s2)))
    Fy3 <- ifelse(y3>0, pnorm((y3-x%*%beta3)/s3),(1-pnorm((x%*%beta3)/s3)))
    norm.cop= ellipCopula("normal", param = c(rho12, rho13, rho23),dim = 3, dispstr = "un")
    u <- as.matrix(cbind(Fy1,Fy2,Fy3))
    dcop<- dCopula(u, norm.cop)
    fy1.lg <- ifelse(fy1==0, log(0.0001), log(fy1))
    fy2.lg <- ifelse(fy2==0, log(0.0001), log(fy2))
    fy3.lg <- ifelse(fy3==0, log(0.0001), log(fy3))
    dcop.lg <- ifelse(dcop==0, log(0.0001), log(dcop)) 
    lglik <- sum(fy1.lg+fy2.lg+fy3.lg) +sum(dcop.lg)
    return(-(lglik))
        }

Vincent Zoonekynd · Accepted Answer · 2013-04-27T19:30:49.573

0

The parameters rho12, rho13, rho23 are apparently unconstrained: there is no guarantee that the correlation matrix is positive semi-definite.

You can force it to be positive semi-definite by "truncating" the negative eigenvalues (but that makes the function you are trying to minimize non-differentiable: this can hamper convergence).

# Sample data
rho12 <- .9
rho13 <- -.9
rho23 <- .9
correlation <- matrix( c( 1,rho12,rho13, rho12,1,rho23, rho13,rho23,1 ), 3, 3 )

# Fix the correlation matrix
e <- eigen(correlation)
e$values    # Not positive semi-definite
correlation <- e$vectors %*% diag(pmax(0,e$values)) %*% t(e$vectors)
correlation <- cov2cor(correlation)

Alternatively, you could parametrize the correlation matrix in a different way (e.g., using angles).

edited Apr 27 '13 at 19:30

answered Apr 27 '13 at 18:59

Vincent Zoonekynd

31,893
5
69
78

if I add this part to my loglikelihood function the rho's wouldn't be parameters but constances, right? I don't seem to understand how to introduce this part in my function. – user2282564 Apr 27 '13 at 19:22
Also, if I just need to constraint the correlations, wouldn't change the if to constraint them make the trick? – user2282564 Apr 27 '13 at 19:22
The first 4 lines were just a numeric example: you should not include them. If would indeed be easier to return NA if the smallest eigenvalue is negative. – Vincent Zoonekynd Apr 27 '13 at 19:29
I tried constraining the rho's in the "if" on the function but gave me singularity problems. I will try your comment and see what happen. Thanks! – user2282564 Apr 27 '13 at 19:40
The constraint is not really on the individual correlations, but on the minimum eigenvalue of the correlation matrix: try `if( min( eigen( matrix( c( 1,rho12,rho13, rho12,1,rho23, rho13,rho23,1 ), 3, 3 ) )$values ) < 0 ) return(NA)`. – Vincent Zoonekynd Apr 27 '13 at 19:50
would that "if" be substituted in place of the if in the individuals correlations? – user2282564 Apr 27 '13 at 19:56
Since the fact that the correlations are between -1 and 1 is a consequence of the fact that the eigenvalues are non-negative (and that there are ones on the diagonal), you can actually replace the conditions on the individual correlations. (But they do not do any harm, so you can keep them and add the conditional on the eigenvalues). – Vincent Zoonekynd Apr 27 '13 at 20:13
Using your recommendation now give me the error: "Error in optim(par = theta0, fn = tobit.ll, method = "BFGS", x = Xi, y1 = Y1, : initial value in 'vmmin' is not finite". – user2282564 Apr 27 '13 at 20:29
This suggests that the initial values for the correlations do not yield a valid correlation matrix. This can happen, for instance, if the initial values are random. – Vincent Zoonekynd Apr 27 '13 at 20:43
Sorry to be asking so maaany questions but how could I overcome this problem? fyi, I didn't choose random number but more as cor(yi,yj). I tried other numbers that should be close based in a similar application with a similar sample and I got the same results. – user2282564 Apr 27 '13 at 20:52
What are the values of rho12, rho13, rho23 before the error? (For correlations, you can probably use 0 as initial values.) – Vincent Zoonekynd Apr 27 '13 at 21:46
The problem actually was mine, I was taking the wrong initial values. Now taking the right ones, rho12= 0.109,rho13=-0.292 and rho3= -0.419, I got a different error: "Error in optim(par = start, fn = tobit.ll, method = "BFGS", x = Xi, y1 = Y1, : non-finite finite-difference value [86]". – user2282564 Apr 27 '13 at 22:23
It probably means that the optimizer managed to evaluate the function at the initial point, but failed to find any nearby point with a finite (non-NA) value -- so it is stuck. This is likely to happen if the scale of your parameters is very different: since the correlations are between -1 and 1, the other parameters should be more or less in the same range (10 or 0.1 is fine, but 1e3 or 1e-3 may be too different). You could try to rescale some of the parameters. Or try other algorithms (`optim`'s `method` argument). – Vincent Zoonekynd Apr 27 '13 at 22:50
using other algorithms, e.g. CG, does give me a solution but I still get some negative eigen values. – user2282564 Apr 28 '13 at 18:36
I also tried scaling with "BFGS" and it result in some eigenvalues being negative too. – user2282564 Apr 28 '13 at 20:34

In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) : NaNs produced

1 Answers1