拟合 3 参数 Weibull 分布

Posted

技术标签:

【中文标题】拟合 3 参数 Weibull 分布【英文标题】:Fitting a 3 parameter Weibull distribution 【发布时间】:2012-08-02 19:26:51 【问题描述】:

我一直在 R 中进行一些数据分析,我试图弄清楚如何将我的数据拟合到 3 参数 Weibull 分布。我找到了如何使用 2 参数 Weibull 来完成它,但是在找到如何使用 3 参数来完成它方面做得很短。

这是我使用 MASS 包中的 fitdistr 函数拟合数据的方法:

y <- fitdistr(x[[6]], 'weibull')

x[[6]] 是我的数据的一个子集,y 是我存储拟合结果的位置。

【问题讨论】:

也许如果您创建了一个 reproducible example 来展示您的问题/问题,人们会发现它更容易回答。具体来说,x[[6]] 是什么样的。至少,发布str(x[[6]] 或者最好是dput(x[[6]]) 的结果。 您不能使用 R 中可用的内置 weibull 分布,因为它是两个参数的 weibull 分布。您必须计算自定义概率密度函数(3 个参数)并改用它。 【参考方案1】:

首先,您可能想查看FAdist package。但是,从rweibull3rweibull 并不是那么难:

> rweibull3
function (n, shape, scale = 1, thres = 0) 
thres + rweibull(n, shape, scale)
<environment: namespace:FAdist>

同样从dweibull3dweibull

> dweibull3
function (x, shape, scale = 1, thres = 0, log = FALSE) 
dweibull(x - thres, shape, scale, log)
<environment: namespace:FAdist>

所以我们有这个

> x <- rweibull3(200, shape = 3, scale = 1, thres = 100)
> fitdistr(x, function(x, shape, scale, thres) 
       dweibull(x-thres, shape, scale), list(shape = 0.1, scale = 1, thres = 0))
      shape          scale          thres    
    2.42498383     0.85074556   100.12372297 
 (  0.26380861) (  0.07235804) (  0.06020083)

编辑:如评论中所述,尝试以这种方式拟合分布时会出现各种警告

Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  non-finite finite-difference value [3]
There were 20 warnings (use warnings() to see them)
Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  L-BFGS-B needs finite values of 'fn'
In dweibull(x, shape, scale, log) : NaNs produced

对我来说,起初它只是NaNs produced,这不是我第一次看到它,所以我认为它没有那么有意义,因为估计很好。经过一番搜索,这似乎是一个非常流行的问题,我既找不到原因也找不到解决方案。一种替代方法是使用stats4 包和mle() 函数,但它似乎也有一些问题。但我可以为您提供使用 danielmedic 修改过的 code 的版本,我已经检查了几次:

thres <- 60
x <- rweibull(200, 3, 1) + thres

EPS = sqrt(.Machine$double.eps) # "epsilon" for very small numbers

llik.weibull <- function(shape, scale, thres, x)
 
  sum(dweibull(x - thres, shape, scale, log=T))


thetahat.weibull <- function(x)
 
  if(any(x <= 0)) stop("x values must be positive")

  toptim <- function(theta) -llik.weibull(theta[1], theta[2], theta[3], x)

  mu = mean(log(x))
  sigma2 = var(log(x))
  shape.guess = 1.2 / sqrt(sigma2)
  scale.guess = exp(mu + (0.572 / shape.guess))
  thres.guess = 1

  res = nlminb(c(shape.guess, scale.guess, thres.guess), toptim, lower=EPS)

  c(shape=res$par[1], scale=res$par[2], thres=res$par[3])


thetahat.weibull(x)
    shape     scale     thres 
 3.325556  1.021171 59.975470 

【讨论】:

我这样做了,我得到一个错误,并显示以下消息:fitdistr(x, function(x, shape, scale, thres) dweibull(x - thres, : 优化失败另外:警告消息:1:在 dweibull(x - thres, shape, scale) 中:产生了 NaN 2:在 dweibull(x - thres, shape, scale) 中:产生了 NaN 3:在 dweibull(x - thres, shape, scale) 中产生了 NaN 4:在 dweibull(x - thres, shape, scale) 中:产生了 NaNs @Wallhood,我编辑了答案,现在它似乎工作得很好,但不幸的是它没有提供有关方差的信息。 哇,我无法告诉你这有多棒,我有多感激。如果你在俄勒冈州的波特兰,我很乐意给你买啤酒。【参考方案2】:

另一种选择:包“lmom”。 L-矩技术估计

library(lmom)
thres <- 60
x <- rweibull(200, 3, 1) + thres
moments = samlmu(x, sort.data = TRUE)
log.moments <- samlmu( log(x), sort.data = TRUE )
weibull_3parml <- pelwei(moments)
weibull_3parml
zeta      beta     delta 
59.993075  1.015128  3.246453  

但我不知道如何在这个包或上面的解决方案中做一些拟合优度统计。其他软件包,您可以轻松地进行拟合优度统计。无论如何,您可以使用以下替代方法:ks.test 或 chisq.test

【讨论】:

以上是关于拟合 3 参数 Weibull 分布的主要内容,如果未能解决你的问题,请参考以下文章

使用 stats.exponweib.fit 在 python 中拟合 Weibull 分布

将 Weibull 累积分布拟合到 R 中的质量传递数据

使用 scipy.stats 将 Weibull 分布拟合到数据是不是表现不佳?

使用 R 中的矩法拟合 Weibull

尝试 MLE 拟合 Weibull 分布时 scipy.optimize.minimize 中的 RuntimeWarning

R语言中一组数据服从威布尔分布,怎么判断拟合的效果