介绍

  数学建模中很多时候，我们有某个随机变量 $X$ 的若干样本 $X_1,X_2,\cdots ,X_n$，想要还原随机变量 $X$ 的概率密度函数 $f(x)$。诚然，高斯核密度估计可以做到这一点，这里我们介绍一个新的 Python 库 Fitter。
  Fitter 的原理是，用指定的某些分布去匹配连续性随机变量 $X$ 的观测值，对于每一种分布都拟合出响应的参数、计算出均方误差。一般选择均方误差最小的分布作为拟合结果。核心方法Fitter的函数原型为：

Fitter(data, xmin=None, xmax=None, bins=100, 
        distributions=None, verbose=True, timeout=10)

  其中：

• data是已知的 $X$ 的观测值集合，或者说是样本点，即上面的 $X_1,\cdots ,X_n$。
• xmin表示忽略观测值data中小于xmin的数据。
• xmax表示忽略观测值data中大于xmax的数据。
• bins是累计直方图的组数（直方图里有几个矩形），值越大拟合曲线越平滑。
• verbose表示是否显示 Fitter 拟合过程中的一些提示信息，设置verbose=False将会使你的控制台更干净。
• timeout表示对于一种分布尝试的最长时间。因为各种各样的原因，给出的数据data可能不适合用某种分布去拟合，当 Fitter 在这种分布上拟合的时间超过timeout秒时，将会自动放弃拟合这种分布。
• distributions是一个列表，用来告诉 Fitter你想尝试哪些概率密度函数（比如正态分布，$t$ 分布，$\beta$ 分布等）。默认情况下会尝试所有的 scipy 分布（$106$ 种）。常用分布distributions=['norm', 't', 'laplace', 'cauchy', 'chi2', 'expon', 'exponpow', 'gamma', ' lognorm', 'uniform']。（比如norm是正态分布，t是 $t$ 分布，laplace是拉普拉斯分布，cauchy是柯西分布……）

参考文献：Python fitter包：拟合数据样本的分布_fitter的80多个分布有哪些-CSDN博客

代码实例

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
from fitter import Fitter

# 直接生成 2000 个 beta 分布的样本点
data = np.array(np.random.beta(a=2, b=5, size=2000))

# Euclid 字体是 latex 公式使用的字体
plt.rcParams['font.family'] = 'Euclid'

# 尝试使用 distributions 里的分布去拟合样本点
f = Fitter(data, distributions=['norm', 't', 'laplace','beta','lognorm'])

f.fit()

print(f.summary())  # 返回排序好的分布拟合质量（拟合效果从好到坏）,并绘制数据分布和Nbest分布
"""
          sumsquare_error         aic  ...  ks_statistic     ks_pvalue
beta           19.103771  164.735384  ...      0.016081  6.791340e-01
lognorm        22.106143  143.794008  ...      0.026787  1.113541e-01
t              30.674133  248.908752  ...      0.060242  9.434693e-07
norm           30.674533  246.915057  ...      0.060244  9.426322e-07
laplace        52.044523  184.287028  ...      0.068199  1.562740e-08
可以看到用 beta 分布去拟合，均方误差最小。
"""

print(f.df_errors)  # 返回这些分布的拟合质量（均方根误差的和）
"""
sumsquare_error         aic  ...  ks_statistic     ks_pvalue
norm           30.674533  246.915057  ...      0.060244  9.426322e-07
laplace        52.044523  184.287028  ...      0.068199  1.562740e-08
beta           19.103771  164.735384  ...      0.016081  6.791340e-01
lognorm        22.106143  143.794008  ...      0.026787  1.113541e-01
t              30.674133  248.908752  ...      0.060242  9.434693e-07
"""

print(f.fitted_param)  # 返回拟合分布的参数
"""
{'norm': (0.2796987317013659, 0.15616306964918625), 
'laplace': (0.26198195036801175, 0.12567921087793762), 
'beta': (1.9506594116888034, 4.9951825006926605, 0.0024740416891275907, 0.9862799904927366), 
'lognorm': (0.31353430150353023, -0.21970151525943354, 0.47572327635046574), 
't': (850256.2207528697, 0.2796983453401272, 0.15616418590110978)}
"""

print(f.fitted_pdf)  # 使用最适合数据分布的分布参数生成的概率密度
"""
结果是返回从分布到 array 的字典。array 是一系列概率密度函数的取值，用来作概率密度曲线。
"""

print(f.get_best(method='sumsquare_error'))  # 返回最佳拟合分布及其参数
"""
{'beta': 
    {'a': 1.9506594116888034, 
    'b': 4.9951825006926605, 
    'loc': 0.0024740416891275907, 
    'scale': 0.9862799904927366}
}
"""

f.hist()  # 绘制组数=bins的标准化直方图
# 下面这个方法绘制分布的概率密度函数，按照方差从小到大取前 Nbest 个。
# 因为前面的 f.summary() 也会画图，所以就不启用这个方法。
# 顺带一提，前面的  f.summary() 也可以设置 Nbest 参数。
# f.plot_pdf(names=None, Nbest=3, lw=2)
plt.show()

  画出来的图如下所示，$\beta$ 分布是均方误差最小的分布，可以认为是最合适的分布。（因为数据就是按照 $\beta$ 分布生成的，所以这个结果很合理。）同时，对数正态分布也比较合适，$t$ 分布次之。

注意：

• Fitter方法的distributions参数决定了需要拟合哪些分布，默认拟合里面所有的 $106$ 种分布。默认情况下会跑比较长的时间，而且往往给出最佳的分布是一个你从来没见过的分布。不建议使用默认情况。
• f.summary()默认参数Nbest=5，给出 $5$ 个最好拟合的结果。如果想要所有分布的结果，请设置Nbest=len(f.distributions)，即将f.summary()改为f.summary(Nbest=len(f.distributions))。