r 中三参数 Weibull 分布的最大似然估计
Maximum Likelihood Estimation for three-parameter Weibull distribution in r
我想估计 3p Weibull 分布的尺度、形状和阈值参数。
到目前为止我所做的如下：
参考这个 post, 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])
}
到 "pre-estimate" 我的 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中，三个参数都未知下面说：
"Define the smallest failure time of ti to be tmin. Then when γ → tmin, ln(tmin - γ) → -∞. If β is less than 1, then (β - 1)ln(tmin - γ) goes to +∞ . For a given solution of β, η and γ, we can always find another set of solutions (for example, by making γ closer to tmin) that will give a larger likelihood value. Therefore, there is no MLE solution for β, η and γ."
这很有道理。出于这个原因，我想按照他们在此页面上描述的方式进行操作。
"In Weibull++, a gradient-based algorithm is used to find the MLE solution for β, η and γ. The upper bound of the range for γ is arbitrarily set to be 0.99 of tmin. Depending on the data set, either a local optimal or 0.99tmin is returned as the MLE solution for γ."
我想为 gamma 设置一个可行的区间（在我的代码中称为 'thres'），这样解决方案就在 (0, .99 * tmin) 之间。
有没有人知道如何解决这个问题？
在函数 fitdistr 中，似乎没有机会进行约束 MLE，约束一个参数。
另一种方法是通过得分向量的外积来估计渐近方差。分数向量可以从上面使用的函数 thetahat.weibul(x) 中获取。但是手动计算外积（没有函数）貌似很费时间，也没有解决约束ML估计的问题。
此致，
蒂姆
设置约束 MLE 并不难。我将在 bbmle::mle2 中执行此操作；您也可以在 stats4::mle 中执行此操作，但 bbmle 具有一些附加功能。
更大的问题是 理论上 很难定义估计值的抽样方差，当它在允许的边界时 space； Wald 方差估计背后的理论崩溃了。您仍然可以通过似然分析来计算置信区间……或者您可以 bootstrap。我运行在做这个的时候陷入了各种优化问题......我还没有真正考虑过是否有具体原因
重新格式化 three-parameter 威布尔函数 mle2 使用（以 x 作为第一个参数，以 log 作为参数）：
dweib3 <- function(x, shape, scale, thres, log=TRUE) {
dweibull(x - thres, shape, scale, log=log)
}
起始函数（稍微重新格式化）：
weib3_start <- function(x) {
mu <- mean(log(x))
sigma2 <- var(log(x))
logshape <- log(1.2 / sqrt(sigma2))
logscale <- mu + (0.572 / logshape)
logthres <- log(0.5*min(x))
list(logshape = logshape, logsc = logscale, logthres = logthres)
}
生成数据：
set.seed(1)
dat <- data.frame(x=rweibull(1000, 2.78, 750) + 450)
拟合模型：为了方便和稳定，我在对数尺度上拟合参数，但您也可以使用零边界。
tmin <- log(0.99*min(dat$x))
library(bbmle)
m1 <- mle2(x~dweib3(exp(logshape),exp(logsc),exp(logthres)),
data=dat,
upper=c(logshape=Inf,logsc=Inf,
logthres=tmin),
start=weib3_start(dat$x),
method="L-BFGS-B")
vcov(m1)，通常应该提供 variance-covariance 估计（除非估计在边界上，这里不是这种情况）给出 NaN 值...不确定为什么不进一步挖掘。
library(emdbook)
tmpf <- function(x,y) m1@minuslogl(logshape=x,
logsc=coef(m1)["logsc"],
logthres=y)
tmpf(1.1,6)
s1 <- curve3d(tmpf,
xlim=c(1,1.2),ylim=c(5.9,tmin),sys3d="image")
with(s1,contour(x,y,z,add=TRUE))
h <- lme4:::hessian(function(x) do.call(m1@minuslogl,as.list(x)),coef(m1))
vv <- solve(h)
diag(vv) ## [1] 0.002672240 0.001703674 0.004674833
(se <- sqrt(diag(vv))) ## standard errors
## [1] 0.05169371 0.04127558 0.06837275
cov2cor(vv)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8852090 -0.8778424
## [2,] 0.8852090 1.0000000 -0.9616941
## [3,] -0.8778424 -0.9616941 1.0000000
这是 log-scaled 变量的 variance-covariance 矩阵。如果要转换为原始比例的 variance-covariance 矩阵，则需要按 (x_i)*(x_j) 进行缩放（即按 t[=88= 的导数） ]信息exp(x)）。
outer(exp(coef(m1)),exp(coef(m1))) * vv
## logshape logsc logthres
## logshape 0.02312803 4.332993 -3.834145
## logsc 4.33299307 1035.966372 -888.980794
## logthres -3.83414498 -888.980794 824.831463
我不知道为什么这不适用于 numDeriv - 非常小心 上面的方差估计。 （也许太靠近理查森外推的边界？）
library(numDeriv)
hessian()
grad(function(x) do.call(m1@minuslogl,as.list(x)),coef(m1)) ## looks OK
vcov(m1)
配置文件看起来不错...（我们必须提供 std.err 因为 Hessian 矩阵不可逆）
pp <- profile(m1,std.err=c(0.01,0.01,0.01))
par(las=1,bty="l",mfcol=c(1,3))
plot(pp,show.points=TRUE)
confint(pp)
## 2.5 % 97.5 %
## logshape 0.9899645 1.193571
## logsc 6.5933070 6.755399
## logthres 5.8508827 6.134346
或者，我们可以在原始尺度上做这个......一种可能是使用log-scaling来适应，然后从这些参数开始重新适应原始比例。
wstart <- as.list(exp(unlist(weib3_start(dat$x))))
names(wstart) <- gsub("log","",names(wstart))
m2 <- mle2(x~dweib3(shape,sc,thres),
data=dat,
lower=c(shape=0.001,sc=0.001,thres=0.001),
upper=c(shape=Inf,sc=Inf,
thres=exp(tmin)),
start=wstart,
method="L-BFGS-B")
vcov(m2)
## shape sc thres
## shape 0.02312399 4.332057 -3.833264
## sc 4.33205658 1035.743511 -888.770787
## thres -3.83326390 -888.770787 824.633714
all.equal(unname(coef(m2)),unname(exp(coef(m1))),tol=1e-4)
与上述值大致相同。
如果我们更小心地限制参数，我们可以适应一个小的形状，但现在我们最终在阈值的边界上，这会给方差计算带来很多问题。
set.seed(1)
dat <- data.frame(x = rweibull(1000, .53, 365) + 100)
tmin <- log(0.99 * min(dat$x))
m1 <- mle2(x ~ dweib3(exp(logshape), exp(logsc), exp(logthres)),
lower=c(logshape=-10,logscale=0,logthres=0),
upper = c(logshape = 20, logsc = 20, logthres = tmin),
data = dat,
start = weib3_start(dat$x), method = "L-BFGS-B")
对于删失数据，您需要将dweibull替换为pweibull；有关提示，请参阅 Errors running Maximum Likelihood Estimation on a three parameter Weibull cdf。
另一种可能的解决方案是进行贝叶斯推理。在形状和尺度参数上使用比例先验，在位置参数上使用统一先验，您可以轻松地 运行 Metropolis-Hastings，如下所示。建议根据 log(shape)、log(scale) 和 log(y_min - location) 重新参数化，因为某些参数的后验变得严重偏斜，尤其是位置参数。请注意，下面的输出显示了反向转换参数的后验。
library(MCMCpack)
logposterior <- function(par,y) {
gamma <- min(y) - exp(par[3])
sum(dweibull(y-gamma,exp(par[1]),exp(par[2]),log=TRUE)) + par[3]
}
y <- rweibull(100,shape=.8,scale=10) + 1
chain0 <- MCMCmetrop1R(logposterior, rep(0,3), y=y, V=.01*diag(3))
chain <- MCMCmetrop1R(logposterior, rep(0,3), y=y, V=var(chain0))
plot(exp(chain))
summary(exp(chain))
这会产生以下输出
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.43717
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Iterations = 501:20500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
[1,] 0.81530 0.06767 0.0004785 0.001668
[2,] 10.59015 1.39636 0.0098738 0.034495
[3,] 0.04236 0.05642 0.0003990 0.001174
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
var1 0.6886083 0.768054 0.81236 0.8608 0.9498
var2 8.0756210 9.637392 10.50210 11.4631 13.5353
var3 0.0003397 0.007525 0.02221 0.0548 0.1939
我想估计 3p Weibull 分布的尺度、形状和阈值参数。
到目前为止我所做的如下：
参考这个 post, 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])
}
到 "pre-estimate" 我的 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中，三个参数都未知下面说：
"Define the smallest failure time of ti to be tmin. Then when γ → tmin, ln(tmin - γ) → -∞. If β is less than 1, then (β - 1)ln(tmin - γ) goes to +∞ . For a given solution of β, η and γ, we can always find another set of solutions (for example, by making γ closer to tmin) that will give a larger likelihood value. Therefore, there is no MLE solution for β, η and γ."
这很有道理。出于这个原因，我想按照他们在此页面上描述的方式进行操作。
"In Weibull++, a gradient-based algorithm is used to find the MLE solution for β, η and γ. The upper bound of the range for γ is arbitrarily set to be 0.99 of tmin. Depending on the data set, either a local optimal or 0.99tmin is returned as the MLE solution for γ."
我想为 gamma 设置一个可行的区间（在我的代码中称为 'thres'），这样解决方案就在 (0, .99 * tmin) 之间。
有没有人知道如何解决这个问题？
在函数 fitdistr 中，似乎没有机会进行约束 MLE，约束一个参数。
另一种方法是通过得分向量的外积来估计渐近方差。分数向量可以从上面使用的函数 thetahat.weibul(x) 中获取。但是手动计算外积（没有函数）貌似很费时间，也没有解决约束ML估计的问题。
此致， 蒂姆
设置约束 MLE 并不难。我将在 bbmle::mle2 中执行此操作；您也可以在 stats4::mle 中执行此操作，但 bbmle 具有一些附加功能。
更大的问题是 理论上 很难定义估计值的抽样方差，当它在允许的边界时 space； Wald 方差估计背后的理论崩溃了。您仍然可以通过似然分析来计算置信区间……或者您可以 bootstrap。我运行在做这个的时候陷入了各种优化问题......我还没有真正考虑过是否有具体原因
重新格式化 three-parameter 威布尔函数 mle2 使用（以 x 作为第一个参数，以 log 作为参数）：
dweib3 <- function(x, shape, scale, thres, log=TRUE) {
dweibull(x - thres, shape, scale, log=log)
}
起始函数（稍微重新格式化）：
weib3_start <- function(x) {
mu <- mean(log(x))
sigma2 <- var(log(x))
logshape <- log(1.2 / sqrt(sigma2))
logscale <- mu + (0.572 / logshape)
logthres <- log(0.5*min(x))
list(logshape = logshape, logsc = logscale, logthres = logthres)
}
生成数据：
set.seed(1)
dat <- data.frame(x=rweibull(1000, 2.78, 750) + 450)
拟合模型：为了方便和稳定，我在对数尺度上拟合参数，但您也可以使用零边界。
tmin <- log(0.99*min(dat$x))
library(bbmle)
m1 <- mle2(x~dweib3(exp(logshape),exp(logsc),exp(logthres)),
data=dat,
upper=c(logshape=Inf,logsc=Inf,
logthres=tmin),
start=weib3_start(dat$x),
method="L-BFGS-B")
vcov(m1)，通常应该提供 variance-covariance 估计（除非估计在边界上，这里不是这种情况）给出 NaN 值...不确定为什么不进一步挖掘。
library(emdbook)
tmpf <- function(x,y) m1@minuslogl(logshape=x,
logsc=coef(m1)["logsc"],
logthres=y)
tmpf(1.1,6)
s1 <- curve3d(tmpf,
xlim=c(1,1.2),ylim=c(5.9,tmin),sys3d="image")
with(s1,contour(x,y,z,add=TRUE))
h <- lme4:::hessian(function(x) do.call(m1@minuslogl,as.list(x)),coef(m1))
vv <- solve(h)
diag(vv) ## [1] 0.002672240 0.001703674 0.004674833
(se <- sqrt(diag(vv))) ## standard errors
## [1] 0.05169371 0.04127558 0.06837275
cov2cor(vv)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8852090 -0.8778424
## [2,] 0.8852090 1.0000000 -0.9616941
## [3,] -0.8778424 -0.9616941 1.0000000
这是 log-scaled 变量的 variance-covariance 矩阵。如果要转换为原始比例的 variance-covariance 矩阵，则需要按 (x_i)*(x_j) 进行缩放（即按 t[=88= 的导数） ]信息exp(x)）。
outer(exp(coef(m1)),exp(coef(m1))) * vv
## logshape logsc logthres
## logshape 0.02312803 4.332993 -3.834145
## logsc 4.33299307 1035.966372 -888.980794
## logthres -3.83414498 -888.980794 824.831463
我不知道为什么这不适用于 numDeriv - 非常小心 上面的方差估计。 （也许太靠近理查森外推的边界？）
library(numDeriv)
hessian()
grad(function(x) do.call(m1@minuslogl,as.list(x)),coef(m1)) ## looks OK
vcov(m1)
配置文件看起来不错...（我们必须提供 std.err 因为 Hessian 矩阵不可逆）
pp <- profile(m1,std.err=c(0.01,0.01,0.01))
par(las=1,bty="l",mfcol=c(1,3))
plot(pp,show.points=TRUE)
confint(pp)
## 2.5 % 97.5 %
## logshape 0.9899645 1.193571
## logsc 6.5933070 6.755399
## logthres 5.8508827 6.134346
或者，我们可以在原始尺度上做这个......一种可能是使用log-scaling来适应，然后从这些参数开始重新适应原始比例。
wstart <- as.list(exp(unlist(weib3_start(dat$x))))
names(wstart) <- gsub("log","",names(wstart))
m2 <- mle2(x~dweib3(shape,sc,thres),
data=dat,
lower=c(shape=0.001,sc=0.001,thres=0.001),
upper=c(shape=Inf,sc=Inf,
thres=exp(tmin)),
start=wstart,
method="L-BFGS-B")
vcov(m2)
## shape sc thres
## shape 0.02312399 4.332057 -3.833264
## sc 4.33205658 1035.743511 -888.770787
## thres -3.83326390 -888.770787 824.633714
all.equal(unname(coef(m2)),unname(exp(coef(m1))),tol=1e-4)
与上述值大致相同。
如果我们更小心地限制参数，我们可以适应一个小的形状，但现在我们最终在阈值的边界上，这会给方差计算带来很多问题。
set.seed(1)
dat <- data.frame(x = rweibull(1000, .53, 365) + 100)
tmin <- log(0.99 * min(dat$x))
m1 <- mle2(x ~ dweib3(exp(logshape), exp(logsc), exp(logthres)),
lower=c(logshape=-10,logscale=0,logthres=0),
upper = c(logshape = 20, logsc = 20, logthres = tmin),
data = dat,
start = weib3_start(dat$x), method = "L-BFGS-B")
对于删失数据，您需要将dweibull替换为pweibull；有关提示，请参阅 Errors running Maximum Likelihood Estimation on a three parameter Weibull cdf。
另一种可能的解决方案是进行贝叶斯推理。在形状和尺度参数上使用比例先验，在位置参数上使用统一先验，您可以轻松地 运行 Metropolis-Hastings，如下所示。建议根据 log(shape)、log(scale) 和 log(y_min - location) 重新参数化，因为某些参数的后验变得严重偏斜，尤其是位置参数。请注意，下面的输出显示了反向转换参数的后验。
library(MCMCpack)
logposterior <- function(par,y) {
gamma <- min(y) - exp(par[3])
sum(dweibull(y-gamma,exp(par[1]),exp(par[2]),log=TRUE)) + par[3]
}
y <- rweibull(100,shape=.8,scale=10) + 1
chain0 <- MCMCmetrop1R(logposterior, rep(0,3), y=y, V=.01*diag(3))
chain <- MCMCmetrop1R(logposterior, rep(0,3), y=y, V=var(chain0))
plot(exp(chain))
summary(exp(chain))
这会产生以下输出
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.43717
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Iterations = 501:20500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
[1,] 0.81530 0.06767 0.0004785 0.001668
[2,] 10.59015 1.39636 0.0098738 0.034495
[3,] 0.04236 0.05642 0.0003990 0.001174
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
var1 0.6886083 0.768054 0.81236 0.8608 0.9498
var2 8.0756210 9.637392 10.50210 11.4631 13.5353
var3 0.0003397 0.007525 0.02221 0.0548 0.1939