statsmodels - 绘制拟合分布答案

【问题标题】：statsmodels - plotting the fitted distributionstatsmodels - 绘制拟合分布
【发布时间】：2014-05-23 09:48:59
【问题描述】：

以下代码使用statsmodels 拟合过度简化的广义线性模型

model = smf.glm('Y ~ 1', family=sm.families.NegativeBinomial(), data=df)
results = model.fit()

这给出了系数和标准错误：

               coef stderr   
Intercept    2.9471  0.120

现在我想以图形方式比较变量 Y（直方图）的真实分布与来自模型的分布。

但我需要两个参数 r 和 p 来评估 stats.nbinom(r,p) 并绘制它。

有没有办法从拟合结果中检索参数？如何绘制 PMF？

【问题讨论】：

glm 不估计 NegativeBinomial 的形状参数。您可以改用discrete_model NegativeBinomial。
能否提供一个代码示例？
statsmodels 和 scipy.stats 的参数化是不同的，我仍然需要弄清楚如何匹配它们。不幸的是，其中一些“额外”在 statsmodels 中仍然缺失。

标签： python statistics scipy modeling statsmodels

【解决方案1】：

statsmodels 中的广义线性模型 GLM 目前不估计负二项分布的额外参数。负二项式属于指数分布族，仅适用于固定形状参数。

但是，statsmodels 在discrete_model 中也将负二项式作为最大似然模型，用于估计所有参数。

计数回归的负二项式的参数化是根据平均值或期望值，这与 scipy.stats.nbinom 中的参数化不同。实际上，负二项式计数回归有两种不同的常用参数化，通常称为nb1和nb2

这是一个快速编写的脚本，可以从估计的参数中恢复 scipy.stats.nbinom 参数 n=size 和 p=prob。获得 scipy.stats.distribution 的参数后，您可以使用所有可用的方法、rvs、pmf 等。

类似的东西应该在 statsmodels 中可用。

在一些示例运行中，我得到了这样的结果

data generating parameters 50 0.25
estimated params           51.7167511571 0.256814610633
estimated params           50.0985814878 0.249989725917

此外，由于底层的指数重新参数化，scipy 优化器有时会出现收敛问题。在这些情况下，提供更好的起始值或使用 Nelder-Mead 作为优化方法通常会有所帮助。

import numpy as np
from scipy import stats
import statsmodels.api as sm

# generate some data to check
nobs = 1000
n, p = 50, 0.25
dist0 = stats.nbinom(n, p)
y = dist0.rvs(size=nobs)
x = np.ones(nobs)

loglike_method = 'nb1'  # or use 'nb2'
res = sm.NegativeBinomial(y, x, loglike_method=loglike_method).fit(start_params=[0.1, 0.1])

print dist0.mean()
print res.params

mu = res.predict()   # use this for mean if not constant
mu = np.exp(res.params[0])   # shortcut, we just regress on a constant
alpha = res.params[1]

if loglike_method == 'nb1':
    Q = 1
elif loglike_method == 'nb2':    
    Q = 0

size = 1. / alpha * mu**Q
prob = size / (size + mu)

print 'data generating parameters', n, p
print 'estimated params          ', size, prob

#estimated distribution
dist_est = stats.nbinom(size, prob)

顺便说一句：我之前遇到过这个但没时间看 https://github.com/statsmodels/statsmodels/issues/106

【讨论】：