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

Question

我想在 R 中拟合三参数对数正态分布（参见 here 以供参考）。

我的 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-Mead，我们可以将问题推迟到计算 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