基于极大似然与优化法的三参数Weibull分布参数估计问题咨询

问题背景 #

使用极大似然估计（ML）结合优化方法拟合直径数据的三参数Weibull分布时，出现部分优化方法得到的参数值看似异常，但绘制的曲线与直方图拟合效果良好的矛盾现象。以下是相关代码与优化结果：

数据与自定义函数代码 #

diameters <- c(8.9,  8.8,  7.9,  8.0,  8.6,  7.0,  8.6,  8.4,  7.9,  9.3,  7.4,  8.0,  6.2,  8.4,  6.8,  8.7,  8.5,  8.1,  7.6,  7.1, 4.9,  8.4,  7.0,  7.2,  9.0,  4.2, 10.0,  7.7,  8.7,  9.2,  8.9,  8.8,  8.0,  6.6,  8.8,  6.5,  8.7,  8.3)

dweibull3 <- function(x,gamma, beta, alpha){
  (gamma/beta)* ((x-alpha)/beta)^(gamma-1)*(exp(-((x-alpha)/beta)^gamma))
}

ll.w3 <- function(p,data) sum(log(dweibull3(data, p[1],p[2],p[3])))

优化结果 #

初始值采用两参数Weibull分布的估计结果（gamma=9.019326, beta=8.381865, alpha=0），多方法优化后结果如下：

矛盾点：BFGS、Nelder-Mead的参数值明显偏离合理范围，但拟合曲线与直方图匹配；bobyqa收敛码为0，结果存疑。

问题原因分析 #

1. 代码缺少参数约束 #

三参数Weibull分布的核心约束是位置参数alpha必须小于所有观测值（仅当x > alpha时概率密度不为0），但当前代码未添加任何约束：

• 当优化算法探索到alpha为负数、gamma极大、beta极大的参数组合时，虽然参数值看似异常，但在观测数据的区间内，(x-alpha)/beta的取值会被调整到与两参数分布（alpha=0）相近的范围，导致PDF曲线拟合效果接近。
• 部分方法（如nlm、spg）出现数值溢出，是因为当alpha接近或大于数据最小值时，x-alpha为非正数，负数的分数次幂会产生复数，导致log()计算出错，最终似然值出现异常。

2. 数据特征导致似然函数存在平坦区域 #

直径数据集中在7-9之间，最小值为4.2，无明显左截断特征：

• 三参数Weibull的似然函数在gamma极大时，分布会趋近于正态分布，此时alpha的微小变化对观测区间内的PDF值影响极小，似然函数进入平坦区域。算法会认为该区域已达到收敛，但对应参数值却极端偏离合理范围。
• 这种情况下，似然函数存在多个局部极大值，部分优化方法容易陷入局部最优而非全局最优。

3. 初始值与优化方法的局限性 #

• 初始值仅采用两参数Weibull的结果（alpha=0），未覆盖三参数分布的合理参数范围，导致部分方法（如BFGS、Nelder-Mead）直接收敛到附近的局部最优。
• Nelder-Mead属于无约束优化方法，对平坦的似然函数区域鲁棒性差，容易在未达到全局最优时停止迭代；BFGS虽满足KKT条件，但找到的是局部最优解。

解决方案 #

1. 添加参数约束 #

使用支持上下界约束的优化方法（如L-BFGS-B），强制参数在合理范围内：

# 修改优化代码，添加约束：alpha < min(diameters)，gamma>0，beta>0
mle.w3_constrained <- optimx(c(gamma=9.019326, beta=8.381865, alpha=min(diameters)-0.1),
                             ll.w3,
                             method="L-BFGS-B",
                             data=diameters,
                             lower=c(1e-6, 1e-6, -Inf),
                             upper=c(Inf, Inf, min(diameters)-1e-6),
                             control=list(fnscale=-1, maxit=100))

2. 优化似然函数的数值稳定性 #

在自定义PDF函数中加入判断，避免非正数的幂运算和log(0)问题：

dweibull3 <- function(x, gamma, beta, alpha){
  ifelse(x <= alpha, 1e-300, 
         (gamma/beta) * ((x - alpha)/beta)^(gamma - 1) * exp(-((x - alpha)/beta)^gamma))
}

3. 尝试多组初始值 #

除了两参数结果，还可以用分位数法生成初始值：

• alpha设为样本最小值减去一个小量（如min(diameters)-0.5）
• 利用样本的0.16分位数和0.5分位数估计gamma和beta

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