python - 有没有比使用 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函数的输入，估计结果会​​更好。有没有人已经在这部分工作或得到了更好的解决方案？