拟合3参数威布尔分布

Question

拟合3参数威布尔分布

rstatisticsdistributionmodel-fittingweibull

10

我一直在使用R进行数据分析，正在试图找出如何将我的数据拟合到三参数威布尔分布中。我已经找到了如何使用二参数威布尔分布进行拟合，但是在使用三参数时遇到了困难。

以下是我使用MASS包中的fitdistr函数拟合数据的方式：

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

x[[6]] 是我的数据的一个子集，而 y 则是我用来存储拟合结果的变量。

- Matthew Crews

也许如果您制作了一个可重现的示例，展示您的问题，人们会更容易回答。具体来说，x[[6]]看起来像什么？至少发布str(x[[6]])或更好的是dput(x[[6]])的结果。 - Andrie

2

你不能使用R中提供的内置的 weibull 分布，因为它是一个两个参数的威布尔分布。你需要计算自定义概率密度函数（3个参数），并使用它代替。 - dickoa

2个回答

1

一种替代方案：使用“ 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等替代方法。

- Marcel

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Julius Vainora · Accepted Answer

首先，您可能想查看FAdist软件包。然而，从rweibull3转换到rweibull并不难：

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

同样地，从 dweibull3 到 dweibull

> 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修改后的代码，我已经检查过几次。

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