我想估计 3p Weibull 分布的尺度、形状和阈值参数。
到目前为止，我所做的如下：
参考这篇文章，Fitting a 3 parameter Weibull distribution in R
我已经使用了这些功能
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])
}
“预先估计”我的 Weibull 参数，以便我可以将它们用作 MASS-Package 的“fitdistr”函数中参数“start”的初始值。
你可能会问我为什么要估计参数两次......原因是我需要估计的方差-协方差矩阵，它也是由 fitdistr 函数估计的。
例子：
set.seed(1)
thres <- 450
dat <- rweibull(1000, 2.78, 750) + thres
pre_mle <- thetahat.weibull(dat)
my_wb <- function(x, shape, scale, thres) {
dweibull(x - thres, shape, scale)
}
ml <- fitdistr(dat, densfun = my_wb, start = list(shape = round(pre_mle[1], digits = 0), scale = round(pre_mle[2], digits = 0),
thres = round(pre_mle[3], digits = 0)))
ml
> ml
shape scale thres
2.942548 779.997177 419.996196 ( 0.152129) ( 32.194294) ( 28.729323)
> ml$vcov
shape scale thres
shape 0.02314322 4.335239 -3.836873
scale 4.33523868 1036.472551 -889.497580
thres -3.83687258 -889.497580 825.374029
这对于形状参数大于 1 的情况非常有效。不幸的是，我的方法应该处理形状参数可能小于 1 的情况。
对于小于 1 的形状参数，这是不可能的原因如下所述：http ://www.weibull.com/hotwire/issue148/hottopics148.htm
在案例 1 中，所有三个参数都是未知的，如下所述：
“定义 ti 的最小失效时间为 tmin。那么当 γ → tmin 时，ln(tmin - γ) → -∞。如果 β 小于 1，则 (β - 1)ln(tmin - γ) 变为 + ∞ . 对于 β、η 和 γ 的给定解，我们总是可以找到另一组解（例如，通过使 γ 更接近 tmin），这将给出更大的似然值。因此，不存在 β 的 MLE 解， η 和 γ。"
这很有意义。出于这个原因，我想按照他们在此页面上描述的方式进行操作。
“在 Weibull++ 中，使用基于梯度的算法来找到 β、η 和 γ 的 MLE 解。γ 范围的上限任意设置为 tmin 的 0.99。根据数据集，局部最优或 0.99tmin 作为 γ 的 MLE 解返回。"
我想为 gamma 设置一个可行的间隔（在我的代码中称为“thres”），以便解决方案介于 (0, .99 * tmin) 之间。
有谁知道如何解决这个问题？
在函数 fitdistr 中似乎没有机会做一个受约束的 MLE，约束一个参数。
另一种方法是通过分数向量的外积来估计渐近方差。得分向量可以取自上述使用的函数 thetahat.weibul(x)。但是手动计算外积（没有函数）似乎非常耗时，并没有解决约束 ML 估计的问题。
最好的问候，蒂姆
