ガンマ回帰モデル（緑本6章）

久保先生の緑本6章には、応答変数が非負であったり、ばらつきを正規分布で仮定できない場合として、ガンマ分布を用いた例がある。ガンマ分布のGLMもRやpythonのstatsmodelsを使った例はよく見るが、pymcで実装した例は見当たらなかっため、実装してみることにした。

ガンマ分布の確率密度関数は、以下と定義される。

$\displaystyle p(y|s,r)=\frac{r^s}{\Gamma(s)} y^{s-1} \exp(-ry)$

$s$ はshapeパラメータ、
$r$ はrateパラメータと呼ばれる。
また、 $1/r$ は、scaleパラメータと呼ばれる。
平均は、 $s/r$ となる。

# 事前にRDataをtxtに変換したファイルを読み込む。
df = pd.read_table('kubo_6_d.txt', sep=' ')
df.head()

	x	y
1	0.001000	0.000887
2	0.017306	0.023465
3	0.033612	0.069876
4	0.049918	0.034340
5	0.066224	0.026520

plt.scatter(df.x,df.y)
plt.xlabel("葉重量$x_i$")
plt.ylabel("花重量$y_i$")

f:id:gaiasky:20180904200742p:plain

xが大きくなるにつれて、yのばらつきが大きくなっており、正規分布のGLMは不適切。

ガンマ回帰モデル

ガンマ回帰モデルでは、データはガンマ分布に従い、またリンク関数には対数リンク関数が用いられる。線形予測子は以下の通り。

$\displaystyle \begin{align} \mu_i &= \exp(a+b \log x_i) \\ \log \mu_i &= a+b \log x_i \end{align}$

まずは、statsmodelsのglmを使って、ガンマ回帰。

# gamma regression model
model = smf.glm('y~np.log(x)',data=df,family=sm.families.Gamma(link=sm.families.links.log))
result = model.fit()
result.summary()

Generalized Linear Model Regression Results
Dep. Variable:	y	No. Observations:	50
Model:	GLM	Df Residuals:	48
Model Family:	Gamma	Df Model:	1
Link Function:	log	Scale:	0.32508
Method:	IRLS	Log-Likelihood:	58.471
Date:	Tue, 04 Sep 2018	Deviance:	17.251
Time:	19:19:34	Pearson chi2:	15.6
No. Iterations:	18	Covariance Type:	nonrobust

	coef	std err	z	P>\|z\|	[0.025	0.975]
Intercept	-1.0403	0.119	-8.759	0.000	-1.273	-0.808
np.log(x)	0.6832	0.068	9.992	0.000	0.549	0.817

グラフ化する。

intercept = result.params["Intercept"]
nplog = result.params["np.log(x)"]
pred = np.exp(intercept + np.log(df.x) * nplog)
plt.scatter(df.x, df.y)
plt.plot(df.x,pred,c="blue")
plt.xlabel("葉重量$x_i$")
plt.ylabel("花重量$y_i$")
plt.title('statsmodels.GLMによるガンマ回帰')

f:id:gaiasky:20180904201029p:plain

pymcによるガンマ回帰。

$\mu = s/r$ であるから、 $r=s/\mu$ となる。
$\mu_i = \exp(a+b \log x_i)$ であるから、各 $x_i$ のrateパラメータは、
$r_i=s/\exp(a+b \log x_i)$ と計算される。

with pm.Model() as model:
    a = pm.Normal('a', mu=0, sd=100)
    b = pm.Normal('b', mu=0, sd=100)
    
    shape = pm.Uniform('shape',0.001, 100)
    rate = pm.Deterministic('rate', shape/pm.math.exp(a+b*np.log(df.x)))
    
    # alpha (float) – Shape parameter (alpha > 0).
    # beta (float) – Rate parameter (beta > 0).
    y = pm.Gamma('y', alpha=shape, beta=rate, observed=df.y)

pm.model_to_graphviz(model)

f:id:gaiasky:20180904201358p:plain

with model:
    trace = pm.sample(2000,start=pm.find_MAP(method='Powell'), step=pm.NUTS())

logp = 42.641, ||grad|| = 3.2989: 100%|██████████| 250/250 [00:00<00:00, 876.42it/s] 
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [shape, b, a]
Sampling 2 chains: 100%|██████████| 5000/5000 [00:07<00:00, 677.39draws/s]

pm.traceplot(trace)

f:id:gaiasky:20180904201442p:plain

グラフ化する。

# 各xの平均値について、ベイズ信用区間95%の範囲を示す。50%を直線で示す。
pred = [pm.quantiles(np.exp(trace['a']+trace['b']*np.log(x)), 
                     qlist=[2.5, 50, 97.5]).values() for x in df.x]
y1,med,y2 = map(list,zip(*pred))
plt.scatter(df.x, df.y)
plt.plot(df.x, med, c="blue")
plt.fill_between(x=df.x, 
                 y1=y1,
                 y2=y2,
                 alpha=0.2, color='blue')
plt.xlabel('葉重量$x_i$')
plt.ylabel('花重量$y_i$')
plt.title('pymcによるガンマ回帰(平均値のベイズ信用区間95%範囲)')

f:id:gaiasky:20180904201531p:plain

各xの予測区間を示す。

from scipy.stats import gamma
# 各xにおける90%予測区間を計算する。パラメータは事後分布の平均値を使用。
pred = []
for x in df.x:
    shape = trace['shape'].mean()
    rate = shape / np.exp(trace['a'].mean()+trace['b'].mean()*np.log(x))
    scale =  1/rate # scipyのgamma分布のパラメータはscaleパラメータ
    lower, upper = gamma.ppf(q=[0.05, 0.95], a=shape, scale=scale)
    pred.append((lower,upper))
# グラフ化
y1,y2 = map(list, zip(*pred))
plt.scatter(df.x, df.y)
plt.plot(df.x, med, c="blue")
plt.fill_between(x=df.x, 
                 y1=y1,
                 y2=y2,
                 alpha=0.2, color='blue')
plt.xlabel('葉重量$x_i$')
plt.ylabel('花重量$y_i$')
plt.title('pymcによるガンマ回帰(90%予測区間)')