R中三参数对数正态分布的拟合

Question

我想在R中拟合三个参数的对数正态分布(参考见这里 )。

我的MWE如下：

set.seed(12345)
library(FAdist)
X <- rlnorm3(n=100, shape = 2, scale = 1.5, thres = 1)

# m: Location Parameter
# s: Scale Parameter
# t: Threshold Parameter
LL3 <- function(X, m, s, t)(1/((X-t)*s*(2*pi)^0.5))*exp(((-(log(X-t)-m)^2)/(2*s^2)))

library(MASS)
fitdistr(x=X, densfun=LL3, start=list(m=2, s=1.5, t=1))

但是，此代码引发以下错误消息：

Error in stats::optim(x = c(30.9012208754183, 223.738029433835,
46.4287558537441,  :    non-finite finite-difference value [3] In addition: Warning message: In log(X - t) : NaNs produced

是否有适合三个参数分布的R包，如三参数Log-normal、Gamma、Weibull和Log-logistic distributions

Answer 1

实际上，似乎dlnorm3 (内置于FAdist包中)在x<=thres时返回的概率为零，因此将dlnorm3直接插入fitdistr似乎很好：

set.seed(12345)
library(FAdist)
library(MASS)
X <- rlnorm3(n=100, shape = 2, scale = 1.5, thres = 1)
fitdistr(X,dlnorm3,start=list(shape = 2, scale = 1.5, thres = 1))

结果：

     shape        scale        thres   
  2.31116615   1.94366899   1.02798643 
 (0.18585476) (0.23426764) (0.01480906)

如果我们使用rllog3函数生成值(我们得到了更多的极值)，这确实失败了：

Y <- rllog3(n=100, shape = 2, scale = 1.5, thres = 1)
fitdistr(Y,dlnorm3,start=list(shape = 2, scale = 1.5, thres = 1),
            method="Nelder-Mead")
## Error in stats::optim(x = c(10.1733112422871, 
##       310.508398424974, 1.08946140904075,  : 
##  non-finite finite-difference value [3]

使用debug(optim)，如果我们切换到Nelder，我们可以将问题推迟到计算Hessian。

如果我们用bbmle::mle2代替，我们至少可以得到系数(警告，Hessian不能倒.)

library(bbmle)
mle2(Y~dlnorm3(m,s,t),
     data=data.frame(Y),
     start=list(m= 2, s = 1.5, t = 1),
         method="Nelder-Mead")


## Call:
## mle2(minuslogl = Y ~ dlnorm3(m, s, t), start = list(m = 2, s = 1.5, 
##     t = 1), method = "Nelder-Mead", data = data.frame(Y))

## Coefficients:
##        m        s        t 
## 4.227529 1.606202 1.001115 

## Log-likelihood: -440.27 
## Warning message:
## In mle2(Y ~ dlnorm3(m, s, t), data = data.frame(Y), start = list(m = 2,  :
##   couldn't invert Hessian

Answer 2

误差信息表明，目标函数的梯度估计存在问题。可能发生这种情况的原因有几个，但最有可能的原因是，在分布拟合/优化过程中，其中一个参数正在变为负值或导致负值(可能阈值参数比对数正态变量更大，在这种情况下，在这些点上的分布应该是0.不幸的是，fitdistr不知道)。

解决这类问题的最好方法是尝试不同的启动参数，或者在fitdistr中找到一种方法使分布在这些情况下为0。

编辑:此外，代码还有其他错误，所以请尝试，正如Jealie所建议的：

LL3 <- function(X, m, s, t)ifelse(t>=X,0,(1/((X-t)*s*(2*pi)^0.5))*exp(((-(log(X-t)-m)^2)/(2*s^2)))) 

rlnorm3而不是rllog3