如何在统计模型中使用Gamma GLM的比例和形状参数

如何解决如何在统计模型中使用Gamma GLM的比例和形状参数

任务

我有如下数据：

我想使用statsmodels从gamma族中拟合出广义线性模型（glm）。使用此模型，对于我的每个观察值，我都希望计算观察到一个小于（或等于）该值的概率。换句话说，我想计算：

P（y

我的问题

• 如何从statsmodels中拟合的glm中获取形状和比例参数？根据{{​​3}}，未按常规方式对statsmodels中的scale参数进行参数化。我可以直接将其用作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'])]
)

但是，当我执行此过程时，会得到以下结果：

this question

目前，预测的概率似乎都很高。图中的红线是预测的平均值。但是，即使对于该线以下的点，预测的累积概率也约为80％。这使我想知道我使用的scale参数是否确实正确。

解决方法

在R中，您可以使用1 /色散作为形状的估计值（请检查此post）。不幸的是，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，如果我们使用形状作为1 /色散，则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)