アイドルの人気投票における多項ロジスティック回帰の適用例
以前、アヒル本の多項ロジスティック回帰についてpymc3版を紹介した。
今回は、多項ロジスティック回帰の例として、「μ's とAqours の人気の差」を題材とした記事があったので、これを紹介したいと思う。 これらの記事ではモデルはStanで実装されていたので、これをpymc3でトレースしてみることにする。
http://mikuhatsune.hatenadiary.com/entry/20170320/1490011326
http://abrahamcow.hatenablog.com/entry/2017/06/24/154902
データ生成のモデルは以下の通り
- 各9人が所属する2つのグループがある。
- メンバーの効果（人気）を
、μ'sグループの効果（人気）を
とする。
- メンバー
が所属するグループを
とする。
はμ's、
はAqousに所属。
- メンバーの効果（人気）を
- 各メンバーの投票数
とする。
- 得票数は多項分布に従う。
- 得票数は多項分布に従う。
まずはデータの作成
dat_str = """name vote group 矢澤にこ 384 muse 西木野真姫 384 muse 南ことり 370 muse 東條希 336 muse 綾瀬絵里 283 muse 園田海未 263 muse 小泉花陽 251 muse 星空凛 226 muse 津島善子 190 aqours 黒澤ルビィ 171 aqours 国木田花丸 155 aqours 渡辺曜 141 aqours 高坂穂乃果 128 muse 松浦果南 110 aqours 黒澤ダイヤ 103 aqours 桜内梨子 93 aqours 小原鞠莉 80 aqours 高海千歌 74 aqours"""
df = pd.read_csv(io.StringIO(dat_str), sep=' ')
名前とグループをコード化する。
name_cate = pd.Categorical(df.name) group_cate = pd.Categorical(df.group) df['name_code'] = name_cate.codes df['group_code'] = group_cate.codes df.sort_values('name_code', inplace=True) # name_code順に並び替える
df.head()
| name | vote | group | name_code | group_code | |
|---|---|---|---|---|---|
| 2 | 南ことり | 370 | muse | 0 | 1 |
| 10 | 国木田花丸 | 155 | aqours | 1 | 0 |
| 5 | 園田海未 | 263 | muse | 2 | 1 |
| 16 | 小原鞠莉 | 80 | aqours | 3 | 0 |
| 6 | 小泉花陽 | 251 | muse | 4 | 1 |
plt.figure(figsize=(9,6)) sns.barplot(x='vote', y='name', data=df, hue='group', orient='h')
pymc3によるモデルの実装とMCMC実行
with pm.Model() as model: # alpha[i] = g*G[i] + b[i] # 各個人の効果 b_mu = pm.Normal('b_mu', mu=0, sd=5) b_sd = pm.HalfCauchy('b_sd', beta=5) b = pm.Normal('b', mu=b_mu, sd=b_sd, shape=len(df)-1) b_f = pm.Deterministic('b_f', tt.concatenate([tt.zeros(1), b])) # 一人目の効果を0で固定 # グループ効果 g = pm.Normal('g', mu=0, sd=10) # 合計効果 alpha = pm.Deterministic('alpha', g*df.group_code+b_f[df.name_code]) # p = pm.Deterministic('p', tt.nnet.softmax(alpha)) y = pm.Multinomial('y', df.vote.sum(), p, observed=df.vote)
with model: trace = pm.sample(5000)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [g, b, b_sd, b_mu]
Sampling 2 chains: 100%|██████████| 11000/11000 [01:05<00:00, 168.16draws/s]
The number of effective samples is smaller than 25% for some parameters.
pm.traceplot(trace, varnames=['g','b','p'])
推定結果の抽出
グループ（μ's）効果の確認。μ'sは、投票にに対してAquosより2.3倍の影響力がある。
pm.plot_posterior(np.exp(trace.get_values('g'))) #muse効果
可視化用のコード。
# 推定結果をグラフ表示するための整形関数 def mk_data4graph(df, trace, vname, do_exp): skip = 20 if do_exp: estimate_params = np.exp(pm.trace_to_dataframe(trace, varnames=[vname]).loc[::skip]).T else: estimate_params = pm.trace_to_dataframe(trace, varnames=[vname]).loc[::skip].T estimate_params.index = name_cate.categories # 平均値が高いキャラでソートする estimate_params['mean'] = estimate_params.mean(axis=1) estimate_params_sorted = estimate_params.sort_values('mean', ascending=False) estimate_params_sorted = estimate_params_sorted.drop(columns='mean') # グラフ用データに変換し、グループ名を追加 data4graph = pd.merge(pd.melt(estimate_params_sorted.T), df[['name','group']], left_on='variable', right_on='name') return data4graph
各メンバーの投票率の可視化。
# 各メンバーの投票率） plt.figure(figsize=(6,9)) sns.violinplot(x='value', y='variable', data=mk_data4graph(df, trace,'p', False), hue='group', palette=sns.color_palette("Set1", 2), orient='h',) plt.xlabel('投票率$p_i$') plt.title('メンバーの投票率$p_i$の推定値')
各メンバーの人気の可視化。
# 各メンバーの人気（exp(b)） plt.figure(figsize=(6,9)) sns.violinplot(x='value', y='variable', data=mk_data4graph(df, trace,'b_f', True), hue='group', palette=sns.color_palette("Set1", 2)[::-1], orient='h') plt.xlabel('人気$\exp(b)$') plt.title('メンバーの人気$\exp(b_i)$の推定値')
モデルを定式化でき、影響の大きさが定量的に把握できるのが、ロジスティック回帰はやはり有用だと思う。
以上
ベイズ生存時間分析（Weibull分布）
データ分析では正規分布を仮定することが多いが、生存時間分析・信頼性工学では、ワイブル分布を仮定することが多い。これはワイブル分布が、形状パラメータ・尺度パラメータによって、所謂バスタブカーブの３要素（初期故障、偶発故障、摩耗）を表現可能であるからと思う。一方、データには打ち切りが存在することが一般的で、そのために、ノンパラメトリックモデルのカプランマイヤー法で推定することが多いが、ベイズ推定を用いるとパラメトリックモデルで打ち切りデータもうまく扱えるようだ。
Bayesian Inference With Stan ~062~ | |
上記の記事はstanで実装されているため、今回はこの例題をpymcで実装してみることにした。
# データセットの読み込み。
df = load_rossi()
rossiデータの概要。
- 再逮捕までの時間（週単位）
- 右側打ち切り（52週）
- 共変量あり（fin=1, 経済的支援あり）
df.head()
| week | arrest | fin | age | race | wexp | mar | paro | prio | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 20 | 1 | 0 | 27 | 1 | 0 | 0 | 1 | 3 |
| 1 | 17 | 1 | 0 | 18 | 1 | 0 | 0 | 1 | 8 |
| 2 | 25 | 1 | 0 | 19 | 0 | 1 | 0 | 1 | 13 |
| 3 | 52 | 0 | 1 | 23 | 1 | 1 | 1 | 1 | 1 |
| 4 | 52 | 0 | 0 | 19 | 0 | 1 | 0 | 1 | 3 |
df_obs = df[df.arrest==1].week.values # 打ち切りなしのデータ df_fin_obs = df[df.arrest==1].fin.values df_cens = df[df.arrest==0].week.values # 右側打ち切りのデータ df_fin_cens = df[df.arrest==0].fin.values
ワイブル分布の関数定義
# pdf def weibProbDist(x, a, b): return (a / b) * (x / b) ** (a - 1) * np.exp(-(x / b) ** a) # cdf def weibCumDist(x, a, b): return 1 - np.exp(-(x / b) ** a)
ワイブル分布への当てはめ
pymc3の公式HPのExampleを参考にした。
Reparameterizing the Weibull Accelerated Failure Time Model — PyMC3 3.5 documentation
打ち切りなしデータと、打ち切りありデータのそれぞれを分けてサンプリングする。
def weibull_lccdf(x, shape, scale): ''' Log complementary cdf of Weibull distribution. ''' return -(x / scale)**shape
with pm.Model() as model_1: # 形状パラメータ alpha_sd = 10.0 alpha_raw = pm.Normal('a0', mu=0, sd=0.1) shape = pm.Deterministic('shape', tt.exp(alpha_sd * alpha_raw)) # 尺度パラメータ mu = pm.Normal('mu', mu=0, sd=100) scale = pm.Deterministic('scale', tt.exp(mu / shape)) y_obs = pm.Weibull('y_obs', alpha=shape, beta=scale, observed=df_obs) y_cens = pm.Potential('y_cens', weibull_lccdf(df_cens, shape, scale))
with model_1: trace_1 = pm.sample(5000, tune=2000, nuts_kwargs={'target_accept': 0.95}, init='adapt_diag')
Auto-assigning NUTS sampler...
Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [mu, a0]
Sampling 2 chains: 100%|██████████| 14000/14000 [00:50<00:00, 277.78draws/s]
The number of effective samples is smaller than 25% for some parameters.
pm.traceplot(trace_1)
pm.summary(trace_1).round(2)
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| a0 | 0.03 | 0.01 | 0.00 | 0.01 | 0.05 | 1680.58 | 1.0 |
| mu | 6.56 | 0.49 | 0.01 | 5.60 | 7.50 | 1671.80 | 1.0 |
| shape | 1.36 | 0.12 | 0.00 | 1.12 | 1.60 | 1671.94 | 1.0 |
| scale | 126.44 | 14.47 | 0.26 | 101.34 | 156.59 | 2817.21 | 1.0 |
推定したパラメータによる生存時間関数をグラフ表示する。
qlist = [2.5, 50, 97.5] shapes = pm.quantiles(trace_1.get_values('shape'), qlist=qlist) scales = pm.quantiles(trace_1.get_values('scale'), qlist=qlist)
# 生存時間関数表示関数 def plot_survival_function(ax, x, shape, scale, color='blue', qlist=[2.5,50,97.5]): low=qlist[0] mid=qlist[1] hig=qlist[2] ax.plot(x, 1-weibCumDist(x, shape[mid], scale[mid]), color=color) ax.fill_between(x=x, y1=1-weibCumDist(x, shape[low], scale[low]), y2=1-weibCumDist(x, shape[hig], scale[hig]), alpha=0.2, color=color) ax.set_xlabel('週') ax.set_ylabel('生存率') return ax
x = np.linspace(0,52,100) fig = plt.figure(figsize=(9,4)) ax = fig.add_subplot(111) ax = plot_survival_function(ax, x, shapes, scales, qlist=qlist)
カプランマイヤー法
- ノンパラメトリックモデルであるカプランマイヤー法による生存率推定結果と比較する。
from lifelines import KaplanMeierFitter kmf = KaplanMeierFitter() kmf.fit(df.week.values, df.arrest.astype(np.bool), alpha=0.975)
x = np.linspace(0,52,100) fig = plt.figure(figsize=(9,4)) ax = fig.add_subplot(111) ax = plot_survival_function(ax, x, shapes, scales, qlist=qlist) kmf.plot(ax=ax)
ベイズ推定はカプランマイヤー法よりも広めに推定されている。 元の記事の推定結果よりも広いため、うまく推定できていないのかも知れない。
比例ハザードモデル
- 形状パラメータではなく、尺度パラメータに共変量の効果を入れる模様。
with pm.Model() as model_2: # 形状パラメータ alpha_sd = 10.0 alpha_raw = pm.Normal('a0', mu=0, sd=0.1) shape = pm.Deterministic('shape', tt.exp(alpha_sd * alpha_raw)) # 尺度パラメータ beta = pm.Normal('beta', mu=0, sd=10**2, shape=2) # 経済的支援ありの場合の尺度パラメータ scale1 = pm.Deterministic('scale1', tt.exp(-(beta[0]+beta[1])/shape)) # 経済的支援なしの場合の尺度パラメータ scale2 = pm.Deterministic('scale2', tt.exp(-(beta[0])/shape)) y_obs = pm.Weibull('y_obs', alpha=shape, beta=tt.exp(-(beta[0]+df_fin_obs*beta[1])/shape), observed=df_obs) y_cens = pm.Potential('y_cens', weibull_lccdf(df_cens, shape, tt.exp(-(beta[0]+df_fin_cens*beta[1])/shape)))
with model_2: trace_2 = pm.sample(5000, tune=2000, nuts_kwargs={'target_accept': 0.9}, init='adapt_diag')
Auto-assigning NUTS sampler...
Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, a0]
Sampling 2 chains: 100%|██████████| 14000/14000 [01:30<00:00, 154.14draws/s]
pm.summary(trace_2).round(2)
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| a0 | 0.03 | 0.01 | 0.00 | 0.01 | 0.05 | 2917.32 | 1.0 |
| beta__0 | -6.42 | 0.49 | 0.01 | -7.38 | -5.49 | 2859.43 | 1.0 |
| beta__1 | -0.37 | 0.19 | 0.00 | -0.73 | 0.01 | 4608.48 | 1.0 |
| shape | 1.37 | 0.12 | 0.00 | 1.14 | 1.61 | 2877.35 | 1.0 |
| scale1 | 146.66 | 21.97 | 0.33 | 108.67 | 190.92 | 4099.34 | 1.0 |
| scale2 | 111.09 | 13.30 | 0.18 | 87.85 | 138.03 | 5590.24 | 1.0 |
pm.traceplot(trace_2)
shapes = pm.quantiles(trace_2.get_values('shape'), qlist=qlist) scales1 = pm.quantiles(trace_2.get_values('scale1'), qlist=qlist) scales2 = pm.quantiles(trace_2.get_values('scale2'), qlist=qlist)
fig = plt.figure(figsize=(9,4)) ax = fig.add_subplot(111) plot_survival_function(ax, x, shapes, scales1) plot_survival_function(ax, x, shapes, scales2, color='red')
Cox比例ハザードモデルとの比較
cph = lifelines.CoxPHFitter() cph.fit(df[['week','arrest','fin']], duration_col='week', event_col='arrest', show_progress=False) cph.print_summary() # access the results using cph.summary
n=432, number of events=114
coef exp(coef) se(coef) z p lower 0.95 upper 0.95
fin -0.3691 0.6914 0.1897 -1.9453 0.0517 -0.7409 0.0028 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Concordance = 0.546
Likelihood ratio test = 3.837 on 1 df, p=0.05013
fig = plt.figure(figsize=(9,4)) ax = fig.add_subplot(111) plot_survival_function(ax, x, shapes, scales1) plot_survival_function(ax, x, shapes, scales2, color='red') cph.plot_covariate_groups('fin',[0,1], ax=ax)
こちらもやはり、ベイズ信頼区間の範囲が元の記事よりも広い。 元の記事（stan）では事前分布が明示されていないため、そこが違いを生んでいるのかも知れない。 とはいえ、pymcでの打ち切りデータを含む場合のweibull分布への当てはめ方法が分かった。
以上
アヒル本10章棋士の強さのモデル化をpymcでやってみる。
この章の例題をpymcで実装してみた。 pymcによる実装の課題は以下の通り。
- 順序付けを行うorderedが、pymcにはない（と思われる）
- 単純にMCMCを実行するとエラーで実行できない
1点目は、以前のドラゴンボールの戦闘力推定の実装を使った。
# データの読み込み players = pd.read_csv('shogi/data-shogi-player-name.txt', sep='\t') battles = pd.read_csv('shogi/data-shogi-player.txt', sep=',')
モデル化は、ドラゴンボールの戦闘力推定モデルとほぼ一緒。
with pm.Model() as model: # prior s_mu = pm.Normal('s_mu', mu=0, sd=100) mu = pm.Normal('mu', mu=0, sd=s_mu, shape=len(players)) sigma = pm.Gamma('sigma', 10, 10, shape=len(players)) loser_perf = pm.Normal('loser_perf', mu=mu[battles.Loser-1], sd=sigma[battles.Loser-1], shape=len(battles)) winner_perf = pm.Normal('winner_perf', mu=mu[battles.Winner-1], sd=sigma[battles.Winner-1], shape=len(battles)) diff = winner_perf-loser_perf p = pm.math.sigmoid(diff) z = pm.Bernoulli('z',p=p, observed=np.ones(len(battles)))
MCMCの実行では、ADVI()による初期化を行うことで、NUTSによる実行が可能となった。 （ADVIによる初期化をしないと、エラーで落ちる）
with model: trace = pm.sample(2000, init='advi', n_init=25000)
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average Loss = 4,742.8: 100%|██████████| 25000/25000 [05:40<00:00, 73.47it/s]
Finished [100%]: Average Loss = 4,743.1
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [winner_perf, loser_perf, sigma, mu, s_mu]
Sampling 2 chains: 100%|██████████| 5000/5000 [44:36<00:00, 2.09draws/s]
The gelman-rubin statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.
pm.traceplot(trace, varnames=['mu','sigma'])
推定結果を表示する際、indexを棋士名に変換する。
def rename_indexes(trace, varname, players): ret = pm.summary(trace, varnames=[varname]) tab = {} for p in range(len(players)): tab.setdefault('%s__%d'%(varname,p), players[players.nid==p+1].kname.values[0]) ret = ret.rename(index=tab) return ret
各棋士の強さの平均を表示する。
ret_mu = rename_indexes(trace, 'mu', players) ret_mu.sort_values('mean', ascending=False).head()
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| 羽生善治 | 2.959619 | 0.268438 | 0.015966 | 2.412412 | 3.475476 | 31.911283 | 1.073706 |
| 渡辺明 | 2.538573 | 0.294463 | 0.017209 | 1.984500 | 3.114538 | 80.892714 | 1.030923 |
| 行方尚史 | 2.116030 | 0.313367 | 0.017155 | 1.464525 | 2.698630 | 223.610725 | 1.019783 |
| 豊島将之 | 2.067839 | 0.286370 | 0.015999 | 1.542303 | 2.653113 | 117.756783 | 1.045872 |
| 郷田真隆 | 1.995864 | 0.285740 | 0.016104 | 1.472823 | 2.532897 | 195.203563 | 1.026106 |
各棋士の勝負ムラを表示する。
ret_sigma = rename_indexes(trace, 'sigma', players) ret_sigma.sort_values('mean', ascending=False).head()
勝負ムラが大きい棋士。
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| 藤森哲也 | 1.126341 | 0.387405 | 0.030618 | 0.353202 | 1.831846 | 65.198876 | 0.999863 |
| 森内俊之 | 1.125223 | 0.344858 | 0.026571 | 0.511552 | 1.838470 | 89.927761 | 1.000393 |
| 深浦康市 | 1.116453 | 0.337526 | 0.025087 | 0.499118 | 1.743568 | 112.326690 | 1.003547 |
| 神崎健二 | 1.113792 | 0.336298 | 0.024777 | 0.483929 | 1.786775 | 59.970475 | 1.000581 |
| 阿久津主税 | 1.113714 | 0.327681 | 0.022573 | 0.472958 | 1.722987 | 119.652929 | 1.023349 |
勝負ムラが小さい棋士。
ret_sigma.sort_values('mean', ascending=False).tail()
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| 菅井竜也 | 0.906430 | 0.282381 | 0.022321 | 0.375220 | 1.422869 | 56.414318 | 1.013031 |
| 先崎学 | 0.904325 | 0.279065 | 0.019915 | 0.380872 | 1.418062 | 103.896750 | 0.999974 |
| 飯島栄治 | 0.899098 | 0.250615 | 0.018413 | 0.408726 | 1.358156 | 72.390061 | 1.004667 |
| 竹内雄悟 | 0.892474 | 0.282656 | 0.019096 | 0.392110 | 1.476494 | 128.417525 | 1.001682 |
| 北島忠雄 | 0.862707 | 0.272092 | 0.018809 | 0.332110 | 1.346679 | 132.050279 | 1.002114 |
書籍の結果とは微妙に異なるが、近い順位が得られた。
以上
アヒル本11章ゼロ過剰ポアソン分布をpymcでやってみる。
名前のとおり、ゼロが多いポアソン分布ということで、ベルヌーイ分布とポアソン分布を組み合わせたZIP分布を使った例題。
![StanとRでベイズ統計モデリング （Wonderful R） [ 松浦 健太郎 ]](https://thumbnail.image.rakuten.co.jp/@0_mall/book/cabinet/2421/9784320112421.jpg?_ex=128x128 "StanとRでベイズ統計モデリング （Wonderful R） [ 松浦 健太郎 ]")
stanではなく、pymcで実装する場合、以下が注意点として挙げられる。
- 切片のデータを追加しておく。
- pymcのZIP分布におけるポアソン分布のパラメータλは、リンク関数(exp)を掛けて値を入れる。 （stanの実装では、poisson(exp(x)）と等価な関数でZIP関数が作られているため、リンク関数を掛けていない。 これに気づかず、中々うまく推定できなかった。）
データを読み込み、可視化する。ゼロが多いカウントデータであることが分かる。
df = pd.read_csv('./data/data-ZIP.txt')
df.Y.plot.hist()
データの整形。切片も加える。
df.Age = df.Age/10 df['Intercept'] = 1 Y = df.Y.values X = df[['Intercept', 'Sex', 'Sake', 'Age']].values n_shape = X.shape[1]
モデルの実装。
with pm.Model() as model: b1 = pm.Normal('b1', mu=0, sd=5, shape=n_shape) b2 = pm.Normal('b2', mu=0, sd=5, shape=n_shape) q_x = pm.math.dot(X,b1) q = pm.invlogit(q_x) lam = pm.math.exp(pm.math.dot(X,b2)) # リンク関数(exp)をかけておく y_pred = pm.ZeroInflatedPoisson('y_pred', q, lam, observed=Y)
with model: trace = pm.sample(2000, start=pm.find_MAP(), step=pm.NUTS())
logp = -416.06, ||grad|| = 0.11749: 100%|██████████| 65/65 [00:00<00:00, 754.49it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [b2, b1]
Sampling 2 chains: 100%|██████████| 5000/5000 [00:27<00:00, 181.93draws/s]
pm.summary(trace)
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| b1__0 | 0.914821 | 0.674350 | 0.012534 | -0.332461 | 2.305058 | 2608.604989 | 0.999940 |
| b1__1 | 1.579079 | 0.415617 | 0.006595 | 0.779755 | 2.401231 | 4781.578735 | 0.999810 |
| b1__2 | 3.258817 | 0.788188 | 0.014349 | 1.781250 | 4.799758 | 2837.311993 | 1.000395 |
| b1__3 | -0.356855 | 0.179263 | 0.003336 | -0.716908 | -0.018116 | 2592.259435 | 1.000172 |
| b2__0 | 1.450101 | 0.134314 | 0.002736 | 1.175606 | 1.699740 | 2785.911268 | 1.000095 |
| b2__1 | -0.745450 | 0.079422 | 0.001370 | -0.893933 | -0.588722 | 3471.970062 | 0.999781 |
| b2__2 | -0.161931 | 0.073769 | 0.001336 | -0.306756 | -0.023078 | 3011.972673 | 0.999774 |
| b2__3 | 0.198165 | 0.033652 | 0.000672 | 0.131085 | 0.262932 | 2823.546198 | 1.000473 |
pm.traceplot(trace)
推測したパラメータでデータを生成し、分布を確認。
pred = pm.sample_ppc(trace, samples=100, model=model)
100%|██████████| 100/100 [00:00<00:00, 682.51it/s]
plt.hist(pred['y_pred'].reshape(-1,))
実際のデータに近い分布が得られた。
以上
アヒル本１０章多項ロジスティック回帰をpymcでやってみる。
8.4章のロジスティック回帰は、ある一つの項目（出席/欠席）についてのみ考えるが、多項ロジスティック回帰は複数の項目について考える。
カテゴリカル分布を使ってモデル化するのだが、各パラメータの１番目の要素を0に固定化するというのが味噌で、これをpymcで実装してみた。
# データの読み込み df = pd.read_csv('data/data-category.txt') # データの整形 X = df[['Age','Sex','Income']].values Y = df.Y - 1 K = len(df.Y.unique()) # カテゴリ数
モデルの実装。
with pm.Model() as model: # パラメータ数は、後で0を加えるため、カテゴリ数-1になる。 alpha = pm.Normal('alpha', mu=0, sd=100, shape=(K-1,)) beta = pm.Normal('beta', mu=0, sd=100, shape=(3,K-1)) # 各パラメータの1つを0に固定化する。 alpha_f = pm.Deterministic('alpha_f', tt.concatenate([tt.zeros(1), alpha])) beta_f = pm.Deterministic('beta_f', tt.concatenate([tt.zeros((3,1)) , beta], axis=1)) # mu = alpha_f + pm.math.dot(X, beta_f) theta = tt.nnet.softmax(mu) y = pm.Categorical('y', p=theta, observed=Y)
with model: start = pm.find_MAP(method='Powell') step = pm.NUTS() trace_sf = pm.sample(3000, step, start)
logp = -inf, ||grad|| = 0: 100%|██████████| 510/510 [00:00<00:00, 637.82it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha]
Sampling 2 chains: 100%|██████████| 7000/7000 [03:05<00:00, 37.79draws/s]
MCMC実行結果。
pm.traceplot(trace_sf)
書籍には推定結果は載っていなかったため、データと推測結果を比較してみた。
# 推定結果の取り出し alpha_f = np.mean(trace_sf['alpha_f'], axis=0) beta_f = np.mean(trace_sf['beta_f'], axis=0)
各データで一番可能性が高いカテゴリを予測する。
preds = [] for uid in df.index: mu = alpha_f + np.dot(X[uid], beta_f) theta = tt.nnet.softmax(mu).eval()[0] argmax = np.argmax(theta) preds.append(argmax)
実際の結果と予測結果を比較する。0が最も多いはず。
(Y - preds).value_counts()
0 101
1 40
2 32
-2 26
4 23
-1 17
3 16
5 14
-4 13
-3 11
-5 7
Name: Y, dtype: int64
実際-予測=0が一番多かった。ということで、うまく推定できたと判断しよう。
以上
アヒル本8.4章ロジスティック回帰階層モデルをpymcでやってみる。
この例題は、学生の出席/欠席を、アルバイト好きか、学問への興味、天気などのデータから予測するもので、 学生差や科目差を考慮した階層モデルになっている。
これを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))
大体同じような結果が得られたと思う。
以上
アヒル本8.3章非線形モデルの階層モデルをpymcでやってみる。
StanとRでベイズ統計モデリングの8.3「非線形モデルの階層モデル」をpymcで実装してみる、というお話。 この本は実践的な例題が多く、非常に参考になる。ただし、タイトルの通りだが、RとStanの実装。
- この節の例題は、経過時間と薬の血中濃度を非線形の関係でモデル化している。
- 各患者のモデルパラメータを階層化してモデル化しており、パラメータの対数値が正規分布に従う、と仮定されている。
- pymcでモデルを実装する際、時系列やパラメータの対数変換はどうするのか？という素朴な疑問からpymcで実装してみた。
実装してみた結果は以下。
- 患者と時間という2軸あるので、少し分かりにくい印象。(Stanの方がコード化しやすい気がする)
- NUTSによるMCMCはエラーで実行できなかった。代わりにADVIで推定した。得られた推定結果は書籍と大体同じであった。
df = pd.read_csv('data/data-conc-2.txt')
df.head()
| PersonID | Time1 | Time2 | Time4 | Time8 | Time12 | Time24 | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 2.4 | 5.0 | 7.5 | 11.9 | 12.5 | 12.7 |
| 1 | 2 | 1.4 | 3.9 | 4.4 | 7.7 | 6.4 | 8.3 |
| 2 | 3 | 5.2 | 9.4 | 19.4 | 20.2 | 22.7 | 24.9 |
| 3 | 4 | 6.7 | 12.6 | 19.1 | 23.4 | 25.8 | 26.1 |
| 4 | 5 | 0.3 | 4.7 | 7.0 | 10.2 | 12.9 | 14.8 |
N = len(df.PersonID.unique()) Times = np.array([1,2,4,8,12,24]) Times = np.array([Times]).T # 列ベクトルに変換する。
with pm.Model() as model: # 階層化部分 a0 = pm.Normal('a0', mu=0, sd=10**2) s_a = pm.HalfCauchy('s_a', beta=10**2) b0 = pm.Normal('b0', mu=0, sd=10**2) s_b = pm.HalfCauchy('s_b', beta=10**2) # 各患者のモデルパラメータ log_a = pm.Normal('log_a', mu=a0, sd=s_a, shape=N) a = pm.Deterministic('a', pm.math.exp(log_a)) log_b = pm.Normal('log_b', mu=b0, sd=s_b, shape=N) b = pm.Deterministic('b', pm.math.exp(log_b)) # 誤差 eps = pm.Normal('eps', mu=0, sd=10**2) # 各時刻における各患者の平均（matrixになるので、少しややこしい） mu = pm.Deterministic('mu', a*(1-pm.math.exp(-b*Times))) y = pm.Normal('y', mu=mu, sd=eps, observed=df.iloc[:,1:].values.T)
パラメータの推定。MCMC（NUTS）はエラーで実行できなかった。
with model: inference = pm.ADVI() approx = pm.fit(100000, method=inference)
Average Loss = 265.72: 100%|██████████| 100000/100000 [02:42<00:00, 614.61it/s]
Finished [100%]: Average Loss = 265.69
trace = approx.sample(draws=3000)
推定結果。書籍の結果と大体同じ値が得られている。
pm.summary(trace, varnames=['a0','b0','s_a','s_b','eps'])
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | |
|---|---|---|---|---|---|
| a0 | 2.859904 | 0.110053 | 0.001996 | 2.636606 | 3.068902 |
| b0 | -1.165525 | 0.089034 | 0.001613 | -1.344461 | -0.993055 |
| s_a | 0.420661 | 0.083319 | 0.001390 | 0.275964 | 0.590157 |
| s_b | 0.344515 | 0.071129 | 0.001492 | 0.210082 | 0.483927 |
| eps | 1.808791 | 0.139476 | 0.002373 | 1.540795 | 2.083439 |
グラフ化して確認。 各患者の平均値の推定値(事後分布の中央値)と、95%予測区間をグラフ化する。
times = Times.T[0] # 事後分布の中央値を使う。 mu = pm.quantiles(trace)['mu'][50] eps = pm.quantiles(trace)['eps'][50] _, axes = plt.subplots(4, 4) for (row, col), ax in np.ndenumerate(axes): id = row * 4+ col # 観測値 ax.scatter(times, df.loc[id,'Time1':].values, color='k') # 平均の推定値（事後分布の中央値） ax.plot(times, mu[:,id], color='k') # 平均の推定値（事後分布の95%最高密度区間） y1,y2 = map(list, zip(*pm.hpd(trace, alpha=0.05)['mu'][:,id])) #ax.fill_between(times, y1=y1, y2=y2, alpha=0.2, color='b') # 平均値の推定値（事後分布の中央値）を使った場合の95%予測区間 y1,y2 = map(list, zip(*norm.ppf(q=[0.025, 0.975], loc=np.array([mu[:,id]]).T, scale=eps))) ax.fill_between(times, y1=y1, y2=y2, alpha=0.2, color='r') # グラフの調整 if row < 3: plt.setp(ax.get_xticklabels(), visible=False) else: plt.setp(ax, xlabel='Time (hour)') if col > 0: plt.setp(ax.get_yticklabels(), visible=False) else: plt.setp(ax, ylabel='Y') plt.setp(ax, title='id:%d'%(id+1), xticks=times, xlim=(0, 24), yticks=np.arange(0, 40, 10), ylim=(-3, 37)) plt.tight_layout()
グラフも大体同じ結果。
以上
