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

编程问答 2022-05-28

问题描述

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

简而言之，我正在尝试进行最大似然估计，我需要估计 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 调用自： .local(u,...) 结束时出错：'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))

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

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

mlenonlinear-optimizationrr