R中三参数反向威布尔模型实现的最大似然估计

Posted 2023-02-19

【中文标题】R中三参数反向威布尔模型实现的最大似然估计

【英文标题】：Maximum-Likelihood Estimation of three parameter reverse Weibull model implementation in R

【发布时间】：2021-11-01 03:15:53

【问题描述】：

我正在 R 中为三参数反向 Weibull 模型实施最大似然估计，但在获得合理结果时遇到了一些麻烦，其中包括： 糟糕的优化结果，不需要的 optimx 行为。除了这些，我想知道如何在这个模型中使用 parscale。

这是我的实现尝试：

为了生成数据，我使用概率积分变换：

#Generate N sigma*RWei(alph)-mu distributed points        
gen.wei <- function(N, theta) 
      alph <- theta[1]
      mu <- theta[2]
      sigma <- theta[3]
      return(
        mu - sigma * (- log (runif(N)))**(1/alph)
      )
    

现在我定义了对数似然和负对数似然以使用 optimx 优化：

#LL----
ll.wei <- function(theta,x) 
  N <- length(x)
  alph <- theta[1]
  mu <- theta[2]
  sigma <- theta[3]
  val <- sum(ifelse(
    x <= mu,
    log(alph/sigma) + (alph-1) * log( (mu-x)/sigma) - ( (mu-x)/sigma)**(alph-1),
    -Inf
  ))
  return(val)

#Negative LL----
nll.wei <- function(theta,x) 
  return(-ll.wei(theta=theta, x=x))
         

然后我定义了负 LL 的分析梯度。备注：有负LL不可微的点（上端点mu）

gradnll.wei <- function(theta,x) 
  N <- length(x)
  alph <- theta[1]
  mu <- theta[2]
  sigma <- theta[3]
  argn <- (mu-x)/sigma
  del.alph <- sum(ifelse(x <= mu,
    1/alph + log(argn) - log(argn) * argn**(alph-1),
    0
  ))
  del.mu <- sum(ifelse(x <= mu,
    (alph-1)/(mu-x) - (alph-1)/sigma * argn**(alph-2),
    0))
  del.sigma <- sum(ifelse(x <= mu,
    ((alph-1)*argn**(alph-1)-alph)/sigma,
    0))
  return (-c(del.alph, del.mu, del.sigma))

最后我尝试使用 optimx 包和方法 Nelder-Mead（无导数）和 BFGS 进行优化（我的 LL 有点平滑，只有一点，这是有问题的）。

      #MLE for Weibull
       mle.wei <- function(start,sample) 
      optimx(
        par=start,
        fn = nll.wei,
        gr = gradnll.wei,
        method = c("BFGS"),
        x = sample
      )
    
    theta.s <- c(4,1,1/2) #test for parameters
    sample <- gen.wei(100, theta.s) #generate 100 data points distributed like theta.s
mle.wei(start=c(8,4, 2), sample) #MLE Estimation

令我惊讶的是，我收到以下错误：

Error in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower,  : 
  Cannot evaluate function at initial parameters

我手动检查：nll 和 gradnll 在初始参数处都是有限的... 如果我切换到 optim 而不是 optimx 我会得到一个结果，但结果非常糟糕：

 $par
[1] 8.178674e-01 9.115766e-01 1.745724e-06

$value
[1] -1072.786

$counts
function gradient 
     574      100 

$convergence
[1] 1

$message
NULL

所以它不会收敛。如果我不向 BFGS 提供梯度，则没有结果。如果我改用 Nelder-Mead：

$par
[1] 1.026393e+00 9.649121e-01 9.865624e-18

$value
[1] -3745.039

$counts
function gradient 
     502       NA 

$convergence
[1] 1

$message
NULL

所以也很糟糕……

我的问题是：

我是否应该将支持之外的 ll 定义为 -Inf 给它一个非常高的负值（如 -1e20）来规避 -Inf 错误，还是不重要？ 与第一个类似，但用于渐变：从技术上讲，ll 不是在支持之外定义的，但由于可能性为 0，尽管在支持之外是常数，所以将 gradnll 定义为外部 0 是否明智？ 3.我检查了 evd 包中的 MLE 估计器 fgev 的实现，发现他们使用了 BFGS 方法，但没有提供梯度，即使梯度确实存在。因此我的问题是，是否存在提供梯度的情况相反，因为它没有在任何地方定义（如 my 和 evd 案例）？ 我在 optimx 中收到“参数 x 匹配多个形式参数”类型的错误，但在 optim 中没有，这让我很吃惊。我在向 optimx 函数提供函数和数据时做错了什么？

非常感谢您！

【参考方案1】：

https://web.ncf.ca/nashjc/optimx202112/ 的软件包版本至少可以处理点 args 中的一些变量冲突。

在 CRAN 上进行之前，需要进行一些单独的清理工作，但是 该软件包目前应该或多或少的健壮。

JN

【讨论】：

你是链接包的作者吗？如果是，我认为您需要披露您的隶属关系。

是的。实际上，我已经在之前的几篇文章中披露了这一点，我问那些遇到麻烦的人是否可以联系我，以便我进行修复。但是这些帖子被删除了，因为删除者声称我实际上并没有提供答案。我猜有些人不仅想要晚餐，他们还想要有人拿着喂他们的勺子。叹。并且 CRAN 风格的包在说明中有完整的归属，所以作者身份是明确的。作为信息，在我重新提交给 CRAN 之前还有一些工作需要完成，但我认为这个线程中提出的问题大部分已经解决了。

【参考方案2】：

Re 3：这是optimx 中的一种错误，但很难避免。它在计算数值梯度时使用x作为变量名；您还可以将其用作函数的“附加参数”。你可以通过重命名你的论点来解决这个问题，例如在你的函数中调用它xdata。

Re 1 & 2：有几种技术可以处理优化中的边界问题。设置一个大的常数值往往不起作用：如果优化器超出范围，它会发现目标函数非常平坦。如果确切的边界是合法的，那么将参数推到边界并添加惩罚有时会起作用。如果确切的边界是非法的，您也许可以反映：例如如果 mu > 0 是一项要求，有时在目标函数中将 mu 替换为 abs(mu) 可以使事情正常进行。有时最好的解决方案是通过转换参数来摆脱边界。

编辑添加更多细节：

对于这个问题，在我看来，参数的转换可能会起作用。我认为alpha 和sigma 都必须是肯定的。设置 alpha <- exp(theta[1]) 和 sigma <- exp(theta[3]) 将保证这一点。对mu 的限制更难，但我认为mu > max(xdata) 是必需的，所以mu <- max(xdata) + exp(theta[2]) 应该保持在界限内。当然，进行这些更改会弄乱您的梯度公式和起始值。

至于资源：恐怕我不知道。此建议基于多年的痛苦经验。

【讨论】：

非常感谢您的回答！回复 3：令人震惊！但是谢谢你告诉我。所以它实际上很容易修复。关于 1 和 2：您能否更详细一点或将我链接到您从哪里获得知识的资源？我从您的回答中得出结论，转换很有帮助。你的意思是什么转变？对 Log-L 添加惩罚可能有效，但将我的 MLE 更改为受惩罚的 MLE 不是吗？

关于您的编辑：我明白。你是对的：alpha 和 sigma 都必须是正数，并且可以通过 max(data) 从下方限制 mu，因为 mu 是分布的右上角。在第 2 点中，我问，是否存在提供梯度的情况会适得其反。你知道他们中的任何一个吗（除了我的 ll 不连续或广泛不顺畅的显而易见的一个。谢谢你的回答，它已经给了我很多见解和新想法来尝试。

提供梯度不好的一种情况是编码错误。 （我已经这样做了好几次了。）梯度的数值估计通常已经足够好了，所以我通常不会打扰它。

好吧，有趣的是我非常怀疑我是否正确地编码了渐变，因为有一条错误消息“错误：渐变函数可能是错误的 - 检查它！”来自 optimx 程序...但是我使用 numderiv 包来检查我的解析导数和数字导数的最大距离。结果是距离最多为 1e-7 阶，由此我得出结论，分析梯度是正确的......我很惊讶一个简单的优化会产生多少问题......

John Nash 在一条现已删除的评论中表示：“OP 在他的函数调用中有一个变量名冲突。¶我相信我已经在 optimx 的主要审查和修订中解决了这个问题，这是一个非常复杂的软件包。但是，由于 CRAN 需要检查，它将在很长一段时间内不会出现在 CRAN 上。感兴趣的用户可以联系 uottawa.ca 上的 nashjc 以获取新的 beta 软件包。”