0

を使用して対数尤度を推定しています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))
        }
4

1 に答える 1