使用scipy实现MLE估计时如何避免数值溢出问题

你当前的代码虽然采用了对数似然的思路，但没有对对数表达式做代数化简，仍然先计算5e5 ** alpha、x ** (alpha+1)这类指数项再取对数，当alpha取值稍大时，指数运算会直接触发数值溢出，导致结果为NaN，等于没有发挥对数优化避免溢出的作用。另外帕累托分布的形状参数alpha必须大于0，优化时如果迭代到alpha<=0的取值，也会触发np.log(alpha)返回NaN。

先对对数PDF做代数展开：
原始帕累托PDF为 f(x) = alpha * 500000^alpha / x^(alpha+1)
取对数后展开可得：
log(f(x)) = log(alpha) + alpha * log(500000) - (alpha + 1) * log(x)
展开后的表达式完全没有高次指数运算，仅涉及乘法、加法，从根源避免溢出。

import numpy as np
from scipy.optimize import minimize

# 导入数据集
data = np.array([1.000, 1.000, 1.000, 1.004, 1.005, 1.008, 1.014, 1.015, 1.023, 1.035, 1.038, 1.046, 1.048, 1.050, 1.050, 1.052, 1.052, 1.057, 1.063, 1.070, 1.070, 1.076, 1.087, 1.090, 1.091, 1.096, 1.101, 1.102, 1.113, 1.114, 1.120, 1.130, 1.131, 1.150, 1.152, 1.154, 1.155, 1.162, 1.170, 1.177, 1.189, 1.191, 1.193, 1.200, 1.200, 1.200, 1.200, 1.205, 1.210, 1.218, 1.238, 1.238, 1.241, 1.250, 1.250, 1.256, 1.257, 1.272, 1.278, 1.289, 1.299, 1.300, 1.316, 1.331, 1.349, 1.374, 1.378, 1.382, 1.396, 1.426, 1.429, 1.439, 1.443, 1.446, 1.473, 1.475, 1.478, 1.499, 1.506, 1.559, 1.568, 1.594, 1.609, 1.626, 1.649, 1.650, 1.669, 1.675, 1.687, 1.715, 1.720, 1.735, 1.750, 1.755, 1.787, 1.797, 1.805, 1.898, 1.908, 1.940, 1.989, 2.010, 2.012, 2.024, 2.047, 2.081, 2.085, 2.097, 2.136, 2.178, 2.181, 2.193, 2.200, 2.220, 2.301, 2.354, 2.359, 2.382, 2.409, 2.418, 2.430, 2.477, 2.500, 2.534, 2.572, 2.588, 2.591, 2.599, 2.660, 2.700, 2.700, 2.744, 2.845, 2.911, 2.952, 3.006, 3.021, 3.048, 3.059, 3.092, 3.152, 3.276, 3.289, 3.440, 3.447, 3.498, 3.705, 3.870, 3.896, 3.969, 4.000, 4.009, 4.196, 4.202, 4.311, 4.467, 4.490, 4.601, 4.697, 5.100, 5.120, 5.136, 5.141, 5.165, 5.260, 5.329, 5.778, 5.794, 6.285, 6.460, 6.917, 7.295, 7.701, 8.032, 8.142, 8.864, 9.263, 9.359, 10.801, 11.037, 11.504, 11.933, 11.998, 12.000, 14.153, 15.000, 15.398, 19.793, 23.150, 27.769, 28.288, 34.325, 42.691, 62.037, 77.839])
n = len(data)
# 提前计算所有x的对数和与阈值对数，避免重复计算
sum_log_x = np.log(data).sum()
log_xm = np.log(5e5)

# 定义负对数似然函数
def neg_ll(alpha):
    # 过滤alpha<=0的非法值
    if alpha <= 0:
        return np.inf
    return -(n * np.log(alpha) + n * alpha * log_xm - (alpha + 1) * sum_log_x)

# 优化时添加alpha>0的边界约束
ret = minimize(neg_ll, x0=np.array([1]), bounds=[(1e-6, None)])
print("估计的alpha值：", ret.x[0])

• 提前计算sum_log_x和log_xm这类常数项，避免优化迭代时重复计算，提升运行速度
• 采用numpy向量化运算替代遍历每个x的循环，大幅降低性能损耗
• 添加了alpha>0的边界约束，避免迭代到非法值触发计算错误

内容的提问来源于stack exchange，提问作者J. Grünwald