2023年fitdistrplus拟合截断对数正态分布的失效问题及解决办法

问题背景 #

尝试为截断数据拟合对数正态分布，参考旧方案后发现失效：truncdist包的dtrunc和ptrunc函数无法通过fitdistrplus的测试函数。

原代码 #

dtruncated_log_normal <- function(x, a,b, meanlog, sdlog)
  dtrunc(x, "lnorm", a=a, b=b, meanlog=meanlog, sdlog=sdlog)

ptruncated_log_normal <- function(q, a,b, meanlog, sdlog)
  ptrunc(q, "lnorm", a=a, b=b, meanlog=meanlog, sdlog=sdlog)

fit <- fitdist(s, "truncated_log_normal", start=list(a=0.001, b=90, meanlog=mean(log(s)), sdlog=sd(log(s))))

错误信息 #

Error in fitdist(s, "truncated_log_normal", start = list(a = 0.001, b = 90, :
the function mle failed to estimate the parameters,
with the error code 100
In addition: Warning messages:
1: In fitdist(s, "truncated_log_normal", start = list(a = 0.001, b = 90, :
The dtruncated_log_normal function should return a vector of with NaN values when input has inconsistent values and not raise an error
2: In fitdist(s, "truncated_log_normal", start = list(a = 0.001, b = 90, :
The ptruncated_log_normal function should return a vector of with NaN values when input has inconsistent parameters and not raise an error

数据样本 #

数据已清理，无异常值、空值或NA，样本如下：

> dput(head(s,150))
c(88.443, 89.296, 89.327, 87.776, 89.405, 89.824, 89.997, 87.678, 
89.665, 88.814, 88.841, 89.728, 89.365, 89.476, 89.189, 88.251, 
88.939, 89.945, 89.567, 89.613, 89.317, 89.622, 87.674, 89.19, 
89.782, 89.891, 89.954, 89.556, 89.093, 89.637, 89.052, 87.395, 
87.835, 89.357, 87.733, 89.459, 88.197, 88.539, 88.564, 87.857, 
88.74, 88.955, 89.691, 88.102, 89.635, 89.116, 89.584, 88.288, 
86.95, 89.182, 89.435, 88.93, 87.567, 89.083, 88.52, 88.897, 
89.54, 88.557, 89.269, 89.854, 89.31, 88.274, 89.126, 89.431, 
88.257, 88.872, 88.978, 89.03, 87.434, 88.305, 89.656, 87.556, 
89.209, 89.508, 87.781, 88.068, 89.933, 87.256, 88.906, 89.067, 
88.92, 87.947, 88.196, 88.951, 89.594, 88.378, 87.482, 88.817, 
89.65, 89.392, 89.932, 87.896, 89.909, 89.265, 89.954, 89.827, 
87.49, 87.786, 89.208, 89.728, 88.905, 87.566, 86.612, 88.363, 
87.457, 87.639, 88.907, 88.425, 87.244, 88.546, 88.221, 89.293, 
87.469, 87.31, 89.107, 88.442, 89.133, 88.812, 88.418, 89.456, 
88.512, 89.514, 87.446, 88.374, 89.282, 87.415, 89.004, 87.627, 
89.107, 89.168, 89.589, 89.288, 88.496, 89.807, 87.518, 88.796, 
88.001, 87.322, 87.353, 88.055, 88.81, 88.456, 87.876, 87.7, 
88.675, 88.996, 89.479, 86.781, 86.928, 87.356)

2023年可行解决办法 #

方案1：自定义稳健版密度与分布函数 #

核心问题是旧的自定义函数在参数/输入异常时直接报错，而非返回NaN。修改函数加入错误捕获机制：

library(fitdistrplus)
library(truncdist)

dtruncated_log_normal <- function(x, a, b, meanlog, sdlog) {
  tryCatch(
    dtrunc(x, "lnorm", a = a, b = b, meanlog = meanlog, sdlog = sdlog),
    error = function(e) {
      rep(NaN, length(x))
    }
  )
}

ptruncated_log_normal <- function(q, a, b, meanlog, sdlog) {
  tryCatch(
    ptrunc(q, "lnorm", a = a, b = b, meanlog = meanlog, sdlog = sdlog),
    error = function(e) {
      rep(NaN, length(q))
    }
  )
}

若截断边界a/b是已知固定值，不要将其放入待估参数，固定后仅拟合meanlog和sdlog能大幅提升稳定性：

fixed_a <- min(s)
fixed_b <- 90

fit <- fitdist(s, "truncated_log_normal", 
               start = list(meanlog = mean(log(s)), sdlog = sd(log(s))),
               fix.arg = list(a = fixed_a, b = fixed_b))

方案2：改用bbmle包直接最大化似然函数 #

跳过fitdistrplus的函数检测，手动构造对数似然函数进行拟合：

library(bbmle)

log_lik_trunc_lnorm <- function(meanlog, sdlog, a, b) {
  sum(dtrunc(s, "lnorm", a = a, b = b, meanlog = meanlog, sdlog = sdlog, log = TRUE))
}

# 若a/b固定，可从参数列表中移除
fit_mle <- mle2(log_lik_trunc_lnorm,
                start = list(meanlog = mean(log(s)), sdlog = sd(log(s)), a = 0.001, b = 90),
                method = "L-BFGS-B")

summary(fit_mle)

方案3：数据转换后拟合标准截断正态分布 #

观察样本数据为上截断（所有值≤90），可将数据转换为对数形式后拟合截断正态分布，再转换回对数正态参数：

library(truncnorm)

y <- log(s)
fit_truncnorm <- fitdistr(y, "truncnorm", 
                          start = list(a = -Inf, b = log(90), mean = mean(y), sd = sd(y)))

# 转换回对数正态分布参数
lnorm_params <- list(meanlog = fit_truncnorm$estimate["mean"],
                     sdlog = fit_truncnorm$estimate["sd"])

内容的提问来源于stack exchange，提问作者Sol1