有没有比使用 Scipy 更好地拟合 beta 素数分布的解决方案？

Question

我试图使用 python 将 beta 素数分布拟合到我的数据中。因为有

scipy.stats.betaprime.fit

，我尝试了这个：

import numpy as np
import math
import scipy.stats as sts
import matplotlib.pyplot as plt

N  = 5000
nb_bin = 100
a = 12; b = 106; scale = 36; loc = -a/(b-1)*scale
y = sts.betaprime.rvs(a,b,loc,scale,N)
a_hat,b_hat,loc_hat,scale_hat = sts.betaprime.fit(y)
print('Estimated parameters: \n a=%.2f, b=%.2f, loc=%.2f, scale=%.2f'%(a_hat,b_hat,loc_hat,scale_hat))

plt.figure()
count, bins, ignored = plt.hist(y, nb_bin, normed=True)
pdf_ini = sts.betaprime.pdf(bins,a,b,loc,scale)
pdf_est  = sts.betaprime.pdf(bins,a_hat,b_hat,loc_hat,scale_hat)
plt.plot(bins,pdf_ini,'g',linewidth=2.0,label='ini');plt.grid()
plt.plot(bins,pdf_est,'y',linewidth=2.0,label='est');plt.legend();plt.show()

它向我显示的结果是：

Estimated parameters:
 a=9935.34, b=10846.64, loc=-90.63, scale=98.93

这与原始图和 PDF 中的图有很大不同：

如果我给出

loc

的真实值和尺度作为

fit

函数的输入，估计结果会更好。有人已经研究过这部分或者有更好的解决方案吗？

Answer 1

原帖中的代码没有修复随机数生成器种子，因此问题不可重现。可复制的版本是：

import numpy as np
import math
import scipy.stats as sts
import matplotlib.pyplot as plt

rng = np.random.default_rng(83492456935469354704)

N  = 5000
nb_bin = 100
a = 12; b = 106; scale = 36; loc = -a/(b-1)*scale
y = sts.betaprime.rvs(a, b, loc, scale, N, random_state=rng)
a_hat, b_hat, loc_hat, scale_hat = sts.betaprime.fit(y)
print('Estimated parameters:\n'
      f'a={a_hat:.2f}, b={b_hat:.2f}, loc={loc_hat:.2f}, scale={scale_hat:.2f}')

plt.figure()
count, bins, ignored = plt.hist(y, nb_bin, density=True)
pdf_ini = sts.betaprime.pdf(bins, a, b, loc, scale)
pdf_est = sts.betaprime.pdf(bins, a_hat, b_hat, loc_hat, scale_hat)
plt.plot(bins, pdf_ini, 'g', linewidth=2.0, label='ini'); plt.grid()
plt.plot(bins, pdf_est, 'y--', linewidth=2.0, label='est'); plt.legend(); plt.show()

拟合后的 PDF 的绘图与原始 PDF 非常匹配。参数不同，但我们可以期望拟合函数做的就是优化参数以最大化或最小化数据的某些函数；这可能会也可能不会导致对基础参数的良好估计。我们可以看到

betaprime.fit

在这里发挥作用，因为当我们沿任一方向扰动参数值时，负对数似然函数会增加。

fig, ax = plt.subplots(2, 2)
ax = ax.ravel()

eps = np.linspace(-5, 5, 300)
params = dict(a=a_hat, b=b_hat, loc=loc_hat, scale=scale_hat)

for hi, axi, param_name in zip(np.eye(4), ax, params):
  param_values = list(params.values())
  nllfs = [stats.betaprime.nnlf(param_values + epsi*hi, y) for epsi in eps]
  axi.plot(eps, nllfs)
  axi.set_xlabel(f'perturbation to `{param_name}`')

fig.suptitle('Effect of parameter perturbations on NLLF')
fig.tight_layout()

似乎多种参数值都可以给出大致相同的负对数似然函数值，即上面的

nllf0 = 8101.751831213405

。例如，使用

stats.fit

（使用随机优化器），如果我们运行几次，我们可以获得几个非常不同的局部最优值，所有 NLLF 值都大致相同。

FitParams(a=16.437726024127105, b=55.11441645100978, loc=-4.464265239903504, scale=14.750853658609275) 8101.030972971992
FitParams(a=67.47536806626265, b=35.31736834983663, loc=-5.97011638305157, scale=3.044776624686649) 8102.577960289996
FitParams(a=39.76533774830455, b=37.594009970143354, loc=-5.467151220796732, scale=5.046462382947629) 8102.258806849436
FitParams(a=18.873466185063798, b=48.85176088955696, loc=-4.6217525946033, scale=11.760288109582378) 8101.294928499379
FitParams(a=14.430145578727348, b=61.61794106850824, loc=-4.300115416751307, scale=18.134272628062437) 8100.849263206925

因此，我不确定是否还有更多工作要做，除非您可以根据对问题的了解来修复某些参数值或提供更严格的界限。例如，如果您可以在参数周围提供合理的界限，

scipy.stats.fit

可以找到比原始参数具有更好（更低）负对数似然函数的参数。

limits = np.asarray([-10, 10])
bounds = dict(a=a+limits, b=b+limits, loc=loc+limits, scale=scale+limits)
res = stats.fit(stats.betaprime, data=y, bounds=bounds)
ref = stats.betaprime.nnlf((a, b, loc, scale), y)
print(res.nllf(), res.params)
# 8100.411619791135 FitParams(a=11.830670229468097, b=98.43447339490054, loc=-4.094123769559834, scale=33.8539065644368)
print(ref, (a, b, loc, scale))
# 8100.927379232642 (12, 106, -4.114285714285714, 36)