任务
我有如下数据:
我想使用statsmodels
将一个广义线性模型(glm)从gamma族中拟合出来。使用这个模型,对于我的每一次观察,我想计算观察到小于(或等于)该值的概率。换句话说,我想计算:
P(y <= y_i | x_i)
我的问题
如何从statsmodels
中拟合的glm中获取形状和比例参数?根据this question,statsmodels中的比例参数不是以正常方式参数化的。我可以直接使用它作为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/分散(检查此post)作为形状的估计来获得。不幸的是,statsmodels中分散估计的命名是
scale
。所以你需要取这个的倒数来得到形状估计。我用下面的一个例子来说明:这是一个仅限截距的模型,我们检查截距和色散(命名比例):
所以平均值是
exp(2.2599) = 9.582131
,如果我们使用形状作为1/色散,shape = 1/0.563667465203953 = 1.774096
,这就是我们模拟的如果我使用模拟数据集,它工作得非常好。这是它的外观,形状为10:
然后,我们像您一样拟合模型:
并绘制:
相关问题 更多 >
编程相关推荐