0

打ち切りデータの ML を使用して、3 パラメーターのワイブル分布のパラメーターを推定しようとしています。

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 を使用して 3 パラメーターのワイブルのパラメーターを推定する別の可能性が必要です。

4

0 に答える 0