作业中使用MLE拟合LogLaplace分布时出现精度损失错误求助

问题场景 #

我在完成学校作业时，尝试用最大似然估计（MLE）拟合LogLaplace分布，步骤如下：

• 通过best_fit得到LogLaplace分布的初始最优拟合参数
• 定义对数似然函数loglaplace_loglikelihood
• 调用minimize函数以初始参数为起点进行优化

运行后出现错误提示：Desired error not necessarily achieved due to precision loss.，且结果中fun为nan，参数值异常。已尝试降低gtol，问题未解决。

原代码 #

#1) Find best fitting distribution parameters
best_fit(macamil2data) 
#     {'loglaplace': {'c': 1.0603278885766216,
#                     'loc': -0.04671203840594998,
#                     'scale': 10.230045114772532}}

#2) Calculate pdf using said parameters
def loglaplace_loglikelihood(params, data):
    c, loc, scale = params
    return stats.loglaplace.logpdf(data, c=c, loc=loc, scale=scale).sum()

#3) Minimize using said parameters
initial_params = [1.0603278885766216, -0.04671203840594998, 10.230045114772532]
results = minimize(loglaplace_loglikelihood, initial_params, args=(macamil2data))

print(results)

错误输出 #

message: Desired error not necessarily achieved due to precision loss.
  success: False
   status: 2
      fun: nan
        x: [ 5.127e+03 -1.765e+05 -1.945e+03]
      nit: 2
      jac: [       nan        nan        nan]
 hess_inv: [[ 6.876e-01 -1.388e+01  5.195e-02]
            [-1.388e+01  5.437e+02  5.473e+00]
            [ 5.195e-02  5.473e+00  1.000e+00]]
     nfev: 464
     njev: 116

数据集 #

macamil2data = [
0.916666, 1, 3.3, 3.38333, 3.68333, 4.16667, 4.2, 6.08333,
6.61667, 7.03333, 7.5, 7.85, 8.15, 9.08333, 9.35, 10.0833,
10.1833, 10.4333, 11.2833, 14.2, 16.5333, 20.0333, 23.8333,
30.35, 30.5167, 32.4667, 37.1, 40.8167, 45.6, 52, 70.0667,
70.5333, 85.2333, 130.967
]

问题原因 #

1. 目标函数方向错误：minimize默认寻找最小值，但MLE需要最大化对数似然。直接返回对数似然和会让优化器往错误方向迭代，导致参数发散。
2. 参数无约束：LogLaplace分布的形状参数c和尺度参数scale必须为正数。优化过程中参数越界会导致logpdf计算返回nan，最终fun变为nan。
3. 参数传递错误：minimize的args参数要求是元组，原代码中args=(macamil2data)会将数组解包为多个参数，导致函数接收的data异常。

修正后的代码 #

import numpy as np
import scipy.optimize as opt
import scipy.stats as stats

# 加载数据集
macamil2data = np.array([
0.916666, 1, 3.3, 3.38333, 3.68333, 4.16667, 4.2, 6.08333,
6.61667, 7.03333, 7.5, 7.85, 8.15, 9.08333, 9.35, 10.0833,
10.1833, 10.4333, 11.2833, 14.2, 16.5333, 20.0333, 23.8333,
30.35, 30.5167, 32.4667, 37.1, 40.8167, 45.6, 52, 70.0667,
70.5333, 85.2333, 130.967
])

# 定义负对数似然函数（转为最小化问题）
def loglaplace_neg_loglikelihood(params, data):
    c, loc, scale = params
    # 提前过滤无效参数，避免nan
    if c <= 1e-6 or scale <= 1e-6:
        return np.inf
    log_likelihood = stats.loglaplace.logpdf(data, c=c, loc=loc, scale=scale).sum()
    return -log_likelihood  # 取负转为最小化问题

initial_params = [1.0603278885766216, -0.04671203840594998, 10.230045114772532]
# 设置参数边界：c>0，scale>0，loc无约束
bounds = [(1e-6, None), (-np.inf, np.inf), (1e-6, None)]
# 修正args为元组（加逗号）
results = opt.minimize(loglaplace_neg_loglikelihood, initial_params, args=(macamil2data,), bounds=bounds)

print(results)

修正说明 #

• 目标函数调整：返回负的对数似然和，将MLE的最大化问题转为minimize可处理的最小化问题。
• 参数约束：通过bounds限制c和scale为正，同时在函数内提前检查参数，返回无穷大让优化器避开无效区域。
• 参数传递修正：args=(macamil2data,)确保数组作为单个参数传递给目标函数。

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