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

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

所以也很不好...

我的问题是：

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

非常感谢您！

Answer 1

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

回复 1 和 2：有几种技术可以处理优化中的边界问题。设置为一个大的常量值往往不起作用：如果优化器超出范围，它会发现 objective 函数真的很平坦。如果确切的边界是合法的，那么将参数推到边界并添加惩罚有时会奏效。如果确切的边界是非法的，您可能会反映：例如如果 mu > 0 是必需的，有时在 objective 函数中用 abs(mu) 替换 mu 会使事情起作用。有时最好的解决办法是通过参数的变换来摆脱边界。

编辑以添加更多详细信息：

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

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

Answer 2

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

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

JN

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

Maximum-Likelihood Estimation of three parameter reverse Weibull model implementation in R

statistics

r

weibull

mle

log-likelihood