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

0

参数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)

或者，您可以以不同的方式（例如，使用角度）对相关矩阵进行参数化。

于 2013-04-27T18:59:13.510 回答