在python中拟合负二项式

Posted 2023-03-12

【中文标题】在python中拟合负二项式

【英文标题】：Fitting negative binomial in python

【发布时间】：2014-07-12 01:31:10

【问题描述】：

在 scipy 中，不支持使用数据拟合负二项分布 （可能是因为 scipy 中的负二项式只是离散的）。

对于正态分布，我会这样做：

from scipy.stats import norm
param = norm.fit(samp)

在任何其他库中是否有类似的“准备使用”功能？

【问题讨论】：

你找到这个问题的答案了吗？

请看@Eran 的回答。 statsmodels.discrete.discrete_model.NegativeBinomial

【参考方案1】：

Statsmodels 有discrete.discrete_model.NegativeBinomial.fit()，见这里： https://www.statsmodels.org/dev/generated/statsmodels.discrete.discrete_model.NegativeBinomial.fit.html#statsmodels.discrete.discrete_model.NegativeBinomial.fit

【参考方案2】：

我知道这个帖子已经很老了，但当前的读者可能想看看为此目的而制作的这个 repo：https://github.com/gokceneraslan/fit_nbinom

这里还有一个实现，虽然是更大包的一部分：https://github.com/ernstlab/ChromTime/blob/master/optimize.py

【参考方案3】：

不仅因为它是离散的，还因为最大似然拟合负二项式可能非常复杂，尤其是在附加位置参数的情况下。这就是为什么没有为它提供.fit() 方法的原因（以及Scipy 中的其他离散分布），这里是一个例子：

In [163]:

import scipy.stats as ss
import scipy.optimize as so
In [164]:
#define a likelihood function
def likelihood_f(P, x, neg=1):
    n=np.round(P[0]) #by definition, it should be an integer 
    p=P[1]
    loc=np.round(P[2])
    return neg*(np.log(ss.nbinom.pmf(x, n, p, loc))).sum()
In [165]:
#generate a random variable
X=ss.nbinom.rvs(n=100, p=0.4, loc=0, size=1000)
In [166]:
#The likelihood
likelihood_f([100,0.4,0], X)
Out[166]:
-4400.3696690513316
In [167]:
#A simple fit, the fit is not good and the parameter estimate is way off
result=so.fmin(likelihood_f, [50, 1, 1], args=(X,-1), full_output=True, disp=False)
P1=result[0]
(result[1], result[0])
Out[167]:
(4418.599495886474, array([ 59.61196161,   0.28650831,   1.15141838]))
In [168]:
#Try a different set of start paramters, the fit is still not good and the parameter estimate is still way off
result=so.fmin(likelihood_f, [50, 0.5, 0], args=(X,-1), full_output=True, disp=False)
P1=result[0]
(result[1], result[0])
Out[168]:
(4417.1495981801972,
 array([  6.24809397e+01,   2.91877405e-01,   6.63343536e-04]))
In [169]:
#In this case we need a loop to get it right
result=[]
for i in range(40, 120): #in fact (80, 120) should probably be enough
    _=so.fmin(likelihood_f, [i, 0.5, 0], args=(X,-1), full_output=True, disp=False)
    result.append((_[1], _[0]))
In [170]:
#get the MLE
P2=sorted(result, key=lambda x: x[0])[0][1]
sorted(result, key=lambda x: x[0])[0]
Out[170]:
(4399.780263084549,
 array([  9.37289361e+01,   3.84587087e-01,   3.36856705e-04]))
In [171]:
#Which one is visually better?
plt.hist(X, bins=20, normed=True)
plt.plot(range(260), ss.nbinom.pmf(range(260), np.round(P1[0]), P1[1], np.round(P1[2])), 'g-')
plt.plot(range(260), ss.nbinom.pmf(range(260), np.round(P2[0]), P2[1], np.round(P2[2])), 'r-')
Out[171]:
[<matplotlib.lines.Line2D at 0x109776c10>]

【讨论】：

我意识到这是一个非常临时的解决方案。优化功能并不总是有效