アヒル本8.4章ロジスティック回帰階層モデルをpymcでやってみる。

この例題は、学生の出席/欠席を、アルバイト好きか、学問への興味、天気などのデータから予測するもので、 学生差や科目差を考慮した階層モデルになっている。

StanとRでベイズ統計モデリング （Wonderful R） [ 松浦 健太郎 ]

これをpymcで実装してみたのだが、比較的複雑なモデルであったためか、収束が遅く、 MCMCサンプリングのパラメータ調整に骨が折れた。

# データの読み込み
df1 = pd.read_csv('data/data-attendance-4-1.txt')
df2 = pd.read_csv('data/data-attendance-4-2.txt')
# 出欠データに学生データをマージする
df = pd.merge(df2, df1, on='PersonID', how='left')

# データの整形
df1.Score = df1.Score/200
df.Weather2 = df.Weather.astype('category').cat.codes
# データ数
I = len(df) #出欠データ数 
N = len(df.PersonID.unique()) # 学生数`
C = len(df.CourseID.unique()) # 科目数

モデルの実装は以下の通り。 パラメータの事前分布も一様分布にするとMCMCの実行でエラーが発生したため、 範囲を狭めた。

with pm.Model() as model:
    # 学生
    b2 = pm.Normal('b2', mu=0, sd=5)
    b3 = pm.Normal('b3', mu=0, sd=5)
    s_p = pm.HalfCauchy('s_p', beta=5)
    x_p = pm.Normal('x_p', mu=b2*df1.A+b3*df1.Score, sd=s_p, shape=N)
    # 科目
    s_c = pm.HalfCauchy('s_c', beta=5)
    x_c = pm.Normal('x_c', mu=0, sd=s_c, shape=C)
    # 出欠
    b4 = pm.Normal('b4', mu=0, sd=5)
    x_w = pm.Deterministic('x_w', b4*df.Weather2)

    # 
    b1 = pm.Normal('b1', mu=0, sd=5)
    x = b1 + x_p[df.PersonID-1] + x_c[df.CourseID-1] + x_w
    q = pm.invlogit(x)
    
    y = pm.Bernoulli('y', q, observed=df.Y)

次はMCMCの実行なのだが、普通に実行すると、 「There were 70 divergences after tuning. Increase target_accept or reparameterize.」 と表示されたため、サンプリングのパラメータを調整した。 これでもまだ、有効サンプル数が少ないと言われているが、ひとまずこれで良しとした。

with model:
    trace = pm.sample(4000, tune=1000, nuts_kwargs=dict(target_accept=.95))

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [b1, b4, x_c, s_c, x_p, s_p, b3, b2]
Sampling 2 chains: 100%|██████████| 10000/10000 [07:30<00:00, 16.01draws/s]
The estimated number of effective samples is smaller than 200 for some parameters.

pm.summary(trace, varnames=['b1','b2','b3','b4','s_p','s_c','x_c'])

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
b1	0.891183	0.597252	0.034681	-0.294326	2.029407	249.037909	0.999991
b2	-0.799847	0.153195	0.003269	-1.097714	-0.500283	2634.308828	1.000426
b3	1.635974	0.594128	0.040074	0.494068	2.894429	159.238272	1.000977
b4	-0.392469	0.069501	0.000932	-0.522982	-0.251189	5431.473522	0.999893
s_p	0.354840	0.086419	0.002680	0.184643	0.525146	949.691355	1.003182
s_c	1.407389	0.417113	0.009598	0.742878	2.230000	1596.695327	1.000975
x_c__0	0.496923	0.528346	0.024537	-0.598436	1.513327	403.326425	1.000081
x_c__1	-1.294490	0.494636	0.025063	-2.297486	-0.304960	325.406740	1.000153
x_c__2	-0.985298	0.501408	0.024779	-2.014632	-0.015746	346.060896	1.000204
x_c__3	1.344727	0.524338	0.024928	0.220888	2.315912	376.628593	1.000045
x_c__4	-0.258379	0.500144	0.024979	-1.259820	0.730618	341.972609	1.000112
x_c__5	-1.931672	0.508829	0.025169	-2.940520	-0.919634	349.604963	0.999921
x_c__6	0.901908	0.512885	0.024878	-0.090079	1.943308	352.460165	1.000018
x_c__7	-0.536358	0.501140	0.024853	-1.523805	0.475292	348.076164	1.000092
x_c__8	1.207217	0.527864	0.025324	0.131012	2.239323	363.558944	1.000101
x_c__9	1.287300	0.530528	0.025280	0.200316	2.323838	371.970131	1.000043

pm.traceplot(trace, varnames=['b1','b2','b3','b4','s_p','s_c','x_c'])

パラメータの可視化。書籍ではviolinplotだが、ここではboxplotにした。

skip=10
estimate_params = pd.DataFrame(
    data = {'b1': trace['b1',::skip],
        'b2': trace['b2',::skip],
        'b3': trace['b3',::skip],
        'b4': trace['b4',::skip],
        's_p': trace['s_p',::skip],
        's_c': trace['s_c',::skip]
    })
sns.boxplot(x='variable', y='value', data=pd.melt(estimate_params))

大体同じような結果が得られたと思う。

以上