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

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

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

• はshapeパラメータ、
• はrateパラメータと呼ばれる。
• また、は、scaleパラメータと呼ばれる。
• 平均は、となる。

# 事前に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$")

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

ガンマ回帰モデル

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

まずは、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によるガンマ回帰')

pymcによるガンマ回帰。

であるから、となる。
であるから、各のrateパラメータは、
と計算される。

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)

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)

グラフ化する。

# 各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%範囲)')

各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%予測区間)')

本とほとんど同じ結果が得られた。

以上

2018年Jリーグの各チームの戦力をモデル化する

前回の記事で、pymc3の公式HPに、ラグビーチームの戦力をモデル化する例題があることに触れた。
https://docs.pymc.io/notebooks/rugby_analytics.html

これは、各チームの攻撃力、防御力にホーム効果の有無を付けたものになっている。 これはラグビーだけでなく、サッカーなどにも使えるのでは？と思い、 Jリーグの各チームの戦力をモデルに当てはめて推定してみた。 公式HPと比べた本記事の新しさは以下の2点。

1. J1リーグの2018年の試合を用いた。
2. 未消化試合の結果を予測をした。

Jリーグの試合結果は、以下のサイトからダウンロードできる。 2018年のJ1全試合の結果の表をコピーし、tsvファイルとして保存した。
https://data.j-league.or.jp/SFMS01/

# データの読み込み
df = pd.read_table('./j_league/j_league_2018.tsv', delim_whitespace=False)
# columnsから空白を除去する。read_tableで、delim_whitespace=Trueにすると、未対戦結果のデータがおかしくなる。
new_col = []
for c in df.columns:
    c2 = c.strip()
    new_col.append(c2)
df.columns = new_col

df.head()

	年度	大会	節	試合日	K/O時刻	ホーム	スコア	アウェイ	スタジアム	入場者数	インターネット中継・TV放送
0	2018	Ｊ１	第１節第１日	02/23(金)	20:03	鳥栖	1-1	神戸	ベアスタ	19,633	ＤＡＺＮ
1	2018	Ｊ１	第１節第２日	02/24(土)	14:03	FC東京	1-1	浦和	味スタ	35,951	ＤＡＺＮ／ＮＨＫ　ＢＳ１／ＴＯＫＹＯ　ＭＸ（録）
2	2018	Ｊ１	第１節第２日	02/24(土)	14:03	広島	1-0	札幌	Ｅスタ	17,026	ＤＡＺＮ／ＮＨＫ広島／ＮＨＫ札幌
3	2018	Ｊ１	第１節第２日	02/24(土)	14:33	Ｇ大阪	2-3	名古屋	吹田Ｓ	28,681	ＤＡＺＮ
4	2018	Ｊ１	第１節第２日	02/24(土)	16:03	湘南	2-1	長崎	ＢＭＷス	12,148	ＤＡＺＮ

データの前処理を行う。

# ヘッダ名の変更
df = df.rename(columns={'ホーム':'home_team', 'アウェイ':'away_team'})

# チーム名をコード化する。
ch = pd.Categorical(df.home_team)
ca = pd.Categorical(df.away_team)
team_dict = dict(zip(ch.tolist(), ch.codes))
team_labes = [ x[0] for x in (sorted(team_dict.items(), key=lambda x:x[1]))]
df['i_home'] = ch.codes
df['i_away'] = ca.codes

df['スコア'].unique()

array(['1-1 ', '1-0 ', '2-3 ', '2-1 ', '0-0 ', '0-3 ', '3-3 ', '2-0 ',
       '0-1 ', '2-4 ', '2-2 ', '1-2 ', '1-3 ', '3-2 ', '3-0 ', '2-5 ',
       '3-1 ', '0-2 ', '4-4 ', '4-1 ', '4-2 ', '5-2 ', '4-0 ', '2-8 ',
       '6-2 ', '1-4 ', '中止 ', '5-1 ', 'vs '], dtype=object)

# 対戦済の結果を抽出する
df_train = df.drop(index=df[(df['スコア']=='vs ')|(df['スコア']=='中止 ')].index)

# スコアをホームチームのスコアとアウェイチームのスコアに分割する。
scores = df_train['スコア'].str.split('-', expand=True)
scores.columns = ['home_score','away_score']
df_train = pd.concat([df_train,scores], axis=1)
# 数値に変換
df_train.home_score = df_train.home_score.astype(int)
df_train.away_score = df_train.away_score.astype(int)

# 未対戦の試合を抽出する。
df_eval = df.drop(index=df[df['スコア']!='vs '].index)

対戦済結果の統計情報。やはりホームの方が有利な模様。

data = {
    'ホームが勝った確率(%)':(100*len(df_train[df_train.home_score>df_train.away_score])/len(df_train)),
    'ホームが負けた確率(%)':(100*len(df_train[df_train.home_score<df_train.away_score])/len(df_train)),
    'ホームが引分る確率(%)':(100*len(df_train[df_train.home_score==df_train.away_score])/len(df_train))
}
pd.Series(data).to_frame().T

	ホームが勝った確率(%)	ホームが負けた確率(%)	ホームが引分る確率(%)
0	42.105263	34.449761	23.444976

戦力のモデル化（公式HPと同じ）

with pm.Model() as model:
    # prior
    home = pm.Flat('home')
    intercept = pm.Flat('intercept')
    
    sd_att = pm.HalfStudentT('sd_att', nu=3, sd=2.5)
    sd_def = pm.HalfStudentT('sd_def', nu=3, sd=2.5)

    # team-specific model parameters
    atts_star = pm.Normal('atts_star', mu=0, sd=sd_att, shape=len(df.home_team.unique()))
    defs_star = pm.Normal('defs_star', mu=0, sd=sd_def, shape=len(df.home_team.unique()))

    atts = pm.Deterministic('atts', atts_star - tt.mean(atts_star))
    defs = pm.Deterministic('defs', defs_star - tt.mean(defs_star))
    
    home_theta = tt.exp(intercept+home+atts[df_train.i_home]+defs[df_train.i_away])
    away_theta = tt.exp(intercept+atts[df_train.i_away]+defs[df_train.i_home])
    
    #
    home_points = pm.Poisson('home_points', mu=home_theta, observed=df_train.home_score)
    away_points = pm.Poisson('away_points', mu=away_theta, observed=df_train.away_score)

実行時間優先のため、ADVIでパラメータ推定。

with model:
    inference = pm.ADVI()
    approx = pm.fit(n=20000, method=inference)
    trace = approx.sample(draws=5000)

Average Loss = 625.17: 100%|██████████| 20000/20000 [00:44<00:00, 452.96it/s]
Finished [100%]: Average Loss = 625.17

pm.traceplot(trace)

各チームの攻撃力と防御力。

まず、攻撃力（大きい方が点が入りやすい）。

plt.figure(figsize=(6,6))
pm.forestplot(trace, varnames=['atts'], ylabels=team_labes)

次に防御力（大きい方が点が入りやすい）。

plt.figure(figsize=(6,6))
pm.forestplot(trace, varnames=['defs'], ylabels=team_labes)

未消化試合の結果予測。

# 対戦結果のシミュレーション
def predict_vs(trace, teams,i_home, i_away):
    df_atts = pd.DataFrame(pm.stats.quantiles(trace['atts'])[50],
                           columns=['hpd_median'],
                           index=teams.i.values)
    df_defs = pd.DataFrame(pm.stats.quantiles(trace['defs'])[50],
                           columns=['hpd_median'],
                           index=teams.i.values)
    intercept = pm.stats.quantiles(trace['intercept'])[50]
    home = pm.stats.quantiles(trace['home'])[50]

    home_theta = tt.exp(intercept+home+df_atts.loc[i_home]+df_defs.loc[i_away])
    away_theta = tt.exp(intercept+df_atts.loc[i_away]+df_defs.loc[i_home])
    
    home_p = pm.Poisson.dist(mu=home_theta).random(size=1000)
    away_p = pm.Poisson.dist(mu=away_theta).random(size=1000)
    return home_p, away_p

# チーム名とindex
#以下コードは間違い。チーム名とindex番号の対応が間違っている。
#teams = df_train.home_team.unique()
#teams = pd.DataFrame(teams, columns=['team'])
#teams['i'] = teams.index
# こっちが正しい。（8/31）
teams = pd.Series(data=team_dict).to_frame(name='i').sort_values(by='i').reset_index().rename(columns={'index':'team'})

vs_pred = []
for ht,at in zip(df_eval.home_team, df_eval.away_team):
    #print('%s vs %s'%(ht,at))
    i_ht = teams[teams.team==ht].i.values[0]
    i_at = teams[teams.team==at].i.values[0]
    # シミュレーション
    hp, ap = predict_vs(trace, teams, i_ht, i_at)
    data = {
        'ホーム':ht,
        'アウェイ':at,
        'ホームの得点(期待値)':hp.mean(),
        'アウェイの得点(期待値)':ap.mean(),
        'ホームが勝つ確率(%)':100*sum(hp>ap)/len(hp),
        'ホームが負ける確率(%)':100*sum(hp<ap)/len(hp),
        'ホームが引分る確率(%)':100*sum(hp==ap)/len(hp),
    }
    vs_pred.append(pd.Series(data))
vs_pred = pd.concat(vs_pred, axis=1).T

vs_pred.head()

	ホーム	アウェイ	ホームの得点(期待値)	アウェイの得点(期待値)	ホームが勝つ確率(%)	ホームが負ける確率(%)	ホームが引分る確率(%)
0	名古屋	浦和	1.041	1.38	29.8	43.7	26.5
1	神戸	横浜FM	1.732	1.146	50.3	23.9	25.8
2	湘南	FC東京	1.135	1.243	35.2	39.6	25.2
3	鳥栖	Ｇ大阪	1.237	0.888	44.2	27	28.8
4	横浜FM	清水	1.74	1.493	43	32.9	24.1

以上

戦闘力のモデル化をpymcでやってみる。

pymc3のExampleには、ラグビーチームの攻撃力や防御力をモデル化する例がある。
https://docs.pymc.io/notebooks/rugby_analytics.html

こういったモデル化もあるのかと思っていると、世の中には、面白いことを考える人がいて、ドラゴンボールのキャラの戦闘力をモデル化した人がいる。

BUGS版（大元の記事）
http://statmodeling.hatenablog.com/entry/bugs-skill-performance

stan版
http://kazufusa1484.hatenablog.com/entry/2017/07/31/043000

BUGS版をベースにpymc3で実装してみた。 大体大元の記事と同じ結果が得られたが、モデルの一部コードを少し変えないと、MCMCがエラーで実行できなかったので、記事に書くことにした。

まず大元記事と同じデータを作る。

# キャラのスペック
data = {
    'CID': [1,2,3,4,5],
    'name':['ベジータ','悟空','ナッパ','クリリン','栽培マン'],
    'power.level':[18000,8000,4000,1770,1200]
}
senshi = pd.DataFrame.from_dict(data)

# 試合結果
data = {
    'winner':[1,1,1,2,3,4],
    'loser':[2,3,5,3,4,5]
}
battles = pd.DataFrame.from_dict(data)

BUGSのstep関数の代わりにswitch関数を使ってみたが、これだとmcmcがエラーで落ちたので、 sigmoid関数を使った。あとは同じ（と思う）。

with pm.Model() as model:
    # prior
    skill = pm.Normal('skill', mu=0, sd=100, shape=len(senshi.CID))
    
    winner_perf = pm.Normal('winner_perf', mu=skill[battles.winner-1], sd=1, shape=len(battles))
    loser_perf = pm.Normal('loser_perf', mu=skill[battles.loser-1], sd=1, shape=len(battles))
    
    diff = winner_perf-loser_perf
    #p = pm.math.switch(diff, 1.0, 0.0)
    p = pm.math.sigmoid(diff)
    z = pm.Bernoulli('z',p=p, observed=np.ones(len(battles)))

with model:
    trace = pm.sample(1000, start=pm.find_MAP(), step=pm.NUTS())

logp = -38.675, ||grad|| = 0.00058994: 100%|██████████| 84/84 [00:00<00:00, 643.83it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [loser_perf, winner_perf, skill]
Sampling 2 chains: 100%|██████████| 3000/3000 [03:08<00:00, 15.91draws/s]
The estimated number of effective samples is smaller than 200 for some parameters.

pm.traceplot(trace)

# 戦闘力部分の推定結果を抽出
result_df = pm.summary(trace, varnames=['skill'])

# パラメータ名をキャラ名に変換
senshi_dict = senshi.to_dict()['name']
senshi_dict2 = {}
for id,name in senshi_dict.items():
    senshi_dict2.setdefault('skill__%s'%(id), name)

result_df.rename(index=senshi_dict2, inplace=True)

真の戦闘力と推定戦闘力をviolinplotで表示する。

# violinplotに入れるデータの生成
data = []
for name in result_df.index:
    mu = result_df.loc[name, 'mean']
    sd = result_df.loc[name, 'sd']
    d = pm.Normal.dist(mu=mu,sd=sd).random(size=100)
    tmp = pd.DataFrame(data=d, columns=['power'])
    tmp['name'] = name
    ap = np.log10(senshi[senshi['name']==name]['power.level'])
    tmp['actual_power'] = np.round(ap.values[0], decimals=2)
    data.append(tmp)
data_df = pd.concat(data)

plt.figure(figsize=(12,9))
sns.violinplot(x='actual_power', y='power', data=data_df, hue='name')

栽培マンの幅が広いが、それ以外は大体同じ結果が得られた。

以上

中古マンション価格における駅・路線の影響分析をpymcでやってみる。

世の中には非常に優秀な人がいるもので、中古マンション価格について、そのスペックだけでなく、駅・路線の効果を分析した方がいる。（すごい！）

http://smrmkt.hatenablog.jp/entry/2014/11/17/234538

原著者の方は幾つかのモデルを試行錯誤的に検討されており、最終的なモデルは以下のようになっている。

こちらもstan+Rで実装されている。これは色々なデータに応用できそうなので、勉強のためにpymc3で実装した。 MCMCの実行においては、サンプリング手法としてNUTSを使うとともに、高速な近似手法であるADVI（自動微分変分法）を適用し、 計算速度と結果を比較した。

結果から言えば、MCMC（NUTS）による推定では、原著者とほぼ同じ推定結果が得られたが、4時間程度掛かった。 一方、ADIVは推定結果は少し（？）異なるが、2分程度で済んだ。従って、試行錯誤中はADVIで行い、 ある程度モデルが固まったらMCMCでしっかり推定する、と言う使い分けをしようと思う。

1. pymc3によるモデル実装

駅効果の部分が難しかったが、いろいろ試した結果、以下のような実装とした。

with pm.Model() as model4:
    # 家のスペック（広さ、駅からの距離など）の効果
    beta = pm.Normal('beta', mu=0, sd=100**2, shape=5)
    # 路線効果
    s_rt = pm.Uniform('s_rt', 0, 100**2)
    rt = pm.Normal('rt', mu=0, sd=s_rt, shape=len(df.train.unique()))
    # 駅単体の効果
    s_rs_base = pm.Uniform('s_rs_base', 0, 100**2)
    rs_base = pm.Normal('rs_base', mu=0, sd=s_rs_base, shape=len(df.station.unique()))
    # 駅は、駅単体と路線の影響を受ける
    s_rs = pm.Uniform('s_rs', 0, 100**2)
    rs = pm.Normal('rs', mu=rs_base+pm.math.dot(rt,data['ST'].T), sd=s_rs, shape=len(df.station.unique()))
    # マンションの平均価格
    price = beta[0] + \
            beta[1]*data['Distance'] + \
            beta[2]*data['From'] + \
            beta[3]*data['Room'] + \
            beta[4]*data['Space'] + \
            rs[data['S']-1]
    # 誤差
    eps = pm.HalfCauchy('eps', beta=100)
    #
    y_pred = pm.Normal('y_pred', mu=price, sd=eps, observed=data['Y'])

グラフィカルモデルは以下の通り。

pm.model_to_graphviz(model4)

2. MCMC（NUTS）による推定

データの読み込み部分から記載すると、以下のようになる。 私のMacbook Air(2011年モデル!)では、4時間ほど掛かった。

# データは原著者が公開しているgithubから。
# https://github.com/smrmkt/reins/tree/master/regression
df = pd.read_csv('mantions.csv')
df.dropna(inplace=True)

st = pd.read_csv('station_train.csv', header=None)

data = {
    'Distance': df.distance.values,
    'From': df['from'].values,
    'Room': df.room.values,
    'Space': df.space.values,
    'ST': st.values,
    'S': df.station.values,
    'T': df.train.values,
    'Y': df.price.values
}

# モデルは上述のため省略。

# MCMCの実行
with model4:
    trace4 = pm.sample(10000, tune=2000, start=pm.find_MAP(method='Powell'), step=pm.NUTS())

logp = -19,899, ||grad|| = 3,581.3: : 5001it [00:10, 489.30it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [eps, rs, s_rs, rs_base, s_rs_base, rt, s_rt, beta]
Sampling 2 chains: 100%|██████████| 24000/24000 [3:53:32<00:00,  1.98draws/s]

3. ADIVによる推定

ADVIは非常に高速。数分で終了。

with model4:
    inference = pm.ADVI()
    approx = pm.fit(n=60000, method=inference)
    trace4_advi = approx.sample(5000)

Average Loss = 18,925: 100%|██████████| 60000/60000 [02:27<00:00, 407.68it/s]
Finished [100%]: Average Loss = 18,925

4. 両者の比較

計算時間は実行パラメータに依るところもあるが、MCMC：4時間、ADVI：2分で、圧倒的にADVIが早い。

MCMCによる推定結果は以下の通り（路線、駅は上位・下位それぞれ5位表示）。 大体、原著者の結果と同じと思う。

print_effectiveness(df, trace4, upper_lower_station_n=5)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
切片	-1657.99	21.13	1.66e-01	-1700.18	-1617.23	12662.80	1.0
駅からの距離	-4.78	0.10	6.65e-04	-4.98	-4.59	20982.55	1.0
築年	0.85	0.01	8.71e-05	0.83	0.87	13096.16	1.0
部屋数	-0.42	0.19	1.51e-03	-0.79	-0.06	16784.44	1.0
床面積	1.93	0.15	1.25e-03	1.65	2.23	15127.33	1.0

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
横浜高速鉄道ＭＭ線	19.66	2.84	0.08	14.06	25.20	1079.64	1.0
東急田園都市線	10.25	2.02	0.07	6.41	14.33	667.58	1.0
東急東横線	9.33	1.81	0.06	5.80	12.89	748.61	1.0
京浜東北・根岸線	5.85	2.55	0.04	0.78	10.79	3999.36	1.0
東急目黒線	5.38	4.53	0.04	-3.62	14.21	10849.40	1.0
...	...	...	...	...	...	...	...
鶴見線	-1.55	3.23	0.05	-7.57	5.15	4981.11	1.0
京王電鉄相模原線	-1.88	3.15	0.06	-8.09	4.29	2647.75	1.0
小田急電鉄多摩線	-4.21	2.87	0.06	-9.85	1.37	2188.33	1.0
京浜急行電鉄大師線	-4.46	2.45	0.06	-9.13	0.52	1514.44	1.0
金沢シーサイドＬ	-7.45	2.28	0.06	-11.95	-2.94	1293.02	1.0

23 rows × 7 columns

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
みなとみらい	35.10	1.68	0.07	31.80	38.39	508.95	1.00
武蔵小杉	25.14	1.43	0.07	22.31	27.92	344.11	1.01
横浜	17.89	1.45	0.07	15.13	20.71	352.20	1.01
日本大通り	17.63	1.92	0.07	13.90	21.42	590.00	1.00
馬車道	17.57	2.24	0.07	13.27	22.05	807.19	1.00
...	...	...	...	...	...	...	...
産業道路	-8.70	3.83	0.07	-16.07	-1.02	2863.14	1.00
小島新田	-9.45	2.61	0.07	-14.62	-4.34	1173.79	1.00
海の公園柴口	-9.95	3.18	0.07	-16.04	-3.57	1940.90	1.00
八景島	-11.39	3.75	0.07	-18.58	-3.92	2534.88	1.00
南部市場	-12.09	3.76	0.07	-19.18	-4.49	2750.36	1.00

133 rows × 7 columns

ADVIの結果は以下。MCMCの推定結果と比較すると、ADVI目線で言うと、

• 家スペック効果では、切片が大きくことなる。また、床面積の効果が少し大きめ。
• 路線効果も順位の傾向は近いが、係数は大分小さい。家スペックの切片が小さいこととの兼ね合いと思われる。東海道本線が上位に来ている。
• 駅効果も上位は傾向が近いが、下位は少し違うようである。また、係数はこちらも小さい。

print_effectiveness(df,trace4_advi, upper_lower_station_n=5)

	mean	sd	mc_error	hpd_2.5	hpd_97.5
切片	-0.39	2.27e-01	3.31e-03	-0.83	4.35e-02
駅からの距離	-5.52	8.87e-02	1.25e-03	-5.69	-5.35e+00
築年	0.02	1.15e-04	1.46e-06	0.02	1.97e-02
部屋数	-0.17	8.60e-02	1.16e-03	-0.34	-1.59e-03
床面積	3.78	5.86e-02	8.64e-04	3.66	3.89e+00

	mean	sd	mc_error	hpd_2.5	hpd_97.5
横浜高速鉄道ＭＭ線	2.85	3.17	0.04	-3.34	9.15
東急東横線	1.46	2.23	0.03	-2.77	5.96
東急田園都市線	1.23	2.40	0.03	-3.42	5.98
東海道本線	1.19	3.61	0.05	-6.12	7.93
東急目黒線	0.91	3.85	0.05	-6.36	8.56
...	...	...	...	...	...
京浜急行電鉄逗子線	-1.33	3.62	0.05	-8.57	5.52
横浜市ブルーライン	-2.04	1.82	0.03	-5.71	1.37
小田急電鉄小田原線	-2.22	2.54	0.03	-7.02	3.01
金沢シーサイドＬ	-2.98	2.74	0.04	-8.48	2.21
京浜急行電鉄本線	-3.17	1.73	0.03	-6.53	0.22

23 rows × 5 columns

	mean	sd	mc_error	hpd_2.5	hpd_97.5
みなとみらい	13.35	3.67	0.05	6.30	20.54
武蔵小杉	12.73	2.16	0.03	8.51	16.98
日本大通り	11.07	3.30	0.05	4.61	17.39
横浜	10.63	1.40	0.02	7.92	13.24
たまプラーザ	10.46	2.26	0.03	6.36	15.10
...	...	...	...	...	...
京急富岡	-10.50	3.89	0.06	-18.09	-3.08
磯子	-10.76	2.27	0.03	-15.08	-6.22
並木北	-11.06	4.93	0.07	-21.17	-1.96
六浦	-11.13	4.24	0.07	-19.72	-3.09
並木中央	-11.40	3.92	0.06	-19.53	-4.00

133 rows × 5 columns

ちなみに、結果を表示している関数は以下。

# MCMCの結果を表示する際、パラメータ変数名ではなく、路線・駅名で表示するための関数。
def print_effectiveness(df, trace, upper_lower_station_n=None):
    # パラメータと家効果の対応表を作成する。
    beta_dict = {
        'beta__0':'切片',
        'beta__1':'駅からの距離',
        'beta__2':'築年',
        'beta__3':'部屋数',
        'beta__4':'床面積'
    }
    # パラメータと路線の対応表を作成する。
    train_dict = (df[['train','train_raw']]).drop_duplicates().set_index('train').to_dict()['train_raw']
    train_rename_dict = {}
    for key, name in train_dict.items():
        train_rename_dict.setdefault('rt__%d'%(key-1))
        train_rename_dict['rt__%d'%(key-1)] = name
    # パラメータと駅の対応表を作成する。
    station_dict = (df[['station','station_raw']]).drop_duplicates().set_index('station').to_dict()['station_raw']
    station_rename_dict = {}
    for key, name in station_dict.items():
        station_rename_dict.setdefault('rs__%d'%(key-1))
        station_rename_dict['rs__%d'%(key-1)] = name
    # traceのパラメータを家スペック・駅・路線名に置き換える
    ret_h = pm.summary(trace, varnames=['beta'])
    ret_s = pm.summary(trace, varnames=['rs'])
    ret_t = pm.summary(trace, varnames=['rt'])
    ret_h.rename(index=beta_dict, inplace=True)
    ret_t.rename(index=train_rename_dict, inplace=True)
    ret_s.rename(index=station_rename_dict, inplace=True)
    # 家スペック効果を表示する
    display(ret_h)
    # 路線効果、駅効果を表示する。
    if upper_lower_station_n is None:
        # 全部表示
        display(ret_t.sort_values('mean',ascending=False))
        display(ret_s.sort_values('mean',ascending=False))
    else:
        # 上位,下位を表示する。
        max_rows = pd.get_option('display.max_rows')
        pd.set_option('display.max_rows', upper_lower_station_n*2)
        display(ret_t.sort_values('mean',ascending=False))
        display(ret_s.sort_values('mean',ascending=False))
        pd.set_option('display.max_rows', max_rows)

5. 結論

トップでも書いたが、ADVIは圧倒的に高速だが、やはり近似であるため、MCMCの結果とは多少異なる。 ただ、傾向はつかめると思うので、検討初期の試行錯誤中はADVIを使い、ある程度モデルが固まってきたらMCMCでしっかり推定する、 というのが良いと思われる。

以上

MCMCにおける標準化処理の効果について

前回紹介した「ベイズ統計モデリング」では、MCMCの効率的な実行のため、データの前処理として、標準化または平均化が推奨されている。

ベイズ統計モデリング R、JAGS、 Stanによるチュートリアル　原著第2版 [ John K. Kruschke ]

今回は、簡単なGLMを例に、その効果（サンプリング時間）を計測し、比較した。

1. 前処理なし
2. 標準化処理あり
3. 中心化処理あり
4. pm.glmモジュール利用（前処理なし）

例題は、pymc3公式HPの例題を使用した。
GLM: Linear regression — PyMC3 3.5 documentation

結果は、いずれもtraceplotやRhatを見ると収束しており、 推定されたパラメータは、それほど大きな差はない。 グラフに表示してもほとんど分からない。（例題にも依ると思うが。）

しかし、MCMCのサンプリングに掛かった時間は、
「標準化あり < 中心化あり < 前処理なし < pm.glmモジュール」
であり、前処理ありの方が早い。また、有効サンプルサイズも前処理ありの方が2,3倍大きい。

	前処理なし	標準化	中心化	pm.glm
実行時間	24.5 s	10.5 s	11.4 s	37.5 s
有効サンプルサイズ(beta)	2523.87	8809.77	8567.02	2916.77
Rhat	1.0	1.0	1.0	1.0

従って、複雑なモデルのパラメータ推定には前処理は有効であると言える。 しかし、前処理されているため、推定されたパラーメータ値が直感的には解釈し難いのが欠点と言える。

中心化は並行移動だけで、標準化と異なりスケールは変わらないため、解釈しやすい。 効果としても標準化を行った場合と大差ないため、うまくMCMCが実行されない場合は 中心化だけでもやった方が良いかもしれない。

今回実行したコードを以下に載せておく。

関数定義

今回はモデルやMCMCパラメータは同じで、データの前処理の有無の差だけなので、関数定義した。

def exec_mcmc(x,y):
    with pm.Model() as model:
        sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)
        alpha = pm.Normal('alpha', 0, sd=20)
        beta = pm.Normal('beta', 0, sd=20)

        y_ = pm.Normal('y', mu=alpha + beta * x, sd=sigma, observed=y)

        trace = pm.sample(3000, tune=500, start=pm.find_MAP(method='Powell'), step=pm.NUTS())
    return model, trace

def plot_figure(x, y, true_regression_line, pred_y=None, label=None):
    fig = plt.figure(figsize=(7, 7))
    ax = fig.add_subplot(111, xlabel='x', ylabel='y', title='Generated data and underlying model')
    ax.plot(x, y, 'x', label='sampled data')
    # 真の回帰線を表示する
    ax.plot(x, true_regression_line, label='true regression line', lw=2.)
    if pred_y is not None:
        # 予測した回帰線を表示する
        if type(pred_y) is list:
            for py,l in zip(pred_y, label):
                ax.plot(x, py, label=l, lw=2.)
        else:
            ax.plot(x, pred_y, label=label, lw=2.)
    plt.legend(loc=0);

データの生成

# データの生成
size = 200
true_intercept = 1
true_slope = 2

x = np.linspace(0, 100, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * x
# add noise
y = true_regression_line + np.random.normal(scale=50.5, size=size)

data = dict(x=x, y=y)
# DataFrame化
df = pd.DataFrame(data=data)

1. 前処理なし

%%time
model1, trace1 = exec_mcmc(df.x,df.y)

logp = -1,077.5, ||grad|| = 5.3893: 100%|██████████| 78/78 [00:00<00:00, 812.06it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:11<00:00, 594.10draws/s]

CPU times: user 8.8 s, sys: 660 ms, total: 9.46 s
Wall time: 24.5 s

2, 標準化

# 標準化処理
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
trans_df = pd.DataFrame(data=scaler.fit_transform(df.values), columns=['x','y'])

%%time
model2, trace2 = exec_mcmc(trans_df.x, trans_df.y)

logp = -207.47, ||grad|| = 0.0020073: 100%|██████████| 83/83 [00:00<00:00, 851.20it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:06<00:00, 1049.27draws/s]


CPU times: user 6.54 s, sys: 304 ms, total: 6.84 s
Wall time: 10.5 s

3. 中心化

# 中心化処理
scaler2 = StandardScaler(with_std=False)
trans_df2 = pd.DataFrame(data=scaler2.fit_transform(df.values), columns=['x','y'])

%%time
model3, trace3 = exec_mcmc(trans_df2.x, trans_df2.y)

logp = -1,077.3, ||grad|| = 0.07435: 100%|██████████| 139/139 [00:00<00:00, 1004.37it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:07<00:00, 915.98draws/s] 

CPU times: user 6.44 s, sys: 240 ms, total: 6.68 s
Wall time: 11.4 s

4. pm.glmモジュールの利用

%%time
with pm.Model() as model4:
    pm.glm.GLM.from_formula('y ~ x', df)
    trace4 = pm.sample(3000, tune=500, start=pm.find_MAP(method='Powell'), step=pm.NUTS())

logp = -1,077.3, ||grad|| = 0.011737: 100%|██████████| 224/224 [00:00<00:00, 620.51it/s] 
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [sd, x, Intercept]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:31<00:00, 222.34draws/s]


CPU times: user 11.9 s, sys: 628 ms, total: 12.6 s
Wall time: 37.5 s

結果の表示

# 前処理なしの結果
pm.summary(trace1)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
alpha	-5.00	6.59	1.40e-01	-17.11	8.61	2686.89	1.0
beta	1.99	0.12	2.41e-03	1.77	2.22	2523.87	1.0
sigma	49.64	2.55	4.61e-02	44.76	54.63	3204.49	1.0

# 標準化ありの結果
pm.summary(trace2)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
alpha	2.45e-04	0.05	4.48e-04	-0.09	0.10	9703.56	1.0
beta	7.62e-01	0.05	5.09e-04	0.67	0.85	8809.77	1.0
sigma	6.55e-01	0.03	3.18e-04	0.59	0.72	8292.82	1.0

# 中心化ありの結果
pm.summary(trace3)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
alpha	0.01	3.48	3.49e-02	-6.79	6.75	8898.01	1.0
beta	2.00	0.12	1.47e-03	1.75	2.23	8567.02	1.0
sigma	49.65	2.49	2.94e-02	45.26	54.97	7896.93	1.0

# pm.glmモジュールの結果
pm.summary(trace4)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
Intercept	-5.72	6.89	1.30e-01	-18.82	7.98	2894.93	1.0
x	2.00	0.12	2.27e-03	1.75	2.22	2916.77	1.0
sd	49.66	2.51	3.97e-02	44.89	54.55	4223.24	1.0

グラフ表示

パラメータの推定結果は、大体一緒。 傾きはうまく推定されているが、切片がノイズが大きいためか、少しズレている。

pred_y_ = []
for tr,d,ss in zip([trace1, trace2, trace3, trace4],
                   [df,trans_df,trans_df2,df],
                   [None,scaler,scaler2,None]):
    try:
        tmp = np.array([a*tr['beta'].mean()+tr['alpha'].mean() for a in d.x])
    except:
        # pm.glmモジュール利用の場合はこっちになる
        tmp = np.array([a*tr['x'].mean()+tr['Intercept'].mean() for a in d.x])
    if ss is not None:
        tmp = ss.inverse_transform(np.array([d.x.values, np.array(tmp)]).T).T[1]
    pred_y_.append(tmp)

plot_figure(x,y,true_regression_line,pred_y=pred_y_,label=['なし','標準化','中心化','glm'])

以上

Kruschkeのモデルダイアグラム

モデルの表記形式

ベイズ推定におけるモデルの表記形式としては、有向グラフが最も一般的であろうと思う。 一方、「ベイズ統計モデリング」の原著者であるKruschke先生は、新しい表記形式でモデルを記載している。

ベイズ統計モデリング R、JAGS、 Stanによるチュートリアル　原著第2版 [ John K. Kruschke ]

800頁近いボリュームと9000円程度の出費と、購入には中々勇気が必要だったが、 内容はベイズの公式やMCMCの仕組みなど基礎的な話から始まり、仮設検定のベイズアプローチ、階層ベイズ、GLMなど 応用に入っていく流れで、非常に勉強になった。費用と時間を掛けた価値は十分にあったと思う。 これからベイズを学ぶ人には、個人的にはお勧め。

このダイアグラムは、その理解のしやすさから、他の書籍（Kruschke先生以外の人の書籍）でも使われ始めているようで、 「Kruschkeのダイアグラム」と呼ばれている。

Kruschkeのダイアグラム

Kruschkeのダイアグラム（以降、単にダイアグラムと記載）は、書籍や先生のブログを見るのが一番早い。 モデルにおいて、どのような確率分布が組み合わされてデータが生成されると仮定しているか、 視覚的に非常に分かりやすい表記になっている。

doingbayesiandataanalysis.blogspot.com

モデルの有向グラフは、pymc3の3.5を使っていれば、pm.model_to_graphviz(model)を使えば、 コードからgraphvizを使って自動的に生成される。（例えば、以前の記事参照）

gaiasky.hatenablog.com

ダイアグラムはどうやって書くのが良いかと調べていたところ、LibreOfficeのDrawで、テンプレートを作って 公開してくれている人がいた。

www.sumsar.net

LibreOfficeはマイクロソフトOfficeに近い（これらのファイルも作成・編集可能）ツールで、 オープンソースなので無料で使える。（Macでも、Linuxでも使える）

Home | LibreOffice - Free Office Suite - Fun Project - Fantastic People

有向グラフは、pythonコードから自動生成できるので楽ちんだが、 複雑なモデルの場合はダイアグラムの方が分かりやすいと思うので、 複雑なモデルはこちらを使おうと思う。

空間統計モデル（CAR model）をpymc3で。（緑本11章）

CAR model（条件つき自己回帰モデル）

いわゆる久保先生の緑本（データ解析のための統計モデリング入門）11章で紹介されている空間構造を考慮したモデル。

データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC （シリーズ確率と情報の科学） [ 久保拓弥 ]

ベイズ統計モデルに興味を持ったきっかけの本。正規分布を仮定する線形モデル（LM）から始まり、正規分布に限定されない一般化線形モデル（GLM）、個体差を考慮する一般化線形混合モデル（GLMM）と発展し、複雑なモデルのパラメータを解析的に求めることの難しさから、ベイズ＋MCMCによる推定への話がつながっており非常に分かりやすかった。

書籍では、モデルはBUGSで記載されており、Webを検索すればstanコードもたくさん出てくる。一方、pythonコードは少なく、特に11章の空間統計モデル（CAR model）の実装は見当たらなかったため、pymc3版を書くことにした。

BUGSではCARモデルの関数が用意されているようだが、pymc3ではそういった関数はない模様。一方で、pymc3の公式HPにはexampleがあったため、これを参考にした。

pymc3によるCARモデル

Intrinsic CARモデルは下式で表される。 wは重みで、場所iにおける場所jから受ける影響の大きさを表す。 例えば、場所iとjが隣接している場合、w_ij=1、隣接していない場合、w_ij=0などとする。 この場合、平均や分散の分母Σw_ijは隣接する場所の数となる。

実装方法としては、公式HPに、CARモデルの例題が紹介されている。①theano.scanを使う方法、②matrix trickを使う方法などが紹介されているが、①は実行速度が遅いと書かれているため、②を使うことにした。②の方法では、CAR2クラスという独自の分布モデルを作成する（ちなみに、CARクラスも例には載っており、これは①に対応したクラス）。パラメータのaが隣接情報（何番目の場所が隣接するか）の行列で、wが重み情報の行列になる。

PyMC3 Modeling tips and heuristic — PyMC3 3.5 documentation

class CAR2(distribution.Continuous):
    """
    Conditional Autoregressive (CAR) distribution

    Parameters
    ----------
    a : adjacency matrix
    w : weight matrix
    tau : precision at each location
    """

    def __init__(self, w, a, tau, *args, **kwargs):
        super(CAR2, self).__init__(*args, **kwargs)
        self.a = a = tt.as_tensor_variable(a)
        self.w = w = tt.as_tensor_variable(w)
        self.tau = tau*tt.sum(w, axis=1)
        self.mode = 0.

    def logp(self, x):
        tau = self.tau
        w = self.w
        a = self.a

        mu_w = tt.sum(x*a, axis=1)/tt.sum(w, axis=1)
        return tt.sum(continuous.Normal.dist(mu=mu_w, tau=tau).logp(x))

まず、データを用意する。

Y = np.array([  0.,   3.,   2.,   5.,   6.,  16.,   8.,  14.,  11.,  10.,  17.,
                19.,  14.,  19.,  19.,  18.,  15.,  13.,  13.,   9.,  11.,  15.,
                18.,  12.,  11.,  17.,  14.,  16.,  15.,   9.,   6.,  15.,  10.,
                11.,  14.,   7.,  14.,  14.,  13.,  17.,   8.,   7.,  10.,   4.,
                 5.,   5.,   7.,   4.,   3.,   1.], np.int64)

次に、隣接情報と重み情報を作成する。これは、自身の両隣で、重みはどちらも同じ大きさ。

# 隣接情報
adj = np.array(
       [[1], [0, 2], [1, 3], [2, 4], [3, 5], [4, 6], [5, 7], [6, 8], [7, 9],
       [8, 10], [9, 11], [10, 12], [11, 13], [12, 14], [13, 15], [14, 16],
       [15, 17], [16, 18], [17, 19], [18, 20], [19, 21], [20, 22],
       [21, 23], [22, 24], [23, 25], [24, 26], [25, 27], [26, 28],
       [27, 29], [28, 30], [29, 31], [30, 32], [31, 33], [32, 34],
       [33, 35], [34, 36], [35, 37], [36, 38], [37, 39], [38, 40],
       [39, 41], [40, 42], [41, 43], [42, 44], [43, 45], [44, 46],
       [45, 47], [46, 48], [47, 49], [49]], dtype=object)

# 重み
weights = np.array(
       [[1.0], [1, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0],
       [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0]], dtype=object)

これらを行列化する。

N = len(weights)
wmat2 = np.zeros((N,N))
amat2 = np.zeros((N,N))
for i, a in enumerate(adj):
    amat2[i,a] = 1
    wmat2[i,a] = weights[i]

wmat2

array([[0., 1., 0., ..., 0., 0., 0.],
       [1., 0., 1., ..., 0., 0., 0.],
       [0., 1., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 1., 0.],
       [0., 0., 0., ..., 1., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.]])

amat2

array([[0., 1., 0., ..., 0., 0., 0.],
       [1., 0., 1., ..., 0., 0., 0.],
       [0., 1., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 1., 0.],
       [0., 0., 0., ..., 1., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.]])

モデルの作成。

with pm.Model() as model:
    # hyper
    s = pm.Uniform('s', lower=0, upper=100)
    # prior
    beta = pm.Normal('beta', mu=0, sd=10)
    tau = 1/(s*s)
    r =CAR2('r', w=wmat2, a=amat2, tau=tau, shape=N)
    #
    y = pm.Poisson('y', mu=np.exp(beta + r[np.arange(50)]), observed=Y)

モデルのグラフィカル表示。

pm.model_to_graphviz(model)

MCMCの実行。NUTSでサンプリングすると、収束に時間が掛かったため、ADVI（自動微分変分法）を使ってみた。

%%time
with model:
    inference = pm.ADVI()
    approx = pm.fit(n=100000, method=inference,random_seed=123, start=pm.find_MAP(method='Powell'))
    trace = approx.sample(draws=5000)

logp = -29.053, ||grad|| = 50.912: : 5001it [00:05, 926.82it/s]                          

Average Loss = 143.95: 100%|██████████| 100000/100000 [02:09<00:00, 771.11it/s]
Finished [100%]: Average Loss = 143.95


CPU times: user 2min 12s, sys: 20.1 s, total: 2min 32s
Wall time: 2min 29s

結果の出力。

pm.traceplot(trace)

データと推定結果をグラフ表示する。青丸がデータ。黒線が推定結果（中央値）、グレー帯が95%HPD区間。

b = np.median(trace.get_values('beta'))
r = np.median(trace.get_values('r'),axis=0)
xx = np.arange(N)
yy = np.exp(b+r)
plt.plot(xx,yy, 'k-')
hpd = pm.hpd(trace.get_values('r'), alpha=0.05)
plt.fill_between(xx, np.exp(b+hpd.T[0]), np.exp(b+hpd.T[1]), color='k', alpha=0.2)

plt.plot(range(len(Y)),Y, 'p')
plt.xlabel('position $j$')
plt.ylabel('number of individuals $y_j$')
plt.ylim(-2.5,30)

推定したパラメータの詳細は下表。 betaもsも書籍とは少し異なった推定値が得られているが、 上図のグラフを見ると、そこそこ良い推定結果が得られていると思う。

pm.summary(trace)

	mean	sd	mc_error	hpd_2.5	hpd_97.5
beta	1.158679	0.047799	0.000675	1.068305	1.254242
r__0	-0.888382	0.301022	0.004542	-1.492071	-0.318896
r__1	-0.653750	0.205648	0.003100	-1.063285	-0.269927
r__2	-0.218527	0.211978	0.003377	-0.651866	0.181730
r__3	0.315092	0.205860	0.003063	-0.103261	0.703068
r__4	0.832849	0.195722	0.002816	0.461069	1.226940
r__5	1.281606	0.176526	0.002431	0.932831	1.627272
r__6	1.260697	0.181086	0.002758	0.902225	1.618805
r__7	1.321297	0.183583	0.002787	0.948448	1.661514
r__8	1.264942	0.180075	0.002479	0.901088	1.611691
r__9	1.315550	0.178816	0.002678	0.956416	1.655632
r__10	1.570310	0.170607	0.002199	1.230681	1.900745
r__11	1.695244	0.165743	0.002195	1.351274	2.006621
r__12	1.612901	0.170895	0.002235	1.284063	1.965944
r__13	1.727920	0.164937	0.002407	1.416415	2.059092
r__14	1.752215	0.168899	0.002160	1.412582	2.071013
r__15	1.705700	0.166617	0.002456	1.397315	2.053247
r__16	1.547344	0.172276	0.002583	1.210396	1.883085
r__17	1.415188	0.177532	0.002456	1.081155	1.770285
r__18	1.296239	0.186362	0.002506	0.938413	1.665728
r__19	1.165170	0.184071	0.002463	0.817927	1.541284
r__20	1.274738	0.179306	0.002607	0.923394	1.621070
r__21	1.502048	0.172138	0.002584	1.168851	1.827205
r__22	1.587426	0.170370	0.002340	1.247124	1.908877
r__23	1.415244	0.171736	0.002149	1.096096	1.761579
r__24	1.373549	0.178739	0.002443	1.022405	1.719808
r__25	1.541359	0.171126	0.002365	1.187936	1.854710
r__26	1.557773	0.173536	0.002299	1.220501	1.894539
r__27	1.557409	0.176749	0.002737	1.210693	1.904542
r__28	1.417466	0.176512	0.002367	1.068018	1.756552
r__29	1.136758	0.184823	0.002607	0.789616	1.508527
r__30	1.020712	0.185264	0.002552	0.671178	1.397461
r__31	1.254099	0.176998	0.002652	0.921095	1.610488
r__32	1.247103	0.179981	0.002625	0.918825	1.624593
r__33	1.256230	0.179677	0.002617	0.903487	1.614004
r__34	1.286939	0.178641	0.002481	0.924141	1.616439
r__35	1.160147	0.188822	0.002600	0.771362	1.514623
r__36	1.340738	0.171982	0.002378	0.998691	1.665439
r__37	1.451340	0.173466	0.002303	1.102961	1.783128
r__38	1.484233	0.174277	0.002155	1.135390	1.825727
r__39	1.461851	0.171625	0.002388	1.121778	1.793317
r__40	1.129038	0.183425	0.002540	0.745838	1.467143
r__41	0.915550	0.190102	0.002935	0.541487	1.276502
r__42	0.825596	0.197854	0.002697	0.431916	1.193682
r__43	0.568216	0.201373	0.002649	0.165221	0.949261
r__44	0.480365	0.195278	0.002611	0.092885	0.859808
r__45	0.523366	0.192335	0.002415	0.158896	0.904041
r__46	0.526652	0.199353	0.002596	0.126299	0.898693
r__47	0.250715	0.199579	0.002248	-0.161358	0.618004
r__48	-0.212703	0.229811	0.003280	-0.637011	0.258963
r__49	-0.801395	0.463723	0.007071	-1.727559	0.096066
s	0.347891	0.040211	0.000537	0.269606	0.424428

以上