首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >威布尔分布参数估计误差

威布尔分布参数估计误差
EN

Stack Overflow用户
提问于 2019-03-27 21:05:29
回答 1查看 250关注 0票数 0

我用下面的函数来估计三参数的威布尔分布。

代码语言:javascript
复制
library(bbmle)
library(FAdist)
set.seed(16)
xl=rweibull3(50, shape = 1,scale=1, thres = 0)
dweib3l <- function(shape, scale, thres) { 
  -sum(dweibull3(xl , shape, scale, thres, log=TRUE))
}
ml <- mle2(dweib3l, start= list(shape = 1, scale = 1, thres=0), data=list(xl))

但是,当我运行上述函数时,将得到以下错误。

代码语言:javascript
复制
Error in optim(par = c(shape = 1, scale = 1, thres = 0), fn = function (p)  : 
  non-finite finite-difference value [3]
In addition: There were 16 warnings (use warnings() to see them)

有没有办法解决这个问题?谢谢!

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2019-03-27 22:00:47

问题是阈值参数是特殊的:它为分布定义了一个锐利的边界,因此任何高于数据最小值的thres值都会给出零似然值(-Inf负对数似然):如果给定的xl值小于指定的阈值,那么根据您定义的统计模型是不可能的。此外,我们已经知道阈值的最大似然值等于数据集中的最小值( MLE estimation of the bounds of a uniform distribution .

我不知道为什么解决这个问题的other questions on SO没有遇到这个特殊的问题--可能是因为他们使用的阈值的起始值远远低于数据集中的最小值……

下面,我使用一个固定的min(xl)-1e-5值作为阈值(向下移动该值可以避免数值问题,当该值正好在边界上时)。我还使用了公式接口,因此我们可以直接调用dweibull3()函数,并对形状和缩放参数设置下界(因此,我需要使用允许约束的method="L-BFGS-B" )。

代码语言:javascript
复制
ml <- mle2(xl ~ dweibull3(shape=shape, scale = scale,
                        thres=min(xl)-1e-5),
           start=list(shape=1, scale=1),
           lower=c(0,0),
           method="L-BFGS-B",
           data=data.frame(xl))

(对于简单的例子来说,公式接口是很方便的:如果您想做一些更复杂的事情,您可能想要显式地定义您自己的日志似然函数。)

如果坚持对阈值参数进行拟合,则可以设置一个上限,即(几乎)等于数据中发生的最小值--任何较大的值都会给出NA值,从而破坏优化。但是,您会发现阈值参数的估计总是收敛到这个上限.因此,这个方法实际上是以困难的方式得到了前面的答案(您还将得到关于参数在边界上的警告,以及关于不能反转Hessian参数的警告)。

代码语言:javascript
复制
eps <- 1e-8
ml3 <- mle2(xl ~ dweibull3(shape=shape, scale = scale, thres = thres),
            start=list(shape=1, scale=1, thres=-5),
            lower=c(shape=0,scale=0,thres=-Inf),
            upper=c(shape=Inf,scale=Inf,thres=min(xl)-eps),
            method="L-BFGS-B",
            data=data.frame(xl))

值得的是,如果您从一个小值开始并使用Nelder优化,那么似乎可以在不确定阈值参数的情况下对模型进行拟合:但是,它似乎提供了不可靠的结果。

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/55386447

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档