python - 如何在 statsmodels 中使用 gamma GLM 的尺度和形状参数
问题描述
任务
我的数据如下所示:
我想使用statsmodels
. 使用这个模型,对于我的每个观察,我想计算观察到小于(或等于)该值的值的概率。换句话说,我想计算:
P(y <= y_i | x_i)
我的问题
如何从拟合的 glm 中获取形状和比例参数
statsmodels
?根据这个问题, statsmodels 中的比例参数没有以正常方式参数化。我可以将它直接用作 gamma 分布的输入scipy
吗?还是我需要先转型?如何使用这些参数(形状和比例)来获得概率?目前我正在使用
scipy
为每个生成分布x_i
并从中获取概率。请参阅下面的实现。
我目前的实现
import scipy.stats as stat
import patsy
import statsmodels.api as sm
# Generate data in correct form
y, X = patsy.dmatrices('y ~ x', data=myData, return_type='dataframe')
# Fit model with gamma family and log link
mod = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.log())).fit()
# Predict mean
myData['mu'] = mod.predict(exog=X)
# Predict probabilities (note that for a gamma distribution mean = shape * scale)
probabilities = np.array(
[stat.gamma(m_i/mod.scale, scale=mod.scale).cdf(y_i) for m_i, y_i in zip(myData['mu'], myData['y'])]
)
但是,当我执行此过程时,我得到以下结果:
目前预测的概率似乎都很高。图中的红线是预测平均值。但即使对于低于这条线的点,预测的累积概率也在 80% 左右。这让我想知道我使用的比例参数是否确实是正确的。
解决方案
在 R 中,您可以使用 1/dispersion 作为形状的估计值(查看这篇文章)。不幸的是,statsmodels 中的色散估计值的命名是scale
. 所以你确实取了这个的倒数来得到形状估计。我用下面的例子来展示它:
values = gamma.rvs(2,scale=5,size=500)
fit = sm.GLM(values, np.repeat(1,500), family=sm.families.Gamma(sm.families.links.log())).fit()
这是一个仅截距模型,我们检查截距和离散度(命名比例):
[fit.params,fit.scale]
[array([2.27875973]), 0.563667465203953]
所以平均值是exp(2.2599) = 9.582131
,如果我们使用 shape 作为 1/dispersion ,shape = 1/0.563667465203953 = 1.774096
这就是我们模拟的。
如果我使用模拟数据集,它工作得很好。这是它的样子,形状为 10:
from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
import patsy
import statsmodels.api as sm
import pandas as pd
_shape = 10
myData = pd.DataFrame({'x':np.random.uniform(0,10,size=500)})
myData['y'] = gamma.rvs(_shape,scale=np.exp(-myData['x']/3 + 0.5)/_shape,size=500)
myData.plot("x","y",kind="scatter")
然后我们像你一样拟合模型:
y, X = patsy.dmatrices('y ~ x', data=myData, return_type='dataframe')
mod = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.log())).fit()
mu = mod.predict(exog=X)
shape_from_model = 1/mod.scale
probabilities = [gamma(shape_from_model, scale=m_i/shape_from_model).cdf(y_i) for m_i, y_i in zip(mu,myData['y'])]
和情节:
fig, ax = plt.subplots()
im = ax.scatter(myData["x"],myData["y"],c=probabilities)
im = ax.scatter(myData['x'],mu,c="r",s=1)
fig.colorbar(im, ax=ax)
推荐阅读
- amazon-web-services - 我的免费套餐结束了吗?
- email - RFC 2822 格式 - 是否包含附件
- ios - tableView(_:heightForHeaderInSection:) 不工作
- haskell - `Functor`不是`Category`的超类是否有原因?
- ionic-framework - 我想在离子选择中选择离子选项后添加离子项目
- xml - 水壶(PDI)8.1读取xml错误
- scala - 在循环 DataFrame 并访问外部范围变量时获取 NullPointerException
- php - 没有这样的主机是已知的
- ibm-watson - Watson 助手中的逻辑确定对话
- c# - 带有 IsBackground 的 Lambda 线程