3p ワイブル分布のスケール、形状、およびしきい値パラメーターを推定したいと考えています。
これまでに行ったことは次のとおりです。
この投稿を参照して、R で 3 パラメーターのワイブル分布をフィッティングします。
関数を使ってみた
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])
}
MASS-Package の「fitdistr」関数の引数「start」の初期値として使用できるように、ワイブル パラメータを「事前推定」します。
パラメータを 2 回推定する理由を尋ねるかもしれません...理由は、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 では、3 つのパラメータすべてが不明であり、次のように述べられています。
「ti の最小の故障時間を tmin と定義します。このとき、γ → tmin のとき、ln(tmin - γ) → -∞。β が 1 未満の場合、(β - 1)ln(tmin - γ) は + ∞ . β, η および γ の与えられた解について, より大きな尤度値を与える別の解のセットをいつでも見つけることができます (たとえば、γ を tmin に近づけることによって). したがって、β の MLE 解はありません。 ηとγ」
これは非常に理にかなっています。まさにこの理由で、私は彼らがこのページで説明した方法でそれをやりたいと思っています.
「Weibull++ では、勾配ベースのアルゴリズムを使用して、β、η、および γ の MLE 解を見つけます。γ の範囲の上限は、tmin の 0.99 に任意に設定されます。データ セットに応じて、局所最適または 0.99tmin が γ の MLE 解として返されます。」
解が (0, .99 * tmin) の間にあるように、ガンマ (「thres」と呼ばれるコード内) の実行可能な間隔を設定したいと考えています。
この問題を解決する方法を知っている人はいますか?
関数 fitdistr では、1 つのパラメーターを制約して、制約付き MLE を実行する機会がないようです。
もう 1 つの方法は、スコア ベクトルの外積による漸近分散の推定です。スコア ベクトルは、上記で使用した関数 thetahat.weibul(x) から取得できます。しかし、外積を手動で (関数を使用せずに) 計算すると、非常に時間がかかり、制約付き ML 推定の問題は解決されません。
敬具、ティム
