r - MLE 在 R 中使用 nlminb - 理解/调试某些错误

Question

这是我在这里的第一个问题，所以我会尽量写得好。如果我犯了一个愚蠢的错误，请霸道。

简而言之，我正在尝试进行最大似然估计，我需要估计 5 个参数。我想解决的问题的一般形式如下：三个copula的加权平均值，每个copula都有一个要估计的参数，其中权重是非负的并且总和为1，也需要估计。

R 中有一些包用于对单个 copula 或具有固定权重的 copula 的加权平均值进行 MLE。但是，据我所知，没有任何软件包可以直接解决我上面概述的问题。因此，我正在尝试自己编写问题代码。有一种特殊类型的错误我无法追踪其来源。下面我试图给出一个最小的可重现示例，其中只需要估计一个参数。

library(copula)

set.seed(150)
x <- rCopula(100, claytonCopula(250)) 

# Copula density
clayton_density <- function(x, theta){
  dCopula(x, claytonCopula(theta))
}

# Negative log-likelihood function
nll.clayton <- function(theta){
  theta_trans <- -1 + exp(theta) # admissible theta values for Clayton copula
  nll <- -sum(log(clayton_density(x, theta_trans)))
  return(nll)
}

# Initial guess for optimization
guess <- function(x){
  
  init <- rep(NA, 1)
  tau.n <- cor(x[,1], x[,2], method = "kendall")
  
  # Guess using method of moments 
  itau <- iTau(claytonCopula(), tau = tau.n)
  
  # In case itau is negative, we need a conditional statement
  # Use log because it is (almost) inverse of theta transformation above
  if (itau <= 0) {
    init[1] <- log(0.1) # Ensures positive initial guess
  }
  else {
    init[1] <- log(itau)
  }
}


estimate <- nlminb(guess(x), nll.clayton)
(parameter <- -1 + exp(estimate$par)) # Retrieve estimated parameter

fitCopula(claytonCopula(), x) # Compare with fitCopula function

这在模拟具有较小 copula 参数值的数据时效果很好，并且每次都给出与 fitCopula() 几乎完全相同的答案。

对于较大的 copula 参数值，例如 250，当我使用 nlminb() 运行该行时会显示以下错误："Error in .local(u, copula, log, ...) : parameter is NA Called from: .local(u, copula, log, ...) 总结时出错：'eval' 中未实现的类型 (29)"

当我运行 fitCopula() 时，优化已完成，但会弹出此消息：“警告消息：在 dlogcdtheta(copula, u) 中：dlogcdtheta() 在第 1 列中为此显式 copula 返回了 NaN；回退到数字这些列的导数"

我已经能够使用 debug() 发现在 nlminb 的优化过程中的某个地方，感兴趣的参数被分配了值 NaN，然后​​在调用 dCopula() 时会产生此错误。但是，我不知道它发生在哪个迭代中，以及当它发生时 nlminb() 正在做什么。我怀疑也许在某些迭代中，目标函数在 Inf/-Inf 处进行评估，但我不知道 nlminb() 接下来会做什么。此外，fitCopula() 似乎也发生了类似的事情，但优化仍然进行到最后，只有上述警告。

在理解正在发生的事情、我如何自己调试它和/或我如何处理这个问题方面，我真的很感激任何帮助。从问题中可以看出，我没有很强的编码背景。非常感谢任何花时间考虑这个问题的人。

更新： 当我运行dCopula(x, claytonCopula(-1+exp(guess(x))))或等效clayton_density(x, -1+exp(guess(x)))时，很明显密度在几个数据点处评估为 0。不幸的是，通过使用创建伪x <- pobs(x)观察并不能解决问题，这可以通过重复看到dCopula(x, claytonCopula(-1+exp(guess(x))))。结果是，当应用对数函数时，我们得到了几个 -Inf 评估，这当然意味着整个负对数似然函数评估为 Inf，如运行所示nll.clayton(guess(x))。因此，除了上述查询之外，欢迎和赞赏在进行数字 MLE 时处理 log(0) 的任何提示。

第二次更新 如下编辑 nll.clayton 中的第二行似乎可以正常工作：

nll <- -sum(log(clayton_density(x, theta_trans) + 1e-8))

但是，我不知道这是否是规避问题的“好”方法，因为它不会引入潜在的大错误（尽管如果这样做会让我感到惊讶）。

标签： rnonlinear-optimizationmle