指数拟合：优化.curve_fit和stats.expon.fit产生不同的结果。

Question

我正试图根据我在这里读到的答案，使用两种不同的方法对指数分布进行直方图拟合。我对获得分布的比例参数的倒数感兴趣。

按照这里给出的答案(用python进行直方图拟合)，我使用的是 fit 的方法 scipy.stats.expon 的分布。

import glob
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sns

fig, ax = plt.subplots(5, 1, sharex = True)
j = 0

for files in glob.glob("data_*"):

    time = []
    hist = []

    with open(files, 'r') as f:
         for line in f:
             line = line.split(' ')
             time.append(float(line[0]))
             H.append(float(line[1]))

    P  = ss.expon.fit(H, floc = 0)
    T  = np.linspace(0,200, 1000)
    rP = ss.expon.pdf(T, *P)

    ax[j].plot(T, rP, lw = 3.0)
    ax[j].hist(H,bins = 30, alpha = 0.6, label = r"$\lambda = $" + str(1/P[1]), density = True, stacked = True)
    ax[j].set_yticks([])
    ax[j].legend()

    j = j +1 

sns.despine(top = True, left = True, right = True)
plt.xlabel("Time")
plt.show()

通过这样做，我得到了下面的图。

拟合效果不错，但我想知道不确定性误差的lambda值。没有任何关于如何在 stats.expon 文件。

这个问题在这里已经问过了 (有没有办法从scipy.stats.norm.fit中获取拟合参数的误差？). 公认的答案建议用curve_fit来代替直方图拟合。因此，按照这里的教程(https:/riptutorial.comscipyexample31081fitting-a-function-to-data-from-histogram。)，我尝试使用curve_fit。以下是修改后的代码（我插入了这些行，而不是使用 scipy.stats.expon）。

    def func(x, a):
        return a*np.exp(-a*x)

    bins = np.linspace(0, 200, 201)
    data_entries, bins = np.histogram(np.array(H), bins = bins)
    binscenters = np.array([0.5 * (bins[i] + bins[i + 1]) for i in range (len(bins)-1)])
    popt, pcov = curve_fit(func, xdata = binscenters, ydata = data_entries)

    ax[j].plot(T, func(T, *popt))
    ax[j].hist(H, bins = 30, alpha = 0.6, label = r"$\lambda = $" + str(popt[0]), density = True, stacked = True)

这个拟合产生的结果和 stats.expon.fit而这似乎（至少在质量上）更适合数据。

我是否用错了 curve_fit？我相信在某种程度上是这样的。curve_fit 和 expon.fit 应该会产生同样的结果。有什么方法可以让我从expon.fit中得到估计lambda的误差吗？我想计算数据的平均值和初始拟合估计的lambda之间的相对误差，但我不知道这是否正确。如果有任何提示，我将非常感激。

Answer 1

我设法解决了我的问题。事实证明，我缺少 density = True 关于 numpy.histogram.

该功能

def func(x, a):
        return a*np.exp(-a*x)

是指数的PDF。由于我的数据没有归一化（因此，不是PDF），所以用 curve_fit 是不正确的。通过这一修改，两个 ss.expon.fit 和 curve_fit 产生相同的lambda值。