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

Question

我正在 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

所以也很糟糕……

我的问题是：

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

非常感谢您！

Answer 1

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

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

JN

Answer 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]) 应该保持在界限内。当然，进行这些更改会弄乱您的梯度公式和起始值。

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