В log (собственные (сигма, симметричные = TRUE, only.values ​​= TRUE) значения $): произведены NaN

Question

Я оцениваю логарифмическое правдоподобие с использованием optim(). У меня возникают некоторые проблемы с собственными значениями, которые не позволяют мне найти действительную гессианную матрицу, и поэтому стандартные ошибки не могут быть вычислены.

Вот «предупредительные» сообщения, которые я получаю:

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 

Может кто-нибудь помочь мне найти проблему?

Вот мой код:

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)) 
     } 

Answer 1

Параметры rho12, rho13, rho23, по-видимому непринужденный: нет никакой гарантии, что корреляционная матрица положительно полуопределенная.

Вы можете заставить его быть неотрицательно определенной по «усечения» отрицательные собственные значения (но это делает функцию, которую вы пытаетесь свести к минимуму недифференцируема: это может помешать сближению).

# 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) 

В качестве альтернативы, можно параметризовать матрицу корреляции другим способом (например, с использованием angles).