0

我正在尝试使用 ML 估计删失数据的三参数 Weibull 分布的参数。

我已经通过使用flexsurv我定义了“自己的”密度函数的包来解决这个问题。

我还按照函数文档中给出的说明flexsurv::flexsurvreg构建了包含所有必需信息的列表，以使用客户密度函数执行 MLE。

在下面你可以看到我到目前为止所做的。

library(FAdist)
library(flexsurv)

set.seed(1)
thres <- 3500
data <- rweibull(n = 1000, shape = 2.2, scale = 25000) + thres
y <- sample(c(0, 1), size = 1000, replace = TRUE)
df1 <- data.frame(x = data, status = y)

dweib3 <- function(x, shape, scale, thres, log = FALSE) {
    dweibull(x - thres, shape, scale, log = log)
}

pweib3 <- function(q, shape, scale, thres, log_p = FALSE) {
    pweibull(q - thres, shape, scale, log.p = log_p)
}

# Not required
#qweib3 <- function(p, shape, scale, thres, log.p = FALSE) {
#   if (log.p == TRUE) {
#        p <- exp(p)
#    }
#    qwei3 <- thres + qweibull(p, shape, scale)
#    return(qwei3)
#}
dweib3 <- Vectorize(dweib3)
pweib3 <- Vectorize(pweib3)

custom.weibull <- list(name = "weib3", 
    pars = c('shape', 'scale', 'thres'), location = 'scale', 
    transforms = c(log, log, log), 
    inv.transforms = c(exp, exp, exp), 
    inits = function(t) { 
        c(1.2 / sqrt((var(log(t)))), exp(mean(log(t)) + (.572 / (1.2 / sqrt((var(log(t))))))), .5 * min(t))
        }
        )
ml <- flexsurvreg(Surv(df1$x, df1$status) ~ 1, data = df1, dist = custom.weibull)

变量 y 应该代表一个单元的状态，其中 1 是失败的，0 是在审查之前未失败的单元。

形状和比例的初始值取自fitdistrplus包中定义的矩。

对于阈值参数，必须有一个约束，因为阈值必须确实小于数据的最小值。因此，阈值的约束在其最大值 .99 * t_min 将是令人满意的（这是我直到现在还没有实现的东西）。

上述 MLE 的输出如下：

> ml
Call:
flexsurvreg(formula = Surv(df1$x, df1$status) ~ 1, data = df1, 
    dist = custom.weibull)

Estimates: 
       est       L95%      U95%      se      
shape  2.37e+00  2.12e+00  2.65e+00  1.33e-01
scale  3.52e+04  3.32e+04  3.74e+04  1.08e+03
thres  2.75e+03  1.51e+03  5.02e+03  8.44e+02

N = 1000,  Events: 481,  Censored: 519
Total time at risk: 25558684
Log-likelihood = -5462.027, df = 3
AIC = 10930.05

即使有审查，估计的参数也不好。我已经用其他随机生成的数据做了几次这个过程......估计总是离“真相”很远。

因此，我需要改进我的代码或另一种可能性来估计具有 MLE 的三参数 Weibull 的参数。

4

0 回答 0