求助：如何用Python的MLE拟合双指数分布？代码不收敛排查

嘿，我看了你的代码，立刻发现问题出在似然函数的定义完全错误——这就是为啥优化算法根本没干活，直接返回你的初始猜测值！咱们一步步来修正：

核心问题分析 #

你原来的MLE函数错误地用了正态分布的PDF来计算似然，但你的数据是直接从双指数混合分布生成的，完全和正态分布无关。更关键的是，你把原始观测数据ydata2当成了“服从均值为拟合曲线yPred的正态分布”，这逻辑完全不符合你的数据生成过程，导致优化算法计算的梯度几乎为零，直接一步就终止了，根本没去寻找最优参数。

而且你加的sd参数完全多余——双指数混合分布的参数只有三个：两个成分的权重，以及两个指数分布的速率（或scale）。

修正后的代码及解释 #

我把你的代码改好了，关键部分都加了注释：

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import scipy.stats as stats

size = 300

def simu_dt():
    ## simulate Exp2 data
    np.random.seed(0)
    x = np.random.rand(size)
    data = []
    for n in x:
        if n < 0.6: # 60%来自第一个指数分布，scale=20 → k1=1/20=0.05
            data.append(np.random.exponential(scale=20))
        else: # 40%来自第二个指数分布，scale=500 → k2=1/500=0.002
            data.append(np.random.exponential(scale=500))
    return np.array(data)

ydata2 = simu_dt() 
# 数据裁剪不变
ydata2 = ydata2[np.where(2 < ydata2)]
ydata2 = ydata2[np.where(ydata2 < 3000)]

# 生成直方图用于绘图（仅可视化，不参与MLE计算）
bins = 10 ** np.linspace(np.log10(np.min(ydata2)), np.log10(np.max(ydata2)), 10)
counts, bin_edges = np.histogram(ydata2, bins=bins)
bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2
bin_width = (bin_edges[1:] - bin_edges[:-1])
counts = counts / bin_width / np.sum(counts)

# 生成用于绘制拟合曲线的x值
x1 = np.logspace(np.log10(np.min(ydata2)), np.log10(np.max(ydata2)), 100)

def MLE(params):
    """ 正确的双指数混合分布负对数似然函数 """
    a1, k1, k2 = params
    # 计算每个数据点在双指数混合分布下的PDF
    pdf = (1 - a1) * k1 * np.exp(-k1 * ydata2) + a1 * k2 * np.exp(-k2 * ydata2)
    # 避免取对数时出现0，加一个极小值
    pdf = np.maximum(pdf, 1e-12)
    # 负对数似然：-sum(log(pdf))
    negLL = -np.sum(np.log(pdf))
    return negLL

# 初始猜测：a1是第二个成分的权重（真实值0.4），k1=1/20=0.05，k2=1/500=0.002
guess = np.array([0.4, 0.05, 0.002])
# 参数边界：a1在(0,1)，k1和k2必须大于0（设极小下限避免数值问题）
bnds = ((1e-6, 0.999), (1e-6, 1), (1e-6, 1))

# 用L-BFGS-B优化（适合带边界的优化，比SLSQP更稳定）
results = optimize.minimize(MLE, guess, method='L-BFGS-B', bounds=bnds)
print("拟合结果：")
print(results)

A1, K1, K2 = results.x
# 计算拟合曲线的PDF值
y_fitted = (1 - A1)*K1*np.exp(-K1*x1) + A1*K2*np.exp(-K2*x1)

# 绘图
plt.plot(bin_centres, counts, 'ko', label="实际数据")
plt.plot(x1, y_fitted, c='r', linestyle='dashed', label="拟合曲线")
plt.xlabel("Dwell Times (s)")
plt.ylabel("Probability")
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.show()

修正后的关键变化 #

1. 移除多余的sd参数：双指数混合分布不需要这个参数，你的模拟数据没有噪声，是纯分布采样。
2. 正确计算似然：直接用原始观测数据ydata2计算每个点的混合PDF，然后取对数求和再取负——这才是MLE的正确计算方式。
3. 优化方法调整：用L-BFGS-B替代SLSQP，对于带边界的连续优化问题，它的稳定性更好。
4. 避免数值问题：给PDF加了极小值1e-12，防止取对数时出现log(0)的错误。

预期结果 #

运行修正后的代码，你会看到优化过程会迭代多次，最终拟合出的参数会非常接近真实值：

• A1（第二个成分的权重）≈0.4
• K1≈0.05（对应scale=20）
• K2≈0.002（对应scale=500）
  拟合曲线也会完美贴合你的直方图数据。

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