时间戳

首页 > 解决方案 > R中三参数反向威布尔模型实现的最大似然估计

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

问题描述

我正在为三参数反向 Weibull 模型在 R 中实现最大似然估计，并且在获得合理结果方面遇到了一些麻烦，其中包括：糟糕的优化结果，不需要的 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 而不是 optim x我会得到一个结果，但是一个非常糟糕的结果：

 $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 函数提供函数和数据时做错了什么？

非常感谢您！

标签： rstatisticsmleweibulllog-likelihood

解决方案

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

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

推荐阅读