StanとPythonでベイズ統計モデリング その2 Chapter5
アヒル本(StanとRでベイズ統計モデリング)のChapter5にPythonで取り組んでいきます。 練習問題を解いて、本文中に書かれてるグラフをPythonで描いてみます。
なおChapter1~3は導入だったのと、Chapter4は練習問題の内容が「はじめての統計データ分析」と被っていたのでパスします。
Chapter5 基礎的な回帰とモデルのチェック
- 重回帰
- 複数の説明変数を用いた回帰のこと
- 重回帰も結局は正規分布を仮定している
- 目的
- 説明変数からの応答変数の予想、及び説明変数の寄与率
- 分布
- 複数の説明変数ならScatterplot matrixを利用すると良い
- MCMCの設定について
- 結果の解釈 : モデルの改善に活かす!
- 複数の説明変数を用いた回帰のこと
- (二項)ロジスティック回帰
- ロジスティック関数をつかった方法
- 応答変数が0から1の間にある場合、説明変数に適用することで値の範囲を変換することが出来る
- 2項分布モデル
- ベルヌーイ分布モデル
- 応答変数がbool値であるような時にはベルヌーイ分布を利用すれば良い
- 生成量
- オッズ
- 二項ロジスティック回帰モデルにおいて、であり、であるならばと変換することが出来る
- はオッズと呼ばれる
- その事象が起こるサンプリング確率の比率を示す
- オッズ
- 可視化
- ロジスティック関数をつかった方法
- ポアソン回帰
- カウントデータに対して予測を行う時に用いる
- でカウントデータYが決定されるものとして、は説明変数の線形結合を0以上に変換したものを使う(がよく使われる)
- パフォーマンスの理由より、が十分に大きいなら、正規分布で近似して良い
- カウントデータに対して予測を行う時に用いる
練習問題
Scatterplot matrixについて
Pythonの実装はpandasマスターのid:sinhrksさんの記事を見るのが良いと思います。
個人的な使った感想は次の通りです。
- 離散値の可視化が上手く行かないことがある(6離散値指定したのに、4つしか表示されないことがある)?
- たぶんこれはseabornのPairPlotのバグです……
- 離散値/連続値判定の初期値が大きすぎると上手くいかないので、離散値の判定閾値を4, 5辺りにいじって使いました。
- 離散値の場合は対角のカウントはcountplotを使ったほうが良い?
- 他の下半分、上半分のグラフと軸を合わせやすくするためです。
- どうせseabornはPairGridで使うことになりますし…
そういうわけで改造させて頂いたバージョンがこちらです。残念ながら1つ目のバグは取れていません。そもそも僕の環境に特有なのかもしれませんし、諦めました。 上三角行列の部分はid:statmodelingさんの記事のまんまです。
#!/usr/bin/env python # vim:fileencoding=utf-8 # Created: 2017-06-06 # Original scripts: http://sinhrks.hatenablog.com/entry/2016/11/01/075527 import matplotlib.pyplot as plt import numpy as np import pandas as pd import seaborn as sns from matplotlib.patches import Ellipse class Dispatcher(object): def __init__(self, fontsize=20, alpha=0.6, cmap='RdBu', threshold=5): self.fontsize = fontsize self.alpha = alpha self.cmap = plt.get_cmap(cmap) # 離散値 / 連続値とみなす閾値 self.threshold = threshold def comb(self, x_series, y_series, label=None, color=None): """ 下三角部分のプロット """ x_nunique = x_series.nunique() y_nunique = y_series.nunique() if x_nunique < self.threshold and y_nunique < self.threshold: # 離散値 x 離散値のプロット return self._dd_plot(x_series, y_series, label=label, color=color) elif x_nunique < self.threshold or y_nunique < self.threshold: # 離散値 x 連続値のプロット return self._dc_plot(x_series, y_series, label=label, color=color) else: # 連続値 x 連続値のプロット return plt.scatter(x_series, y_series, label=label, color=color) def _dd_plot(self, x_series, y_series, label=None, color=None): """ 離散値 x 離散値のプロット """ # x, y 各組み合わせの個数を集計 total = y_series.groupby([x_series, y_series]).count() # x, y軸をプロットする位置を取得 xloc = total.index.labels[0] yloc = total.index.labels[1] values = total.values ax = plt.gca() for xp, yp, vp in zip(xloc, yloc, values): ax.annotate(vp, (xp, yp), fontsize=self.fontsize, ha='center', va='center') # 組み合わせの個数を散布図としてプロット size = values / (values.max() * 1.1) * 100 * 20 ax.scatter(xloc, yloc, s=size, label=label, color=color) ax.set_ylim(yloc[0] - 0.5, yloc[-1] + 0.5) def _dc_plot(self, x_series, y_series, label=None, color=None): """ 離散値 x 連続値のプロット """ if y_series.nunique() < x_series.nunique(): # y軸が離散値の場合は、x, yを入替 # 水平方向に箱ひげ図をプロット orient = "h" else: orient = "v" ax = sns.boxplot(x_series, y_series, orient=orient, color=color) ax = sns.swarmplot(x_series, y_series, orient=orient, color=color) def diag(self, x_series, label=None, color=None): """ 対角部分のプロット """ x_nunique = x_series.nunique() if x_nunique < self.threshold: ax = sns.countplot(x_series, color=color) else: ax = sns.distplot(x_series, kde=False, color=color) ax.yaxis.set_visible(False) def ellipse(self, x_series, y_series, label=None, color=None): """ 上三角部分のプロット """ # 相関係数を楕円としてプロット r = x_series.corr(y_series) c = self.cmap(0.5 * (r + 1)) ax = plt.gca() ax.axis('off') artist = Ellipse(xy=[.5, .5], width=np.sqrt(1+r), height=np.sqrt(1-r), angle=45, facecolor=c, edgecolor='none', transform=ax.transAxes) ax.add_artist(artist) ax.text(.5, .5, '{:.0f}'.format(r*100), fontsize=28, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes) d = Dispatcher() # データを作る x = np.random.normal(100, 5, 30) y = np.random.choice([1,2,3], 30, p=[0.1,0.2,0.7]) z = np.random.choice([1,2,3,4], 30, p=[0.2,0.2,0.2,0.4]) w = np.random.normal(4, 0.1, 30) df = pd.DataFrame(np.array([x,y,z,w]).T, columns=["x", "y", "z", "w"]) # 可視化する g = sns.PairGrid(df, diag_sharey=False) g.map_diag(d.diag) g.map_lower(d.comb) g.map_upper(d.ellipse) plt.show()
これとは別に、勉強として自分でも実装してみましたが、離散×離散の可視化が難しいところです。僕はヒートマップを使おうかと思ったのですが、軸を表示しようとする関係でズレます。つらい。バブルプロットを使うのがテキスト通りの適切なアイディアです。 一応載せておきます。
import matplotlib.pyplot as plt import numpy as np import pandas as pd import seaborn as sns from matplotlib.patches import Ellipse from scipy import stats def value_type(v): if v.unique().size < 5: v_type = "discrete" else: v_type = "continuous" return v_type def correlation(x, y, **kws): r, _ = stats.spearmanr(x, y) ax = plt.gca() ax.axis('off') ellcolor = plt.cm.RdBu(0.5*(r+1)) txtcolor = 'black' if np.absolute(r) < 0.5 else 'white' artist = Ellipse(xy=[.5, .5], width=np.sqrt(1+r), height=np.sqrt(1-r), angle=45, facecolor=ellcolor, edgecolor='none', transform=ax.transAxes) ax.add_artist(artist) ax.text(.5, .5, '{:.0f}'.format(r*100), color=txtcolor, fontsize=28, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes) def diagonal(x, **kws): ax = plt.gca() x_type = value_type(x) if x_type == "discrete": sns.countplot(x) else: sns.distplot(x, kde=False) def distribution(x, y, **kws): ax = plt.gca() x_type = value_type(x) y_type = value_type(y) if x_type == "discrete": if y_type == "discrete": df = pd.DataFrame([x, y]).T df.columns=[x.name, y.name] pdf = df.pivot_table(index=x.name, columns=y.name, aggfunc=lambda x: len(x)).fillna(0).T plt.pcolor(pdf) else: sns.boxplot(x, y) sns.swarmplot(x, y) else: if y_type == "discrete": sns.boxplot(x, y, orient="h") sns.swarmplot(x, y, orient="h") else: plt.scatter(x, y) def scatterplot_matrix(df): g = sns.PairGrid(df) g.map_upper(correlation) g.map_diag(diagonal) g.map_lower(distribution) return g x = np.random.normal(100, 5, 30) y = np.random.choice([1,2,3], 30, p=[0.1,0.2,0.7]) z = np.random.choice([1,2,3,4], 30, p=[0.4,0.3,0.2,0.1]) w = np.random.normal(4, 0.1, 30) df = pd.DataFrame(np.array([x,y,z,w]).T, columns=["x", "y", "z", "w"]) g = scatterplot_matrix(df) plt.show(g)
練習問題(1)
まぁほとんど弄る部分もありません。 いつも通りPyStanを回すだけです。
import pystan import pandas as pd import os import pickle from IPython.core.display import display stan_file = os.path.join("RStanBook", "chap05", "model", "model5-3.stan") model = pystan.StanModel(stan_file) df = pd.read_csv("./RStanBook/chap05/input/data-attendance-1.txt") N = len(df) A = df["A"].values Score = df["Score"].values / 200 Y = df["Y"].values stan_data = { "N": N, "A": A, "Score": Score, "Y": Y } fit = model.sampling( data=stan_data, iter=2100, chains=4, warmup=100, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df) mu = fit.extract()["mu"] diff = np.apply_along_axis(lambda x: Y - x, axis=1, arr=mu)
モデルはgithubで公開されているものをそのまま使いました。
練習問題(2)
Stanの内部の生成量を出すところで残差を計算します。 個人的には汚い部分を全部Stanの中に隠せるのでこっちのほうが好きです。
import pystan import pandas as pd import os import pickle from IPython.core.display import display stan_file = os.path.join("RStanBook", "chap05", "model", "model5-3_kai.stan") stan_file_c = os.path.join("RStanBook", "chap05", "model", "model5-3_kai.pkl") # コンパイルチェック if os.path.exists(stan_file_c) == False: # コンパイル済みのバイナリが無かったらstanファイルからコンパイル model = pystan.StanModel(stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) else: # バイナリがあったらそれを利用 with open(stan_file_c, "rb") as f: model = pickle.load(f) df = pd.read_csv("./RStanBook/chap05/input/data-attendance-1.txt") N = len(df) A = df["A"].values Score = df["Score"].values / 200 Y = df["Y"].values stan_data = { "N": N, "A": A, "Score": Score, "Y": Y } fit = model.sampling( data=stan_data, iter=2100, chains=4, warmup=100, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df)
モデルは次の通りです。
data { int N; int<lower=0, upper=1> A[N]; real<lower=0, upper=1> Score[N]; real<lower=0, upper=1> Y[N]; } parameters { real b1; real b2; real b3; real<lower=0> sigma; } transformed parameters { real mu[N]; for (n in 1:N) mu[n] = b1 + b2*A[n] + b3*Score[n]; } model { for (n in 1:N) Y[n] ~ normal(mu[n], sigma); } generated quantities { real y_pred[N]; real diff[N] ; for (n in 1:N){ y_pred[n] = normal_rng(mu[n], sigma) ; diff[n] = Y[n] - mu[n] ; } }
つづいてpythonでグラフを描いて、目視で比較をします。
まず実測値とモデルによる予測値を比較してみます。 テキストでは予測分布を使っていますが、試しにMCMCで推定されたパラメータ分布の中央値を使って値を予測してみます。この目的は、後述する事後予測分布との比較です。これは「はじめての統計データ分析」で言うところの、条件付き予測分布を使った手法です。こちらの手法だとStanファイル中で予測をしなくていい代わりに、エラーバーをつけるのがめんどくさいです。何故ならYを導出するためのにもにも分布があるので、その両方を考慮する必要があります。今回はパスさせて下さい。
import numpy as np import seaborn as sns import matplotlib.pyplot as plt # Aの値でindexを分類する A0_samples = df[df["A"] == 0].index A1_samples = df[df["A"] == 1].index # summaryのtableからmuの中央値とsigmaを取得する summary_df_A0 = summary_df.loc[["mu[{0}]".format(i) for i in A0_samples], :] summary_df_A1 = summary_df.loc[["mu[{0}]".format(i) for i in A1_samples], :] map_A0 = summary_df_A0["50%"] map_A1 = summary_df_A1["50%"] sigma = summary_df.loc["sigma", "50%"] # パラメータを使って条件付き予測分布を推定する Y_A0 = Y[A0_samples] Y_A1 = Y[A1_samples] Y_A0A1 = np.hstack([Y_A0, Y_A1]) Y_ast_A0 = np.random.normal(map_A0, sigma) Y_ast_A1 = np.random.normal(map_A1, sigma) Y_ast_A0A1 = np.hstack([Y_ast_A0, Y_ast_A1]) # グラフのために相関係数を導出する r = np.corrcoef(Y_A0A1, Y_ast_A0A1)[0][1] # 可視化 fig, ax = plt.subplots() ax.plot(Y_ast_A0, Y_A0, "o", label="A=0") ax.plot(Y_ast_A1, Y_A1, "^", label="A=1") ax.plot([0, 0.5], [0, 0.5], "k-") ax.legend() # 体裁を整える ax.set_aspect('equal') ax.set_xlabel("Predicted value(by EAP)") ax.set_ylabel("Observed value") ax.set_title("r={0}".format(r)) plt.show()
Aについてバラけているわけではないですが、あまり精度が良く無さそうです。モデルの精度が良くないのでしょうか(後述しますが、これは違います)?
続いて、MCMCサンプリングの中でyの予測値を推定します。これは「はじめての統計データ分析」で言うところの、事後予測分布を使った手法です。こちらのほうが分布全体を考慮して行っているので、正確な予測値を計算することが出来ます。エラーバーについてもMCMCの分布を使えばいいので簡単です。今回は教科書の80%区間ではなく、デフォルトで計算される75%区間を利用しています。
# Aの値でindexを分類する A0_samples = df[df["A"] == 0].index A1_samples = df[df["A"] == 1].index Y_A0 = Y[A0_samples] Y_A1 = Y[A1_samples] Y_A0A1 = np.hstack([Y_A0, Y_A1]) # プロット点の計算 summary_df_A0 = summary_df.loc[["y_pred[{0}]".format(i) for i in A0_samples], :] summary_df_A1 = summary_df.loc[["y_pred[{0}]".format(i) for i in A1_samples], :] Y_pred_A0 = summary_df_A0["50%"] Y_pred_A1 = summary_df_A1["50%"] Y_pred = np.hstack([Y_pred_A0, Y_pred_A1]) # 誤差の計算 Y_pred_A0_err_down = (Y_pred_A0 - Y_pred_summary_A0["25%"]).values Y_pred_A1_err_down = (Y_pred_A1 - Y_pred_summary_A1["25%"]).values Y_pred_A0_err_up = (Y_pred_summary_A0["75%"] - Y_pred_A0).values Y_pred_A1_err_up = (Y_pred_summary_A1["75%"] - Y_pred_A1).values Y_pred_A0_err = np.vstack([Y_pred_A0_err_down, Y_pred_A0_err_up]) Y_pred_A1_err = np.vstack([Y_pred_A1_err_down, Y_pred_A1_err_up]) # 相関係数の計算 r = np.corrcoef(Y_A0A1, Y_pred)[0][1] # 可視化 fig, ax = plt.subplots() ax.errorbar(Y_pred_A0, Y_A0, yerr=Y_pred_A0_err, fmt="o", label="A=0") ax.errorbar(Y_pred_A1, Y_A1, yerr=Y_pred_A1_err, fmt="^", label="A=1") # y=xの線を書く ax.plot([0, 0.5], [0, 0.5], "k-") # 凡例をつける ax.legend() # 体裁を整える ax.set_aspect('equal') ax.set_xlabel("Predicted value") ax.set_ylabel("Observed value") ax.set_title("r={0}".format(r)) plt.show()
結果がわりと異なります。相関係数に如実に現れました。条件付き予測分布ではダメかと思ったモデルですが、こちらでは問題がなさげなことが推定されます。以上のことから、条件付き予測分布はオンラインで利用できるなどの利点がありますが、その分結果の浮動性があることを頭に入れたほうが良いと思います。とはいえ、どちらの方法でもAに関する偏ったバラつきが無いことはわかりますので、何を見たいかの使い分けが結局のところ肝心です。
最後に実測値とパラメータの誤差を調べます。
summary_df_diff = summary_df.loc[["diff[{0}]".format(i) for i in range(0, 50)], :] fig, ax = plt.subplots() sns.distplot(summary_df_diff["50%"], bins=np.arange(-0.15, 0.15, 0.03), ax=ax) ax.set_xlabel("Epsilon EAP") ax.set_ylabel("Count") plt.show()
教科書の結果から分かっていましたが、きちんとあてはめた正規分布に従っています。
練習問題(3)
離散値と離散値の集計です。pandasのpivot_tableメソッドを使えばいいのですが、芸がないのでヒートマップで可視化しておきます。AとYの各組み合わせで足していけばいいので、ラムダ式を読み出してlenでカウントしました。
import pandas as pd import seaborn as sns import matplotlib.pyplot as plt df = pd.read_csv("./RStanBook/chap05/input/data-attendance-3.txt") pdf = df[["A", "Y"]].pivot_table(index="A", columns="Y", aggfunc=lambda x:len(x)) g = sns.heatmap(pdf, annot=True) plt.show(g)
練習問題(4)
晴れがもたらす影響を基準値として、曇と雨が持つ影響についてもMCMCから推定してみます。これだけだと芸が無いので、WAICによるモデルのフィッティングについても評価してみます。
import pystan import pandas as pd import os import pickle import numpy as np from IPython.core.display import display from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger logger = get_logger() stan_file = os.path.join("RStanBook", "chap05", "model", "model5-5_kai.stan") model = pystan.StanModel(stan_file) trans = {"A": 1, "B": 2, "C": 3} df = pd.read_csv("./RStanBook/chap05/input/data-attendance-3.txt") I = len(df) A = df["A"].values Score = df["Score"].values / 200 W = df["Weather"].values W = np.array([trans[w] for w in W]) Y = df["Y"].values stan_data = { "I": I, "A": A, "Score": Score, "W": W, "Y": Y } fit = model.sampling( data=stan_data, iter=1200, chains=4, warmup=200, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df) log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
モデルは以下の通りです。
data { int I; int<lower=0, upper=1> A[I]; real<lower=0, upper=1> Score[I]; int<lower=1, upper=3> W[I]; int<lower=0, upper=1> Y[I]; } parameters { real b[3]; real b_w[2] ; } transformed parameters { real q[I]; real beta_w[3] ; beta_w[1] = 0 ; beta_w[2] = b_w[1] ; beta_w[3] = b_w[2] ; for (i in 1:I) q[i] = inv_logit(b[1] + b[2]*A[i] + b[3]*Score[i] + beta_w[W[i]]); } model { for (i in 1:I) Y[i] ~ bernoulli(q[i]); } generated quantities { real log_lik[I] ; for (i in 1:I){ log_lik[i] = bernoulli_lpmf(Y[i]|q[i]) ; } }
さて経験的にパラメータを決定した時と比べて、WAICはどうなったでしょうか。
- 経験的に決めた場合 : WAIC=2769
- MCMCからサンプリングした場合: WAIC=2764
ちょっとだけ、MCMCからサンプリングした場合のほうが値が良くなります。EAPを見てモデルを書いてみると
- 経験的に決めた場合
- MCMCからサンプリングした場合
となります。Weatherの部分も係数の-0.45を考慮すると、そこまで変わりません。僕は経験的に得られた0:0.2:1が結構いいものであったと解釈しました。
つづいてpythonでグラフを描いて、目視で結果を比較します。 テキストに描いてあるように、ロジスティック回帰の観測値は0/1になるので、可視化が難しいところです。qの値が、データからどのように変化するかを調べます。テキストの例では連続値のScoreと離散値のAが可視化されています。
まずScoreの影響を可視化します。
- 推測したqの中央値をデータと対応付けるために、indexを振り直します。
- カラムを複数使ったデータの選択なので、queryメソッドで取得します。
- qの中央値の線をきれいに書くためにソートを挟んでいます
# qの値を取得 summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :] q_range = summary_df_q[["25%", "50%", "75%"]] q_range.index = range(len(q_eap)) # データと結合(順番はこの時点では保存されている) df_concat = pd.concat([df, q_range], axis=1) # 一部のデータだけ選択 df_filtered = df_concat.query("A == 0 & Weather == 'A'") df_filtered = df_filtered.sort_values("Score") # Yの点の位置をそれぞれちょっとずらす y_randomized = [y + np.random.normal(0, 0.05) for y in df_filtered["Y"]] # 可視化 fig, ax = plt.subplots() ax.plot(df_concat_filtered["Score"], df_concat_filtered["50%"], "k-") ax.fill_between(df_concat_filtered["Score"], df_concat_filtered["25%"], df_concat_filtered["75%"], facecolor="pink") ax.plot(df_concat_filtered["Score"], y_randomized, "o", label="Y") ax.legend() ax.set_xlabel("Score") ax.set_ylabel("q") plt.show()
次にAに対するプロットを書きましょう。matplotlibで頑張るのが面倒だったのでseabornで書きました(なんかごめんなさい)。適材適所です!
# qの値を取得 summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :] q_range = summary_df_q[["25%", "50%", "75%"]] q_range.index = range(len(q_eap)) # データと結合(順番はこの時点では保存されている) df_concat = pd.concat([df, q_range], axis=1) # 可視化 fig, ax = plt.subplots() sns.violinplot(data=df_concat, x="50%", y="Y", ax=ax, orient="h") sns.swarmplot(data=df_concat, x="50%", y="Y", hue="A", ax=ax, orient="h", palette=sns.color_palette("Set2", n_colors=2)) ax.set_xlabel("q") ax.set_ylabel("Y") plt.show()
最後にROCカーブを書いて、AUCを計算します。これについてはScikit-learnを使うのが簡単でしょう。
from sklearn.metrics import roc_curve, roc_auc_score summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :] q = summary_df_q["50%"] # AUCの計算 auc = roc_auc_score(Y, q) # ROCカーブの座標を計算 false_positive_rate, true_positive_rate, _ = roc_curve(Y, q) fig, ax = plt.subplots() ax.plot(false_positive_rate, true_positive_rate) # y=xを描く ax.plot([0, 1], [0, 1], "k-") # 体裁を整える ax.set_aspect('equal') ax.set_xlabel("False Positive Rate") ax.set_ylabel("True Positive Rate") ax.set_title("AUC = {0}".format(auc)) plt.show()
練習問題(5)
import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from IPython.core.display import display stan_file = os.path.join("RStanBook", "chap05", "model", "model5-6b.stan") model = pystan.StanModel(stan_file) trans = {"A": 1, "B": 2, "C": 3} df = pd.read_csv("./RStanBook/chap05/input/data-attendance-2.txt") N = len(df) A = df["A"].values Score = df["Score"].values / 200 M = df["M"].values stan_data = { "N": N, "A": A, "Score": Score, "M": M } fit = model.sampling( data=stan_data, iter=2100, chains=4, warmup=100, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df) Mpred = np.median(fit.extract()["m_pred"], axis=0) fig = plt.figure(figsize=(5,5)) ax = fig.gca() ax.plot(Mpred, M, "o") ax.plot(range(0, 100), range(0, 100)) ax.set_xlim(0, 100) ax.set_ylim(0, 100) ax.set_xlabel("Predicted") ax.set_ylabel("Measured") plt.show(ax)
全然y=xに乗っていません…(合っているのかめちゃめちゃ不安です)。
この結果からは、予測が出来ていないことが分かります。説明変数が足りてないということですね。よくあることですが、計測が終わった後にこうなった場合は、新しくメタデータを足して計測し直しになるわけです。つらい。
練習問題(6)
練習問題(5)と(6)はちょっと気をつける必要があります。アヒル本だけ読んだ時には「処理」の違いというものが
- 両方に別々の処理を施した
- 片方には処理を施して、片方には何もしなかった
という2つのパターンのどちらなのか分かりません。前者の場合にはベースラインとなるものに別々の効果が生じることを予測できますが、後者の場合には一方だけはベースラインのままで、もう一方に何か効果が現れることになります。
データの元となった緑本を読んでみると、後者であることが分かります。前者を踏まえてモデルを作ってみると、MCMCが収束せずに頭を悩まるハメになります(僕は緑本を読んだはずなのに忘れていて喘ぎました)。
前者の場合でも相対的な効果を踏まえると同じじゃないの?と思っていたのですが収束しないですし、そうではないんですね。ちょっとした発見でした。これはChapter 10の一番初めに出てくる識別可能性の問題に引っかかるからですね。
import pystan import pandas as pd import os import pickle from IPython.core.display import display stan_file = os.path.join("RStanBook", "chap05", "model", "practice6.stan") model = pystan.StanModel(stan_file) trans = {"C": 0, "T": 1} df = pd.read_csv("./RStanBook/chap05/input/data3a.csv") N = len(df) X = df["x"].values F = df["f"].values F = np.array([trans[x] for x in F]) Y = df["y"].values stan_data = { "N": N, "X": X, "F": F, "Y": Y } fit = model.sampling( data=stan_data, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df)
モデルは次です。
data { int N; real<lower=0> X[N]; int<lower=0, upper=1> F[N]; int<lower=0> Y[N]; } parameters { real b[3] ; } transformed parameters { real<lower=0> lambda[N] ; for (n in 1:N){ lambda[n] = exp(b[1] + b[2]*X[n] + b[3] * F[n]) ; } } model { for (n in 1:N){ Y[n] ~ poisson(lambda[n]); } }
練習問題(7)
(6)と同様に処理の効果にだけ気をつけて下さい。
import pystan import pandas as pd import os import pickle from IPython.core.display import display stan_file = os.path.join("RStanBook", "chap05", "model", "practice7.stan") model = pystan.StanModel(stan_file) trans = {"C": 0, "T": 1} df = pd.read_csv("./RStanBook/chap05/input/data4a.csv") I = len(df) N = df["N"].values X = df["x"].values F = df["f"].values F = np.array([trans[x] for x in F]) Y = df["y"].values stan_data = { "I": I, "N": N, "X": X, "F": F, "Y": Y } fit = model.sampling( data=stan_data, iter=1500, chains=4, warmup=500, seed=1234, ) summary = fit.summary() summary_df = pd.DataFrame( summary["summary"], index =summary["summary_rownames"], columns=summary["summary_colnames"] ) display(summary_df)
モデルは以下の通りです。
data { int I; int<lower=0> N[I]; real<lower=0> X[I]; int<lower=0, upper=1> F[I]; int<lower=0> Y[I]; } parameters { real b[3] ; } transformed parameters { real<lower=0> q[I] ; for (i in 1:I){ q[i] = inv_logit(b[1] + b[2]*X[i] + b[3]*F[i]) ; } } model { for (i in 1:I){ Y[i] ~ binomial(N[i], q[i]) ; } }
練習問題(4), (6)と(7)は初めに取り組んだ時に、モデルが全然収束しなくて足踏みしました。そういう時にはMCMCサンプリングが無理をして変な部分を探索しようとするのでめちゃめちゃ時間がかかります。なので待機時間で「失敗したな」と悟れたりするのはちょっと面白いところです。
StanとPythonでベイズ統計モデリング その1
StanとRでベイズ統計モデリング(通称アヒル本)をだいたい読みました。
StanとRでベイズ統計モデリング (Wonderful R)
- 作者: 松浦健太郎,石田基広
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
- この商品を含むブログ (8件) を見る
本の紹介
既に様々な書評もありますし、方々から賛辞の声を挙げられている本です。僕としても非常に分かりやすく、使える本だと感じました。著者の松浦さんがウリを書いてくださっているので、まずそれを読むのが良いと思います。様々な方の書評も纏められています。
読んでみての感想を、良い点と改訂版に期待する点(笑)で書いてみたいと思います。
良い点
使えそうな分布が結構紹介されている
これは僕の専門分野が生命情報解析、著者の松浦さんの専門が医療統計で近いということもあるのですが、使えそうでありつつ統計入門レベルでは出てこない分布が実例を交えて紹介されています。例えば
- コーシー分布
- ゼロ過剰ポアソン分布
- ディリクレ分布
などです。ここらへんは数式だけではよく分かりませんし、とっつきづらいなところです。Stanと交えながら活用例を紹介してくださっているのは、イメージがつかみやすかったです。
階層モデル以外のベイジアンモデリングがふんだんに登場する
ベイジアンモデリングに挑戦してみよう、となると久保先生の通称緑本が第一の選択肢にあがります。僕も1年前に読みました。緑本はベイズ化する前の線形モデルから始まり、ベイズまで展開していくことで個人差が生じる過程をどう解析するのか非常に分かりやすく理解することが出来ます。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (29件) を見る
ただ、場所差、個体差などのカテゴリーの違いで本が終わってしまうので「ベイズ化すればパラメータ推定が楽になるんだな」ぐらいの感覚におちついてしまいました。これは緑本が「ベイジアンモデル」の本ではなく、「モデリング」の本ということです。
アヒル本はベイズモデルに特化しており、階層モデル以外にも
- 状態空間モデル
- トピックモデル
- 打ち切りや外れ値への対応
が紹介されています。そのためベイジアンモデルの「あいまいさ自体をパラメータ化してしまう」という特徴が他にどのように活かせるのか理解しやすく、応用範囲の広さにも気が付くことが出来ました。
ベイズモデルの難しいところは、「パラメータを生み出す分布がある」という感覚の理解かと思っています。トピックモデルなんかは、これがないと数式がごっついので難解ですし。アヒル本はその感覚の養成に役立ちます。
有用なMCMCを使ったベイジアンモデルの参考書
Stanを使ったMCMCでのモデルフィッティングにおいて、想定される問題点の解決策が述べられています。
ここらへんはStanのマニュアルを読めば書いてある事象ではあります。ただ僕は読んでみても行間が抜けすぎていて、わからない部分が結構出てきました。 なんたって、600ページありますから…。
アヒル本はマニュアルの時には行間として省略されていたところを初学者向けに補足してくれているので、緑本を読んだ後ぐらいだったら殆ど詰まること無く読めました。現状だと確率的プログラミング言語としてEdwardがアツい印象がありますが、アヒル本はベイジアンモデルの内容を包括的に書いてくれているのでStanが使えなくなったから無駄になる知識ではなさそうです。
文法が新しい
Stanは2.10.0辺りで新しい文法が提案されたようです。そのためマニュアルを読んでも、本中で紹介されているtarget記法が出てこなかったりします。これは2.9.0を参考に和訳されている日本語マニュアルでも同様です。
アヒル本の文法は、2017-06-14時点では問題なく新しいバージョンで書かれています。なのでエラーメッセージにうんざりすることもありません。
改訂版に期待する点
紹介した分布を全部使って欲しい
6章で分布を幾つか紹介しています。先述したコーシー分布やディリクレ分布などはここで登場します。
ただ、ここでラプラス分布は紹介されるものの、本中ではそれっきりです。 ぜひこれらの分布についても実例を挙げて解説していただき、活かせる点を見せてくれたらいいなと思いました。きっと何かしらの有用性があると思いますので。
練習問題を面白くして欲しい
4章以降は各章の最後に練習問題が載っています。学習した内容を復習・確認することが出来ますし、答えについてもレポジトリに R + Stanの実装で配布されています。最高です。
ただあくまで確認が目的なのか、他の参考書のデータを引用していたり、本文中の内容の焼き直しだったりしたので退屈に感じました。練習問題に+alphaの内容を盛り込んで、更なる学習が出来ればよりモチベーションアップに繋がるのではないでしょうか。
書評はこんな感じです。 「はじめての統計データ分析」みたいに全部の問題をやったりはしませんが、Python + PyStanで幾つか解いてみたいと思います。タイトルはそういう意味です笑。
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その7
その7です。今回は第6章の章末問題に取り組んでいきます。
6章 比率とクロス表の推測
内容
- 離散的な値をとるカウントデータの解析
カテゴリカル分布 : カテゴリカルな分類においてカウントとして得られるデータに対応する分布
- ベルヌイ試行
- 2項分布
- 多項分布
- 各事象が起こる確率が独立であることに注意
1つの2項分布の比率の推測:「きのこの山とたけのこの山の好きな人の比率」など
- 確率の事前分布は無情報であるので、を仮定してMCMCを行う
- (二項分布を仮定した時の)事象が生起する確率を推測することが可能
- 生成量
- オッズ :
- 2つの選択肢のどちらの確率が大きいかを測ることが出来る
- 仮説(基準とした確率値に対して、pが下回る確率など; 即ち)の検証
- オッズ :
1つの多項分布の比率の推測
- 同様に事前分布は無情報であるので、番目のカテゴリの事前分布にを仮定してMCMCを行う
- 全ての確率の和が1となるので、カテゴリが全部でn個なら n個目は として決定される制限がある
- 0からn番目の事象が独立してそれぞれ起こるのであれば、母数が発生する確率は1つ1つの事象をとして である
- 生成量
- カテゴリ間の比較( など)
- 連言命題が正しい確率
- 同様に事前分布は無情報であるので、番目のカテゴリの事前分布にを仮定してMCMCを行う
複数の2項分布の推測(独立したクロス表) : カテゴリが個あって、それぞれで2つの観測結果が有った際にカテゴリ毎に or 全体でどのような差があるかを調べたい場合
- (独立した)2 × 2のクロス表
- 興味があるベルヌイ試行の確率だけ見る
- 男と女でベルヌイ試行の確率が異なる(独立した2項分布が2つある)と考える
- 生成量
- 比率の差: , 性質の比:
- 性質の違いを考えることが出来る
- オッズ比 :
- 正反応は他方の反応の何倍生じやすいか、を表す
- 比率の差: , 性質の比:
- g × 2 のクロス表
- 各f(theta)が得られる確率については何もわからないので、の無情報事前分布を採用
- (独立した)2 × 2のクロス表
対応あるクロス表の推測(構造のある多項分布): 1つの標本から2回測定が有った場合
- 2 × 2のクロス表: 2つの事象の発生にどれだけ関係があるか?
- 2つの分布が独立ではなく、対応がある場合
- あるいはベルヌイ試行の興味が2つの分布の同時性にある場合
- ;
- 独立と連関
- 変数Aのカテゴリ, 変数Bのカテゴリが独立である場合、
- ピアソン残差 :
- セルごとに計算され、独立性/連関性を示す
- 独立である場合0となる
- 正の値は独立の場合より観測されやすく、負の場合は独立の場合より観測されづらいことを示す
- クラメルの連関係数 : ;はカテゴリの総数
- クロス表全体での連関の大きさを示す
- 小さいほど独立、大きいほど連関である
- 2 × 2のクロス表: 2つの事象の発生にどれだけ関係があるか?
章末問題
前準備
import os from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger logger = get_logger() # 出力用ディレクトリの用意 result_dir = os.path.join("result", "chapter6") if os.path.exists(result_dir) is False: os.makedirs(result_dir) logger.info("{0} is made".format(result_dir))
2項分布によるコインの予測
大学入試みたいに表と裏のでる確率が異なるイカサマコインではなく、普通の等確率なコインを想定しました。
from scipy.stats import binom prob = binom.pmf(3, 5, 0.5) logger.info("コインを5回投げて3回表が出る確率は{0}".format(prob))
多項分布によるじゃんけんの予測
from scipy.stats import multinomial prob = multinomial.pmf([2,2,1], 5, [0.3, 0.3, 0.4]) logger.info("グー:チョキ:パー = 2:2:1となる確率は{0}".format(prob))
MCMCによるベルヌイ試行の確率推定
Stanによる統計モデルの構築
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "binomial.stan") stan_file_c = os.path.join("stan", "binomial.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) logger.info("Stan model is compiled to {0}".format(stan_file_c))
- binomial.stan
data { int<lower=0> g ; int<lower=0> x[g] ; int<lower=0> n[g] ; } parameters { real<lower=0, upper=1> p[g] ; } model { for(i in 1:g){ x[i] ~ binomial(n[i], p[i]) ; } } generated quantities { vector[g] log_lik ; vector[g] prob_p_upper_05 ; for(i in 1:g){ prob_p_upper_05[i] = p[i] > 0.5 ? 1 : 0 ; log_lik[i] = binomial_lpmf(x[i] | n[i], p[i]) ; } }
統計モデルによる事後分布のサンプリング
from IPython.core.display import display import numpy as np import matplotlib import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 20) # Stanで使うデータの用意 stan_data = {"g": 1, "x": np.array([55]), "n": np.array([100])} # 興味のあるパラメータの設定 par = ["p", "prob_p_upper_05", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "binomial.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_1.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_1.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と図)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
p[0] | 0.548981 | 0.000259346 | 0.0491617 | 0.451546 | 0.467666 | 0.515895 | 0.549386 | 0.582306 | 0.629074 | 0.644074 | 35933 | 1.00001 |
prob_p_upper_05[0] | 0.83937 | 0.00175169 | 0.367191 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 43941 | 1.00003 |
log_lik[0] | -3.02129 | 0.0034005 | 0.708006 | -5.02757 | -4.429 | -3.17868 | -2.74975 | -2.57516 | -2.52779 | -2.52636 | 43350 | 1.00014 |
比率の差の推測
相談相手問題(1つの2項分布)
Stanによる統計モデルの構築
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "multinomial.stan") stan_file_c = os.path.join("stan", "multinomial.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) logger.info("Stan model is compiled to {0}".format(stan_file_c))
- multinomial.stan
data { int<lower=0> a; int<lower=0> b; int<lower=0> x[a, b] ; } transformed data { # セル数の和 int<lower=0> ab ; # データのベクトル形式 int<lower=0> xy[a*b] ; # 値の和 int<lower=0> N ; ab = a * b ; for(i in 1:a){ for(j in 1:b){ xy[(i - 1) * b + j] = x[i, j] ; } } # sum関数はint[,]では計算できないので、ベクトル形式に変換が必要 N = sum(xy) ; } parameters { # simplex型は和が1となる simplex[ab] p ; } transformed parameters { # simplex型はベクトルにしか使えないので、matrix形式に変換する real<lower=0, upper=1> pm[a, b] ; for(i in 1:a){ for(j in 1:b){ pm[i, j] = p[(i - 1) * b + j] ; } } } model { xy ~ multinomial(p) ; } generated quantities{ int<lower=1, upper=ab> index ; vector[ab] log_lik ; real pa[a] ; real pb[b] ; real pr[ab] ; real pr2[ab] ; real cac ; # 周辺確率の計算 for(i in 1:a){ pa[i] = sum(pm[i, :]) ; } for(i in 1:b){ pb[i] = sum(pm[:, i]) ; } for(i in 1:a){ for(j in 1:b){ index = (i - 1) * b + j ; pr[index] = (pm[i, j] - pa[i] * pb[j]) / sqrt(pa[i] * pb[j]) ; pr2[index] = pow(pr[index], 2) ; log_lik[index] = multinomial_lpmf(xy | p) ; } } cac = sqrt(sum(pr2) / (min(a, b) - 1)) ; }
統計モデルによる事後分布のサンプリング
from IPython.core.display import display import numpy as np import matplotlib import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 10) # Stanで使うデータの用意 stan_data = {"a": 1, "b": 6, "x": np.array([[30,12,4,20,22,8]])} # 興味のあるパラメータの設定 par = ["p", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "multinomial.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_2.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_2.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と図)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
p[0] | 0.303771 | 0.000143133 | 0.0452625 | 0.219397 | 0.231827 | 0.272348 | 0.302484 | 0.333784 | 0.380652 | 0.395746 | 100000 | 0.999974 |
p[1] | 0.127472 | 0.000103832 | 0.0328345 | 0.0703212 | 0.077966 | 0.104117 | 0.125192 | 0.148024 | 0.185467 | 0.198181 | 100000 | 0.999978 |
p[2] | 0.0490361 | 6.7299e-05 | 0.0212818 | 0.0162597 | 0.019695 | 0.0335044 | 0.0460529 | 0.0616122 | 0.0882118 | 0.098237 | 100000 | 0.999981 |
p[3] | 0.206034 | 0.000126221 | 0.0399147 | 0.133482 | 0.143911 | 0.178019 | 0.204084 | 0.231992 | 0.275117 | 0.289446 | 100000 | 0.999987 |
p[4] | 0.225439 | 0.000129709 | 0.0410175 | 0.150289 | 0.161048 | 0.196605 | 0.223729 | 0.252287 | 0.295423 | 0.310788 | 100000 | 0.999972 |
p[5] | 0.0882487 | 8.84886e-05 | 0.0279826 | 0.0416277 | 0.0471727 | 0.068063 | 0.0855951 | 0.105674 | 0.138238 | 0.149799 | 100000 | 1.00002 |
log_lik[0] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
log_lik[1] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
log_lik[2] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
log_lik[3] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
log_lik[4] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
log_lik[5] | -12.4762 | 0.00655559 | 1.51627 | -16.25 | -15.3943 | -13.2568 | -12.1647 | -11.3654 | -10.6288 | -10.4778 | 53497 | 1 |
治療法問題(2つの2項分布)
Stanによる統計モデルの構築
- binomial.stanをそのまま利用
統計モデルによる事後分布のサンプリング
from IPython.core.display import display import numpy as np import matplotlib import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 10) # Stanで使うデータの用意 stan_data = {"g": 2, "x": np.array([90, 78]), "n": np.array([115, 117])} # 興味のあるパラメータの設定 par = ["p", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "binomial.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_3.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_3.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と図)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
p[0] | 0.777804 | 0.000131177 | 0.0381837 | 0.698637 | 0.712412 | 0.752983 | 0.779475 | 0.804413 | 0.837762 | 0.847901 | 84731 | 1.00002 |
p[1] | 0.663806 | 0.000142741 | 0.0431943 | 0.576869 | 0.591607 | 0.635059 | 0.664743 | 0.693503 | 0.733634 | 0.746176 | 91571 | 1.00004 |
log_lik[0] | -2.9017 | 0.00319529 | 0.699761 | -4.89146 | -4.30242 | -3.06149 | -2.63201 | -2.45849 | -2.41124 | -2.40982 | 47960 | 1.00004 |
log_lik[1] | -3.04626 | 0.00319586 | 0.700447 | -5.04176 | -4.45784 | -3.20754 | -2.77576 | -2.60114 | -2.55245 | -2.55097 | 48037 | 0.999998 |
連関の推測
広告効果(2×2のクロス表)
Stanによる統計モデルの構築
- multinomial.stanをそのまま利用
統計モデルによる事後分布のサンプリング
from IPython.core.display import display import numpy as np import matplotlib import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 20) # Stanで使うデータの用意 stan_data = {"a": 2, "b": 2, "x": np.array([[33, 18], [84, 23]])} # 興味のあるパラメータの設定 par = ["p", "pr", "cac", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "multinomial.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_4.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_4.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と図)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
p[0] | 0.209943 | 0.000100962 | 0.031927 | 0.151078 | 0.159612 | 0.187721 | 0.208703 | 0.230906 | 0.264185 | 0.275663 | 100000 | 0.999996 |
p[1] | 0.117181 | 7.95044e-05 | 0.0251415 | 0.072712 | 0.0787064 | 0.0994485 | 0.115658 | 0.133121 | 0.161086 | 0.170502 | 100000 | 0.999968 |
p[2] | 0.524635 | 0.000123124 | 0.0389353 | 0.448127 | 0.46035 | 0.498185 | 0.524872 | 0.550989 | 0.588352 | 0.600344 | 100000 | 0.999988 |
p[3] | 0.148241 | 8.82448e-05 | 0.0279055 | 0.0976095 | 0.104724 | 0.128717 | 0.146867 | 0.166217 | 0.196445 | 0.206956 | 100000 | 1.00002 |
pr[0] | -0.0618101 | 0.000109123 | 0.0345078 | -0.13111 | -0.119605 | -0.0846981 | -0.061201 | -0.0382136 | -0.00607329 | 0.00411544 | 100000 | 0.999975 |
pr[1] | 0.102741 | 0.000178103 | 0.0563212 | -0.00695254 | 0.0102769 | 0.0642157 | 0.102485 | 0.141141 | 0.195901 | 0.214191 | 100000 | 0.999981 |
pr[2] | 0.0431673 | 7.75096e-05 | 0.0245107 | -0.0028464 | 0.00422608 | 0.0263984 | 0.0423908 | 0.0590213 | 0.0846067 | 0.0937759 | 100000 | 0.999969 |
pr[3] | -0.0716558 | 0.000125491 | 0.0396838 | -0.151139 | -0.137544 | -0.0983243 | -0.0712307 | -0.0444078 | -0.00712286 | 0.00476227 | 100000 | 0.999974 |
cac | 0.148376 | 0.00024284 | 0.0767926 | 0.0130583 | 0.0253949 | 0.0916138 | 0.145924 | 0.200898 | 0.279509 | 0.305541 | 100000 | 0.999981 |
log_lik[0] | -8.68599 | 0.00539516 | 1.20283 | -11.8121 | -11.0697 | -9.22119 | -8.38001 | -7.80783 | -7.38651 | -7.31759 | 49705 | 1.00009 |
log_lik[1] | -8.68599 | 0.00539516 | 1.20283 | -11.8121 | -11.0697 | -9.22119 | -8.38001 | -7.80783 | -7.38651 | -7.31759 | 49705 | 1.00009 |
log_lik[2] | -8.68599 | 0.00539516 | 1.20283 | -11.8121 | -11.0697 | -9.22119 | -8.38001 | -7.80783 | -7.38651 | -7.31759 | 49705 | 1.00009 |
log_lik[3] | -8.68599 | 0.00539516 | 1.20283 | -11.8121 | -11.0697 | -9.22119 | -8.38001 | -7.80783 | -7.38651 | -7.31759 | 49705 | 1.00009 |
ワイン(3×3のクロス表)
Stanによる統計モデルの構築
- multinomial.stanをそのまま利用
統計モデルによる事後分布のサンプリング
from IPython.core.display import display import numpy as np import matplotlib import matplotlib.pyplot as plt import os import pandas as pd import pickle import pystan from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 20) # Stanで使うデータの用意 stan_data = {"a": 3, "b": 3, "x": np.array([[13, 6, 21], [7, 17, 7], [13, 6, 13]])} # 興味のあるパラメータの設定 par = ["p", "pr", "cac", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "multinomial.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_5.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_5.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と図)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
p[0] | 0.125075 | 9.88918e-05 | 0.0312723 | 0.0706493 | 0.0778646 | 0.102857 | 0.12275 | 0.14495 | 0.180202 | 0.19217 | 100000 | 0.999969 |
p[1] | 0.0625879 | 7.2072e-05 | 0.0227912 | 0.0257446 | 0.0299953 | 0.0460044 | 0.0600517 | 0.076303 | 0.104002 | 0.114101 | 100000 | 0.999979 |
p[2] | 0.1964 | 0.000118333 | 0.0374201 | 0.128527 | 0.137895 | 0.170376 | 0.194506 | 0.220619 | 0.260754 | 0.27471 | 100000 | 0.999964 |
p[3] | 0.071492 | 7.68005e-05 | 0.0242864 | 0.0314877 | 0.0362726 | 0.053982 | 0.0688749 | 0.0863221 | 0.115217 | 0.12588 | 100000 | 0.99997 |
p[4] | 0.160706 | 0.000109287 | 0.0345596 | 0.0990021 | 0.107436 | 0.136436 | 0.158562 | 0.182849 | 0.220719 | 0.233882 | 100000 | 0.999966 |
p[5] | 0.0713786 | 7.66519e-05 | 0.0242395 | 0.0317293 | 0.0362675 | 0.0539307 | 0.0688132 | 0.086041 | 0.115415 | 0.125865 | 100000 | 0.999966 |
p[6] | 0.124919 | 9.84871e-05 | 0.0311444 | 0.0707794 | 0.0779037 | 0.102729 | 0.122671 | 0.144522 | 0.179687 | 0.191926 | 100000 | 1 |
p[7] | 0.0625321 | 7.237e-05 | 0.0228854 | 0.025732 | 0.029808 | 0.0459515 | 0.0599421 | 0.0762803 | 0.104219 | 0.114268 | 100000 | 0.999984 |
p[8] | 0.12491 | 9.85523e-05 | 0.031165 | 0.0709911 | 0.0779175 | 0.102493 | 0.122675 | 0.144806 | 0.179601 | 0.191908 | 100000 | 0.999981 |
pr[0] | 0.00455472 | 0.000191147 | 0.0604459 | -0.111419 | -0.0932988 | -0.0370003 | 0.00381112 | 0.0454738 | 0.104957 | 0.123828 | 100000 | 0.99997 |
pr[1] | -0.141598 | 0.000173091 | 0.0547361 | -0.240132 | -0.226578 | -0.180169 | -0.144331 | -0.106119 | -0.046784 | -0.026671 | 100000 | 0.999976 |
pr[2] | 0.117052 | 0.000180329 | 0.0570251 | 0.00421428 | 0.0226175 | 0.0784596 | 0.117427 | 0.155714 | 0.210575 | 0.22766 | 100000 | 0.999961 |
pr[3] | -0.0830821 | 0.00018874 | 0.059685 | -0.19199 | -0.176676 | -0.125273 | -0.0857259 | -0.0437753 | 0.019607 | 0.0406338 | 100000 | 0.99999 |
pr[4] | 0.250214 | 0.000209959 | 0.066395 | 0.116994 | 0.139032 | 0.20578 | 0.251421 | 0.295633 | 0.357572 | 0.377051 | 100000 | 0.999979 |
pr[5] | -0.137782 | 0.000173256 | 0.0547885 | -0.237001 | -0.223096 | -0.176572 | -0.140576 | -0.102 | -0.0427316 | -0.0240325 | 100000 | 0.999979 |
pr[6] | 0.077024 | 0.000205633 | 0.065027 | -0.048242 | -0.0291772 | 0.0322187 | 0.0765335 | 0.121479 | 0.184749 | 0.204917 | 100000 | 1.00001 |
pr[7] | -0.0890021 | 0.000190123 | 0.0601222 | -0.197285 | -0.182537 | -0.131403 | -0.0920298 | -0.0497219 | 0.0150218 | 0.0371228 | 100000 | 0.999977 |
pr[8] | 0.00640808 | 0.000190191 | 0.0601436 | -0.108758 | -0.0916682 | -0.0350148 | 0.0061111 | 0.0472022 | 0.106246 | 0.125763 | 100000 | 0.999969 |
cac | 0.283669 | 0.000194402 | 0.0614752 | 0.162228 | 0.181986 | 0.242123 | 0.28426 | 0.325619 | 0.384129 | 0.402499 | 100000 | 0.999977 |
log_lik[0] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[1] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[2] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[3] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[4] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[5] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[6] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[7] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
log_lik[8] | -19.3975 | 0.00872255 | 1.87604 | -23.8484 | -22.9175 | -20.4394 | -19.0998 | -18.0235 | -16.9287 | -16.6699 | 46259 | 1.00004 |
備考
個人的には、クロス集計はよく使いますのでここで扱ってくれたのはいいものだと思いました。
- ただクロス集計では次のような問題には対応できない
- A群とB群の人間が居て、それぞれa人とb人存在する
- 各人は好きな色(全部でn種類)の飴を好きなだけピックアップする
- 各人の飴の本数、m本は任意の値をとる
- この時、A群とB群で特定の色の飴が取得されやすいか、されづらいかはどう判定すればよい?
- こういう問題ですと、より階層的なモデルが必要になってきます。ぱっと思いつくのはディリクレ分布を使って、トピックの混合率をモデル化するトピックモデルを使う手法でしょうか。
- ただクロス集計では次のような問題には対応できない
「はじめての統計データ分析」はこれで最終回ですね。大分Stanへの恐怖心もなくなってきましたので、次は「実践ベイズモデリング」とか「StanとRでベイズ統計モデリング」あたりを流していこうかと思います。
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その6
その6です。今回は第5章の章末問題に取り組んでいきます。
5章 実験計画による多群の差の推測
内容
- 分散分析によるF検定へのオルタナ
- 要因から生じる影響の調査
- 要因 : 離散的なカテゴリーによる影響のこと
- 水準: 要因が取りうる様々な状態のこと
- 水準数: 水準のパターン数。以下要因の水準数をで表す。
独立した1要因モデル
- モデル
- 要因の番目の水準における番目の測定値が、番目の水準がもつ平均的な効果と正規分布から生じる乱数による現れる
- 以上を纏めて、事後確率はと表される
- 水準内の測定値は独立に観測されてると考える
- 同時確率は
- 即ち、尤度は
- モデル
生成量
モデルを使った要因の効果の方法例
独立した2要因モデル
- 要因を複数導入
- 水準の組み合わせ1単位のことをセルと以下では呼称する
- 交互作用: 一方の要因の水準毎に、他方の要因の平均の高低パターンが異なる現象
- モデル
- 水準の番目の測定値が、全体の平均に要因の番目の水準の効果 + 要因の番目の水準の効果 + との交互作用の効果 + 乱数から構成されることを示す
- セルの平均 + 乱数
- それぞれ全平均(水準の効果)の性質から次の制約を持つ(水準の効果を全水準で足すと0になる性質を定式化)
- 尤度
- 要因を複数導入
章末問題
1要因の推測
データの取得
import numpy as np import os from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger logger = get_logger() # データの準備 ld = np.array([5.02, 6.67, 8.17, 2.79, 8.13, 6.34, 6.32, 3.97]) ll = np.array([9.89, 9.58, 11.20, 9.05, 12.33, 9.39, 10.88, 9.37, 17.40]) dm = np.array([10.20, 7.29, 7.57, 3.42, 5.82, 10.92, 5.21, 13.47, 8.64, 6.05]) y = np.concatenate([a, b, c]) A_ld = np.full(ld.size, 1, np.int) A_ll = np.full(ll.size, 2, np.int) A_dm = np.full(dm.size, 3, np.int) A = np.concatenate([A_ld, A_ll, A_dm]) # 出力用ディレクトリの用意 result_dir = os.path.join("result", "chapter5") if os.path.exists(result_dir) is False: os.makedirs(result_dir) logger.info("{0} is made".format(result_dir))
Stanによる統計モデルの構築
特に難しいところは無いような気がします。説明量を計算するために、分散値を計算する手間があるぐらいでしょうか?
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "e1_ind.stan") stan_file_c = os.path.join("stan", "e1_ind.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) logger.info("Stan model is compiled to {0}".format(stan_file_c))
e1_ind.stan
data { int<lower=0> n ; int<lower=0> a ; real<lower=0> y[n] ; int<lower=1, upper=3> A[n] ; } parameters { vector[a] mu_a ; real<lower=0> sigma_e ; } model { for(i in 1:n){ y[i] ~ normal(mu_a[A[i]], sigma_e) ; } } generated quantities { real mu ; vector[a] eol ; vector[a] eol2 ; vector[a] prob_eol_upper_0 ; real var_e ; real var_a ; real eta2 ; real delta ; int<lower=0, upper=1> prob_ll_upper_dm ; int<lower=0, upper=1> prob_ll_upper_ld ; int<lower=0, upper=1> prob_dm_upper_ld ; vector[n] log_lik ; mu = mean(mu_a) ; for(i in 1:a){ eol[i] = mu_a[i] - mu ; eol2[i] = pow(eol[i], 2) ; prob_eol_upper_0[i] = eol[i] > 0 ? 1 : 0 ; } var_a = mean(eol2) ; var_e = pow(sigma_e, 2) ; eta2 = var_a / (var_a + var_e) ; delta = sqrt(var_a) / sigma_e ; prob_ll_upper_dm = mu_a[2] - mu_a[3] > 0 ? 1 : 0 ; prob_ll_upper_ld = mu_a[2] - mu_a[1] > 0 ? 1 : 0 ; prob_dm_upper_ld = mu_a[3] - mu_a[1] > 0 ? 1 : 0 ; for(i in 1:n){ log_lik[i] = normal_lpdf(y[i] | mu_a[A[i]], sigma_e) ; } }
統計モデルによる事後分布のサンプリング
import pandas as pd import pickle import pystan import matplotlib import matplotlib.pyplot as plt import os from IPython.core.display import display from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 50) # Stanで使うデータの用意 stan_data = {"n": y.size, "a": 3, "y": y, "A": A} # 興味のあるパラメータの設定 par = ["mu_a", "sigma_e", "eol", "eta2", "delta", "prob_eol_upper_0", "prob_ll_upper_dm", "prob_ll_upper_ld", "prob_dm_upper_ld", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "e1_ind.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_1.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot_1.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu_a[0] | 5.92465 | 0.0030795 | 0.973823 | 3.99253 | 4.32772 | 5.28835 | 5.92998 | 6.56208 | 7.51919 | 7.84787 | 100000 | 1.00007 |
mu_a[1] | 11.0072 | 0.00292346 | 0.924478 | 9.17396 | 9.49254 | 10.4026 | 11.0076 | 11.6092 | 12.5191 | 12.8357 | 100000 | 0.999985 |
mu_a[2] | 7.85955 | 0.00276635 | 0.874797 | 6.12594 | 6.42644 | 7.28474 | 7.86113 | 8.43353 | 9.29189 | 9.58305 | 100000 | 1 |
sigma_e | 2.74451 | 0.00151073 | 0.42629 | 2.06146 | 2.14526 | 2.44184 | 2.69385 | 2.9907 | 3.51953 | 3.72462 | 79623 | 0.999968 |
eol[0] | -2.33916 | 0.00245045 | 0.7749 | -3.88038 | -3.61524 | -2.8447 | -2.33923 | -1.83011 | -1.07629 | -0.817259 | 100000 | 1.00008 |
eol[1] | 2.74342 | 0.00239336 | 0.756846 | 1.24128 | 1.50518 | 2.24837 | 2.74077 | 3.23758 | 3.98729 | 4.24112 | 100000 | 0.999999 |
eol[2] | -0.404261 | 0.00233393 | 0.738053 | -1.86101 | -1.60891 | -0.887243 | -0.405952 | 0.0805197 | 0.804455 | 1.05156 | 100000 | 1.00003 |
eta2 | 0.384362 | 0.000399327 | 0.126278 | 0.122444 | 0.162176 | 0.299146 | 0.391642 | 0.476528 | 0.579842 | 0.608479 | 100000 | 0.999981 |
delta | 0.804616 | 0.000704594 | 0.222812 | 0.373535 | 0.439964 | 0.653323 | 0.802353 | 0.954108 | 1.17476 | 1.24665 | 100000 | 0.999984 |
prob_eol_upper_0[0] | 0.00215 | 0.000165199 | 0.0463185 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 78613 | 1.00006 |
prob_eol_upper_0[1] | 0.99963 | 6.2415e-05 | 0.0192319 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 94944 | 1.00003 |
prob_eol_upper_0[2] | 0.28713 | 0.00152318 | 0.452425 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 88225 | 1.00002 |
prob_ll_upper_dm | 0.99168 | 0.000333967 | 0.0908342 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 73976 | 0.999978 |
prob_ll_upper_ld | 0.99975 | 5.19414e-05 | 0.0158095 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 92642 | 0.999991 |
prob_dm_upper_ld | 0.93277 | 0.000930336 | 0.250421 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 72454 | 1.00007 |
log_lik[0] | -2.03678 | 0.000758605 | 0.199372 | -2.49373 | -2.39575 | -2.14793 | -2.01548 | -1.89931 | -1.75258 | -1.7086 | 69071 | 0.999994 |
log_lik[1] | -2.01783 | 0.000754794 | 0.191828 | -2.4532 | -2.35815 | -2.12696 | -1.99886 | -1.88578 | -1.74113 | -1.69764 | 64590 | 1.00002 |
log_lik[2] | -2.33612 | 0.000992587 | 0.313883 | -3.08375 | -2.92534 | -2.50565 | -2.2841 | -2.11391 | -1.92039 | -1.86752 | 100000 | 1.00004 |
log_lik[3] | -2.67679 | 0.00135735 | 0.429231 | -3.69657 | -3.47937 | -2.91461 | -2.60959 | -2.36613 | -2.10742 | -2.04192 | 100000 | 1.00004 |
log_lik[4] | -2.32348 | 0.000978475 | 0.309421 | -3.06092 | -2.90399 | -2.49026 | -2.27211 | -2.10467 | -1.91389 | -1.86127 | 100000 | 1.00004 |
log_lik[5] | -1.99069 | 0.000735103 | 0.179753 | -2.39033 | -2.3073 | -2.0972 | -1.97534 | -1.86648 | -1.72545 | -1.6824 | 59794 | 1.00001 |
log_lik[6] | -1.98954 | 0.00073424 | 0.17923 | -2.38764 | -2.30472 | -2.09579 | -1.9742 | -1.86567 | -1.72457 | -1.6816 | 59586 | 1.00001 |
log_lik[7] | -2.25015 | 0.000893585 | 0.282577 | -2.92974 | -2.78158 | -2.40036 | -2.20673 | -2.05154 | -1.87364 | -1.82393 | 100000 | 1.00003 |
log_lik[8] | -2.06094 | 0.000756641 | 0.203617 | -2.52871 | -2.42581 | -2.17552 | -2.0385 | -1.91954 | -1.77099 | -1.72765 | 72418 | 1 |
log_lik[9] | -2.11696 | 0.000790546 | 0.223723 | -2.63793 | -2.52526 | -2.2399 | -2.08855 | -1.96137 | -1.80487 | -1.76051 | 80088 | 0.999999 |
log_lik[10] | -1.9749 | 0.000705525 | 0.171411 | -2.3497 | -2.27598 | -2.07886 | -1.96221 | -1.85611 | -1.71601 | -1.67376 | 59027 | 1.00001 |
log_lik[11] | -2.24435 | 0.000846443 | 0.267669 | -2.87634 | -2.74126 | -2.38938 | -2.20473 | -2.05568 | -1.88209 | -1.83201 | 100000 | 0.999994 |
log_lik[12] | -2.09647 | 0.000768147 | 0.216744 | -2.60167 | -2.49056 | -2.21654 | -2.06987 | -1.94642 | -1.79251 | -1.74675 | 79617 | 1 |
log_lik[13] | -2.15804 | 0.000816763 | 0.238127 | -2.71684 | -2.59706 | -2.2877 | -2.12586 | -1.9917 | -1.83017 | -1.78421 | 85001 | 0.999997 |
log_lik[14] | -1.97342 | 0.000706486 | 0.170781 | -2.34448 | -2.27222 | -2.07722 | -1.96171 | -1.85505 | -1.71519 | -1.67338 | 58435 | 1.00001 |
log_lik[15] | -2.16266 | 0.000819785 | 0.239732 | -2.72544 | -2.60451 | -2.29303 | -2.13016 | -1.99503 | -1.83272 | -1.78664 | 85517 | 0.999997 |
log_lik[16] | -4.87403 | 0.00339268 | 1.07286 | -7.32982 | -6.8243 | -5.5138 | -4.74535 | -4.09518 | -3.35255 | -3.15931 | 100000 | 0.999977 |
log_lik[17] | -2.35558 | 0.000919091 | 0.290642 | -3.03811 | -2.89536 | -2.51793 | -2.31435 | -2.14877 | -1.95975 | -1.90793 | 100000 | 1 |
log_lik[18] | -1.98969 | 0.000680232 | 0.173859 | -2.37117 | -2.2949 | -2.09513 | -1.97689 | -1.86882 | -1.72784 | -1.68611 | 65325 | 0.999965 |
log_lik[19] | -1.97261 | 0.000674084 | 0.168196 | -2.33564 | -2.26748 | -2.07642 | -1.96134 | -1.85564 | -1.7164 | -1.67546 | 62259 | 0.999965 |
log_lik[20] | -3.36627 | 0.00188208 | 0.595165 | -4.73815 | -4.46265 | -3.7144 | -3.2838 | -2.93321 | -2.54231 | -2.44062 | 100000 | 0.999985 |
log_lik[21] | -2.26206 | 0.000825123 | 0.260927 | -2.87464 | -2.74387 | -2.40867 | -2.22729 | -2.07788 | -1.90249 | -1.85212 | 100000 | 0.999979 |
log_lik[22] | -2.6317 | 0.00118696 | 0.37535 | -3.50938 | -3.32991 | -2.84462 | -2.57612 | -2.36118 | -2.12327 | -2.05933 | 100000 | 1 |
log_lik[23] | -2.46518 | 0.00102346 | 0.323646 | -3.2244 | -3.0644 | -2.64582 | -2.41881 | -2.23303 | -2.02502 | -1.96827 | 100000 | 0.999983 |
log_lik[24] | -4.20173 | 0.00267471 | 0.845818 | -6.14049 | -5.74926 | -4.70266 | -4.09605 | -3.58666 | -3.0153 | -2.86741 | 100000 | 0.99999 |
log_lik[25] | -2.00989 | 0.000689887 | 0.180727 | -2.40828 | -2.32687 | -2.11763 | -1.99477 | -1.88481 | -1.74166 | -1.69933 | 68626 | 0.999981 |
log_lik[26] | -2.19919 | 0.000762851 | 0.241235 | -2.76133 | -2.64207 | -2.33492 | -2.16884 | -2.03002 | -1.86336 | -1.81479 | 100000 | 0.999977 |
2要因の推測
データの取得
データの取得がめちゃめちゃ汚い気がしますが、許してください…
import numpy as np import os from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger logger = get_logger() # データの準備 y = np.array([ 140,146,149,136,147,147,143,143,143,141, 139,136,136,140,135,132,140,134, 123,127,131,130,138,128,129, 115,120,118,118,121,124,129,119,128, 128,124,123,121,122,126,131,122, 121,121,120,116,117,113,118, 143,141,142,145,149,145,143,141,142,155, 138,134,142,136,135,136,131,133, 131,128,128,128,127,130,130, 117,125,132,122,119,122,129,117,127, 117,120,124,122,122,122,118,122, 119,125,122,116,119,113,122]) A_loaded = np.full(49, 1, np.int) A_empty = np.full(y.size - A_loaded.size, 2, np.int) A = np.concatenate([A_loaded, A_empty]) B_loaded_st = np.full(10, 1, np.int) B_loaded_cut = np.full(8, 2, np.int) B_loaded_f = np.full(7, 3, np.int) B_loaded_ch = np.full(9, 4, np.int) B_loaded_sl = np.full(8, 5, np.int) B_loaded_cur = np.full(7, 6, np.int) B_empty_st = np.full(10, 1, np.int) B_empty_cut = np.full(8, 2, np.int) B_empty_f = np.full(7, 3, np.int) B_empty_ch = np.full(9, 4, np.int) B_empty_sl = np.full(8, 5, np.int) B_empty_cur = np.full(7, 6, np.int) B = np.concatenate([B_loaded_st, B_loaded_cut, B_loaded_f, B_loaded_ch, B_loaded_sl, B_loaded_cur, B_empty_st, B_empty_cut, B_empty_f, B_empty_ch, B_empty_sl, B_empty_cur]) # 出力用ディレクトリの用意 result_dir = os.path.join("result", "chapter5") if os.path.exists(result_dir) is False: os.makedirs(result_dir) logger.info("{0} is made".format(result_dir))
Stanによる統計モデルの構築
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "e2_ind.stan") stan_file_c = os.path.join("stan", "e2_ind.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) logger.info("Stan model is compiled to {0}".format(stan_file_c))
e2_ind.stan
- PyStanがmatrixを可視化する際にうまくうけとってくれないので、ベクトルに変換してからPythonへ渡しています
- 本の付録のStanコードではtransformed parametersの部分が冗長だったので、forとif文で対応しています。
data { int<lower=1> n ; int<lower=1> n_a ; int<lower=1> n_b ; int<lower=0> y[n] ; int<lower=1, upper=n_a> A[n] ; int<lower=1, upper=n_b> B[n] ; } parameters { real<lower=0> mu ; vector[n_a-1] a ; vector[n_b-1] b ; matrix[n_a-1, n_b-1] ab ; real<lower=0> sigma_e ; } transformed parameters { # 水準の効果の和が0になる制約条件を使う vector[n_a] a_r ; vector[n_b] b_r ; matrix[n_a, n_b] ab_r ; vector[n_a*n_b] ab_v ; for(i in 1:n_a-1){ a_r[i] = a[i] ; } a_r[n_a] = - sum(a_r[1:n_a - 1]) ; for(i in 1:(n_b-1)){ b_r[i] = b[i] ; } b_r[n_b] = - sum(b_r[1:n_b - 1]) ; for (j in 1:n_a) { for (k in 1:n_b) { if (j == n_a && k == n_b) { ab_r[j, k] = sum(ab[1:j-1, 1:k-1]) ; } else if (j == n_a){ ab_r[j, k] = -sum(ab[1:j-1, k]) ; } else if (k == n_b){ ab_r[j, k] = -sum(ab[j, 1:k-1]) ; } else { ab_r[j, k] = ab[j, k] ; } ab_v[(j-1) * n_b + k] = ab_r[j, k] ; } } } model { for (i in 1:n) { y[i] ~ normal(mu + a_r[A[i]] + b_r[B[i]] + ab_r[A[i], B[i]], sigma_e) ; } } generated quantities { vector[n] log_lik ; vector[n_b] prob_b_upper_0 ; vector[n_a] a_r2 ; vector[n_b] b_r2 ; vector[n_a*n_b] ab_r2 ; real var_a ; real var_b ; real var_ab ; real var_e ; real var_y ; real eta_a2 ; real eta_b2 ; real eta_ab2 ; real eta_t2 ; real delta_a ; real delta_b ; real delta_ab ; matrix[n_b, n_b] prob_b_comparison ; for (i in 1:n) { log_lik[i] = normal_lpdf(y[i] | mu + a_r[A[i]] + b_r[B[i]] + ab_r[A[i], B[i]], sigma_e) ; } for (i in 1:n_b) { prob_b_upper_0[i] = b_r[i] > 0 ? 1 : 0 ; } for (i in 1:n_a) { a_r2[i] = pow(a_r[i], 2) ; } for (i in 1:n_b) { b_r2[i] = pow(b_r[i], 2) ; } for (i in 1:n_a*n_b) { ab_r2[i] = pow(ab_v[i], 2) ; } var_a = mean(a_r2) ; var_b = mean(b_r2) ; var_ab = mean(ab_r2) ; var_e = pow(sigma_e, 2) ; var_y = var_a + var_b + var_ab + var_e ; eta_a2 = var_a / var_y ; eta_b2 = var_b / var_y ; eta_ab2 = var_ab / var_y ; eta_t2 = (var_a + var_b + var_ab) / var_y ; delta_a = sqrt(var_a) / sigma_e ; delta_b = sqrt(var_b) / sigma_e ; delta_ab = sqrt(var_ab) / sigma_e ; for (j in 1:n_b){ for (k in 1:n_b){ prob_b_comparison[j, k] = b_r[j] > b_r[k] ? 1 : 0 ; } } }
統計モデルによる事後分布のサンプリング
- データ型が行列のパラメータは、可視化する前に除いています
- 水準の効果の差については、別途ピボットテーブルを作成して可視化しました。
import pandas as pd import pickle import pystan import matplotlib import matplotlib.pyplot as plt import os from IPython.core.display import display from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 50) # Stanで使うデータの用意 stan_data = {"n": y.size, "n_a": np.unique(A).size, "n_b": np.unique(B).size, "y": y, "A": A, "B": B} # 興味のあるパラメータの設定 pars = ["mu", "sigma_e", "a_r", "b_r", "ab_v", "prob_b_upper_0", "eta_a2", "eta_b2", "eta_ab2", "eta_t2", "delta_a", "delta_b", "delta_ab", "prob_b_comparison", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "e2_ind.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=pars, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=pars, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary_2.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 水準毎に上回る確率を調べる comparison = fit.summary(pars=["prob_b_comparison"], probs=prob) comparison_df = pd.DataFrame(comparison["summary"], index=comparison["summary_rownames"], columns=comparison["summary_colnames"]) comparison_index = pd.Series(comparison_df.index.str.strip("prob_b_comparison[").str.strip("]").str.split(",")) b_dict = {"0":"straight", "1":"cutball", "2":"fork", "3":"changeup", "4":"slider", "5":"curve"} upper_index = comparison_index.apply(lambda x:b_dict[x[0]]) lower_index = comparison_index.apply(lambda x:b_dict[x[1]]) comparison_triple = pd.DataFrame([upper_index, lower_index, comparison_df["mean"].values]).T comparison_triple.columns = ["upper", "lower", "possibility"] comparison_pivot = comparison_triple.pivot(index="upper", columns="lower", values="possibility") display(comparison_pivot) b_df_path = os.path.join(result_dir, "b_comparison.md") with open(b_df_path, "w") as f: f.write(tabulate(comparison_pivot, comparison_pivot.columns, tablefmt="pipe")) logger.info("comparison table is saved at {0}".format(b_df_path)) # 事後分布の可視化 pars.remove("prob_b_upper_0") pars.remove("prob_b_comparison") fit.traceplot(pars) traceplot_path = os.path.join(result_dir, "traceplot_2.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu | 128.842 | 0.00125739 | 0.397622 | 128.059 | 128.188 | 128.577 | 128.842 | 129.107 | 129.495 | 129.622 | 100000 | 0.999982 |
sigma_e | 3.89006 | 0.000955725 | 0.302227 | 3.35328 | 3.42584 | 3.67831 | 3.87124 | 4.08093 | 4.41686 | 4.53412 | 100000 | 1.00008 |
a_r[0] | 0.0566725 | 0.00125659 | 0.39737 | -0.722588 | -0.595841 | -0.208101 | 0.0562403 | 0.321384 | 0.710725 | 0.840036 | 100000 | 0.999989 |
a_r[1] | -0.0566725 | 0.00125659 | 0.39737 | -0.840036 | -0.710725 | -0.321384 | -0.0562403 | 0.208101 | 0.595841 | 0.722588 | 100000 | 0.999989 |
b_r[0] | 15.2089 | 0.00256537 | 0.811242 | 13.6165 | 13.8763 | 14.6635 | 15.2074 | 15.7524 | 16.5397 | 16.8083 | 100000 | 0.999994 |
b_r[1] | 7.22147 | 0.00280979 | 0.888533 | 5.47369 | 5.76311 | 6.62748 | 7.22389 | 7.81345 | 8.68103 | 8.97568 | 100000 | 1 |
b_r[2] | 0.303162 | 0.00295644 | 0.93491 | -1.53609 | -1.2345 | -0.321006 | 0.298546 | 0.928776 | 1.83405 | 2.15155 | 100000 | 0.999977 |
b_r[3] | -6.50771 | 0.00270552 | 0.85556 | -8.17839 | -7.90278 | -7.08083 | -6.51271 | -5.93491 | -5.10512 | -4.82134 | 100000 | 0.999979 |
b_r[4] | -6.09233 | 0.00281005 | 0.888616 | -7.83205 | -7.54945 | -6.68611 | -6.09106 | -5.49716 | -4.62802 | -4.3491 | 100000 | 0.999976 |
b_r[5] | -10.1335 | 0.00296533 | 0.93772 | -11.9753 | -11.6719 | -10.7628 | -10.1342 | -9.49813 | -8.59766 | -8.30235 | 100000 | 0.999964 |
ab_v[0] | -0.607631 | 0.00258641 | 0.817893 | -2.21712 | -1.9559 | -1.15031 | -0.607124 | -0.0654354 | 0.744995 | 1.00949 | 100000 | 0.999985 |
ab_v[1] | 0.381475 | 0.00281664 | 0.8907 | -1.37495 | -1.08112 | -0.212969 | 0.382057 | 0.976252 | 1.84594 | 2.13356 | 100000 | 0.999998 |
ab_v[2] | 0.23221 | 0.00297549 | 0.940934 | -1.60767 | -1.31047 | -0.399151 | 0.229995 | 0.866498 | 1.77931 | 2.08061 | 100000 | 0.999969 |
ab_v[3] | -1.05319 | 0.00268952 | 0.8505 | -2.72909 | -2.45038 | -1.62572 | -1.05216 | -0.482988 | 0.345127 | 0.618895 | 100000 | 0.99998 |
ab_v[4] | 1.81474 | 0.00281049 | 0.888756 | 0.0653521 | 0.350221 | 1.22523 | 1.81277 | 2.40855 | 3.27038 | 3.56064 | 100000 | 0.999986 |
ab_v[5] | -0.767603 | 0.00297959 | 0.94223 | -2.60335 | -2.31877 | -1.40147 | -0.769324 | -0.130914 | 0.778092 | 1.07769 | 100000 | 0.999986 |
ab_v[6] | 0.607631 | 0.00258641 | 0.817893 | -1.00949 | -0.744995 | 0.0654354 | 0.607124 | 1.15031 | 1.9559 | 2.21712 | 100000 | 0.999985 |
ab_v[7] | -0.381475 | 0.00281664 | 0.8907 | -2.13356 | -1.84594 | -0.976252 | -0.382057 | 0.212969 | 1.08112 | 1.37495 | 100000 | 0.999998 |
ab_v[8] | -0.23221 | 0.00297549 | 0.940934 | -2.08061 | -1.77931 | -0.866498 | -0.229995 | 0.399151 | 1.31047 | 1.60767 | 100000 | 0.999969 |
ab_v[9] | 1.05319 | 0.00268952 | 0.8505 | -0.618895 | -0.345127 | 0.482988 | 1.05216 | 1.62572 | 2.45038 | 2.72909 | 100000 | 0.99998 |
ab_v[10] | -1.81474 | 0.00281049 | 0.888756 | -3.56064 | -3.27038 | -2.40855 | -1.81277 | -1.22523 | -0.350221 | -0.0653521 | 100000 | 0.999986 |
ab_v[11] | 0.767603 | 0.00297959 | 0.94223 | -1.07769 | -0.778092 | 0.130914 | 0.769324 | 1.40147 | 2.31877 | 2.60335 | 100000 | 0.999986 |
prob_b_upper_0[0] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_upper_0[1] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_upper_0[2] | 0.62685 | 0.00152942 | 0.483644 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 100000 | 1 |
prob_b_upper_0[3] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_upper_0[4] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_upper_0[5] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
eta_a2 | 0.00168409 | 1.10406e-05 | 0.00242859 | 1.69709e-06 | 6.70442e-06 | 0.000168269 | 0.000753203 | 0.00220135 | 0.00650296 | 0.00855088 | 48386 | 1.00007 |
eta_b2 | 0.820698 | 8.31697e-05 | 0.0263006 | 0.762582 | 0.77391 | 0.804681 | 0.822937 | 0.83927 | 0.859708 | 0.865594 | 100000 | 1.00006 |
eta_ab2 | 0.0179692 | 3.28066e-05 | 0.00948945 | 0.00401563 | 0.00533333 | 0.0109958 | 0.0164943 | 0.0233222 | 0.0357468 | 0.0404685 | 83668 | 1.00012 |
eta_t2 | 0.840352 | 7.44251e-05 | 0.0235353 | 0.788162 | 0.798497 | 0.826116 | 0.842443 | 0.857014 | 0.875117 | 0.880387 | 100000 | 1.00003 |
delta_a | 0.0820056 | 0.000277709 | 0.0621218 | 0.0032616 | 0.00655807 | 0.032731 | 0.0691728 | 0.117879 | 0.201804 | 0.231203 | 50039 | 1.00011 |
delta_b | 2.28756 | 0.000631375 | 0.199658 | 1.90128 | 1.96222 | 2.15243 | 2.28528 | 2.42122 | 2.62019 | 2.68599 | 100000 | 1.00004 |
delta_ab | 0.325882 | 0.000285382 | 0.0902456 | 0.157327 | 0.181951 | 0.262588 | 0.3238 | 0.386283 | 0.478362 | 0.50898 | 100000 | 1.00008 |
prob_b_comparison[0,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[1,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[2,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[3,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[4,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[5,0] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[0,1] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[1,1] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[2,1] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[3,1] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[4,1] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[5,1] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[0,2] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[1,2] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[2,2] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[3,2] | 1e-05 | 1e-05 | 0.00316228 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | 1 |
prob_b_comparison[4,2] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[5,2] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[0,3] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[1,3] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[2,3] | 0.99999 | 1e-05 | 0.00316228 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | 1 |
prob_b_comparison[3,3] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[4,3] | 0.62121 | 0.00153398 | 0.485088 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 100000 | 0.999983 |
prob_b_comparison[5,3] | 0.00526 | 0.000240106 | 0.0723352 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 90760 | 0.999992 |
prob_b_comparison[0,4] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[1,4] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[2,4] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[3,4] | 0.37879 | 0.00153398 | 0.485088 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 100000 | 0.999983 |
prob_b_comparison[4,4] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_b_comparison[5,4] | 0.00256 | 0.000167342 | 0.0505319 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 91184 | 1 |
prob_b_comparison[0,5] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[1,5] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[2,5] | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_b_comparison[3,5] | 0.99474 | 0.000240106 | 0.0723352 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 90760 | 0.999992 |
prob_b_comparison[4,5] | 0.99744 | 0.000167342 | 0.0505319 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 91184 | 1 |
prob_b_comparison[5,5] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
log_lik[0] | -2.73633 | 0.000935194 | 0.295734 | -3.43849 | -3.28893 | -2.90328 | -2.68784 | -2.5182 | -2.34766 | -2.30554 | 100000 | 1.00002 |
log_lik[1] | -2.53453 | 0.00070078 | 0.221606 | -3.08005 | -2.95548 | -2.65052 | -2.49078 | -2.37503 | -2.25745 | -2.22527 | 100000 | 0.999975 |
log_lik[2] | -3.34177 | 0.00146612 | 0.463629 | -4.39499 | -4.1848 | -3.6203 | -3.28603 | -3.00338 | -2.6859 | -2.6025 | 100000 | 0.999994 |
log_lik[3] | -4.21624 | 0.00207097 | 0.654899 | -5.66336 | -5.38642 | -4.62001 | -4.15503 | -3.74499 | -3.25356 | -3.11148 | 100000 | 0.999974 |
log_lik[4] | -2.73634 | 0.000936502 | 0.296148 | -3.44194 | -3.28963 | -2.90404 | -2.68719 | -2.51799 | -2.34749 | -2.30663 | 100000 | 0.999975 |
log_lik[5] | -2.73634 | 0.000936502 | 0.296148 | -3.44194 | -3.28963 | -2.90404 | -2.68719 | -2.51799 | -2.34749 | -2.30663 | 100000 | 0.999975 |
log_lik[6] | -2.33272 | 0.000426818 | 0.111143 | -2.59392 | -2.52945 | -2.39027 | -2.31979 | -2.25852 | -2.17872 | -2.15482 | 67807 | 1.0002 |
log_lik[7] | -2.33272 | 0.000426818 | 0.111143 | -2.59392 | -2.52945 | -2.39027 | -2.31979 | -2.25852 | -2.17872 | -2.15482 | 67807 | 1.0002 |
log_lik[8] | -2.33272 | 0.000426818 | 0.111143 | -2.59392 | -2.52945 | -2.39027 | -2.31979 | -2.25852 | -2.17872 | -2.15482 | 67807 | 1.0002 |
log_lik[9] | -2.53452 | 0.00069946 | 0.221189 | -3.07629 | -2.95565 | -2.65031 | -2.49165 | -2.37471 | -2.25725 | -2.22585 | 100000 | 1.00006 |
log_lik[10] | -2.54659 | 0.000785481 | 0.248391 | -3.16542 | -3.02514 | -2.67084 | -2.4917 | -2.36823 | -2.24892 | -2.2172 | 100000 | 1.00005 |
log_lik[11] | -2.34498 | 0.000496465 | 0.125425 | -2.64993 | -2.57005 | -2.40342 | -2.3269 | -2.26265 | -2.18144 | -2.15687 | 63825 | 1.00009 |
log_lik[12] | -2.34498 | 0.000496465 | 0.125425 | -2.64993 | -2.57005 | -2.40342 | -2.3269 | -2.26265 | -2.18144 | -2.15687 | 63825 | 1.00009 |
log_lik[13] | -2.74833 | 0.00104876 | 0.331647 | -3.546 | -3.37344 | -2.93064 | -2.68727 | -2.50315 | -2.32712 | -2.28515 | 100000 | 1.00003 |
log_lik[14] | -2.41232 | 0.000613607 | 0.176642 | -2.86322 | -2.75229 | -2.4868 | -2.3749 | -2.29388 | -2.20328 | -2.1765 | 82872 | 1.00003 |
log_lik[15] | -3.01794 | 0.00133705 | 0.422813 | -4.00776 | -3.80561 | -3.25894 | -2.95459 | -2.70599 | -2.4528 | -2.39164 | 100000 | 1.00002 |
log_lik[16] | -2.74833 | 0.00104876 | 0.331647 | -3.546 | -3.37344 | -2.93064 | -2.68727 | -2.50315 | -2.32712 | -2.28515 | 100000 | 1.00003 |
log_lik[17] | -2.54692 | 0.000791085 | 0.250163 | -3.16951 | -3.02844 | -2.67035 | -2.49157 | -2.36771 | -2.24974 | -2.21829 | 100000 | 1.00002 |
log_lik[18] | -3.73765 | 0.00205657 | 0.650345 | -5.21599 | -4.92322 | -4.1275 | -3.66109 | -3.2632 | -2.81517 | -2.69676 | 100000 | 0.999991 |
log_lik[19] | -2.5447 | 0.000828546 | 0.262009 | -3.21331 | -3.05791 | -2.67094 | -2.48239 | -2.35768 | -2.24115 | -2.21039 | 100000 | 1.00001 |
log_lik[20] | -2.42805 | 0.000693641 | 0.195391 | -2.93941 | -2.809 | -2.50779 | -2.38289 | -2.29826 | -2.2048 | -2.17837 | 79349 | 1.00005 |
log_lik[21] | -2.3563 | 0.000560908 | 0.138402 | -2.70641 | -2.61165 | -2.4159 | -2.33316 | -2.26661 | -2.1839 | -2.15899 | 60884 | 1.00014 |
log_lik[22] | -4.81377 | 0.00283885 | 0.897723 | -6.79217 | -6.42019 | -5.36376 | -4.72578 | -4.16988 | -3.49737 | -3.30988 | 100000 | 0.999977 |
log_lik[23] | -2.41463 | 0.000656152 | 0.184738 | -2.90055 | -2.776 | -2.48901 | -2.37315 | -2.29235 | -2.20083 | -2.17451 | 79269 | 1.00006 |
log_lik[24] | -2.35183 | 0.000548798 | 0.133508 | -2.68773 | -2.59816 | -2.41028 | -2.33023 | -2.26506 | -2.18237 | -2.15774 | 59182 | 1.00015 |
log_lik[25] | -3.68089 | 0.00179969 | 0.569113 | -4.95951 | -4.71565 | -4.02176 | -3.61904 | -3.26868 | -2.86285 | -2.7564 | 100000 | 0.99997 |
log_lik[26] | -2.39015 | 0.000531975 | 0.156079 | -2.78711 | -2.69127 | -2.45852 | -2.36012 | -2.28619 | -2.19753 | -2.17168 | 86081 | 1.00001 |
log_lik[27] | -2.70464 | 0.000949753 | 0.300338 | -3.42914 | -3.27577 | -2.86761 | -2.65075 | -2.48359 | -2.32096 | -2.28113 | 100000 | 0.999965 |
log_lik[28] | -2.70464 | 0.000949753 | 0.300338 | -3.42914 | -3.27577 | -2.86761 | -2.65075 | -2.48359 | -2.32096 | -2.28113 | 100000 | 0.999965 |
log_lik[29] | -2.3338 | 0.000448643 | 0.114007 | -2.60405 | -2.53915 | -2.39156 | -2.31964 | -2.25791 | -2.17838 | -2.15378 | 64574 | 1.00008 |
log_lik[30] | -2.56839 | 0.000779171 | 0.246395 | -3.17826 | -3.04225 | -2.69622 | -2.51769 | -2.39041 | -2.2649 | -2.23265 | 100000 | 0.999982 |
log_lik[31] | -4.30475 | 0.00222789 | 0.70452 | -5.86668 | -5.56526 | -4.73925 | -4.23866 | -3.80034 | -3.26551 | -3.12075 | 100000 | 0.999971 |
log_lik[32] | -2.51376 | 0.000705027 | 0.222949 | -3.07066 | -2.94619 | -2.62376 | -2.46574 | -2.35573 | -2.24338 | -2.21269 | 100000 | 0.999976 |
log_lik[33] | -3.82294 | 0.00190039 | 0.600955 | -5.17464 | -4.90752 | -4.18988 | -3.75953 | -3.39001 | -2.94919 | -2.83484 | 100000 | 0.999968 |
log_lik[34] | -2.72106 | 0.00102137 | 0.322986 | -3.50706 | -3.33946 | -2.89458 | -2.66028 | -2.48187 | -2.31426 | -2.27516 | 100000 | 0.99998 |
log_lik[35] | -2.34996 | 0.000499418 | 0.129235 | -2.67219 | -2.58849 | -2.40913 | -2.32994 | -2.2649 | -2.18353 | -2.15836 | 66963 | 1.00008 |
log_lik[36] | -2.42536 | 0.000583436 | 0.184499 | -2.9031 | -2.7866 | -2.50546 | -2.38443 | -2.3008 | -2.20712 | -2.18051 | 100000 | 1.00001 |
log_lik[37] | -2.77796 | 0.00108888 | 0.344334 | -3.61273 | -3.42961 | -2.96609 | -2.71626 | -2.52376 | -2.33709 | -2.29548 | 100000 | 0.99997 |
log_lik[38] | -2.56802 | 0.000821847 | 0.259891 | -3.21964 | -3.07226 | -2.69821 | -2.51079 | -2.37974 | -2.25574 | -2.22447 | 100000 | 0.999978 |
log_lik[39] | -2.40097 | 0.000587076 | 0.168242 | -2.83704 | -2.72648 | -2.47175 | -2.36656 | -2.28945 | -2.19963 | -2.17254 | 82126 | 1.00004 |
log_lik[40] | -3.7057 | 0.00191992 | 0.607131 | -5.07782 | -4.81287 | -4.07006 | -3.63623 | -3.26529 | -2.83716 | -2.72761 | 100000 | 0.999971 |
log_lik[41] | -2.56802 | 0.000821847 | 0.259891 | -3.21964 | -3.07226 | -2.69821 | -2.51079 | -2.37974 | -2.25574 | -2.22447 | 100000 | 0.999978 |
log_lik[42] | -2.64915 | 0.000989698 | 0.31297 | -3.42632 | -3.25123 | -2.80993 | -2.57989 | -2.41906 | -2.27508 | -2.24037 | 100000 | 0.999983 |
log_lik[43] | -2.64915 | 0.000989698 | 0.31297 | -3.42632 | -3.25123 | -2.80993 | -2.57989 | -2.41906 | -2.27508 | -2.24037 | 100000 | 0.999983 |
log_lik[44] | -2.48078 | 0.000720873 | 0.22796 | -3.07134 | -2.92787 | -2.58129 | -2.42439 | -2.32413 | -2.22039 | -2.19265 | 100000 | 0.999998 |
log_lik[45] | -2.47998 | 0.0007191 | 0.227399 | -3.06829 | -2.92811 | -2.58068 | -2.42439 | -2.32295 | -2.2199 | -2.19207 | 100000 | 0.999991 |
log_lik[46] | -2.37927 | 0.000572994 | 0.1578 | -2.78904 | -2.68228 | -2.44339 | -2.34819 | -2.27626 | -2.19009 | -2.16486 | 75843 | 1.00002 |
log_lik[47] | -3.1857 | 0.00158996 | 0.502788 | -4.36319 | -4.12359 | -3.47549 | -3.11292 | -2.81382 | -2.50661 | -2.43336 | 100000 | 0.999983 |
log_lik[48] | -2.34584 | 0.00051513 | 0.126888 | -2.65972 | -2.57614 | -2.40448 | -2.32682 | -2.26275 | -2.18111 | -2.15621 | 60675 | 1.00006 |
log_lik[49] | -2.41088 | 0.000514388 | 0.162664 | -2.82099 | -2.72195 | -2.48676 | -2.37974 | -2.29996 | -2.20902 | -2.18182 | 100000 | 0.999999 |
log_lik[50] | -2.76093 | 0.000962221 | 0.304281 | -3.48081 | -3.32983 | -2.932 | -2.71175 | -2.53689 | -2.35856 | -2.31605 | 100000 | 0.999983 |
log_lik[51] | -2.55227 | 0.000724273 | 0.229035 | -3.11222 | -2.98866 | -2.67227 | -2.5078 | -2.38647 | -2.26473 | -2.23413 | 100000 | 0.999982 |
log_lik[52] | -2.32991 | 0.000421662 | 0.109997 | -2.58737 | -2.52469 | -2.38765 | -2.31737 | -2.25642 | -2.17761 | -2.15289 | 68051 | 1.00014 |
log_lik[53] | -2.97519 | 0.001168 | 0.369353 | -3.83217 | -3.65693 | -3.18955 | -2.92343 | -2.70504 | -2.47099 | -2.41237 | 100000 | 1.00001 |
log_lik[54] | -2.32991 | 0.000421662 | 0.109997 | -2.58737 | -2.52469 | -2.38765 | -2.31737 | -2.25642 | -2.17761 | -2.15289 | 68051 | 1.00014 |
log_lik[55] | -2.41088 | 0.000514388 | 0.162664 | -2.82099 | -2.72195 | -2.48676 | -2.37974 | -2.29996 | -2.20902 | -2.18182 | 100000 | 0.999999 |
log_lik[56] | -2.76093 | 0.000962221 | 0.304281 | -3.48081 | -3.32983 | -2.932 | -2.71175 | -2.53689 | -2.35856 | -2.31605 | 100000 | 0.999983 |
log_lik[57] | -2.55227 | 0.000724273 | 0.229035 | -3.11222 | -2.98866 | -2.67227 | -2.5078 | -2.38647 | -2.26473 | -2.23413 | 100000 | 0.999982 |
log_lik[58] | -5.96118 | 0.00310773 | 0.982751 | -8.08284 | -7.69382 | -6.58038 | -5.88692 | -5.26594 | -4.47592 | -4.2481 | 100000 | 1 |
log_lik[59] | -2.52657 | 0.000757293 | 0.239477 | -3.12965 | -2.99042 | -2.64358 | -2.47232 | -2.35648 | -2.24149 | -2.21098 | 100000 | 0.999994 |
log_lik[60] | -2.42542 | 0.000581973 | 0.184036 | -2.89803 | -2.78347 | -2.50676 | -2.38427 | -2.30052 | -2.2072 | -2.17959 | 100000 | 1.00003 |
log_lik[61] | -3.70403 | 0.00191809 | 0.606552 | -5.07289 | -4.8054 | -4.06898 | -3.63492 | -3.26146 | -2.84054 | -2.73181 | 100000 | 1 |
log_lik[62] | -2.34145 | 0.000491539 | 0.121126 | -2.63625 | -2.55886 | -2.40019 | -2.32462 | -2.26172 | -2.18046 | -2.15654 | 60724 | 1.00006 |
log_lik[63] | -2.3498 | 0.000495617 | 0.128819 | -2.66905 | -2.58578 | -2.40933 | -2.33002 | -2.26458 | -2.18315 | -2.158 | 67557 | 1.00007 |
log_lik[64] | -2.34145 | 0.000491539 | 0.121126 | -2.63625 | -2.55886 | -2.40019 | -2.32462 | -2.26172 | -2.18046 | -2.15654 | 60724 | 1.00006 |
log_lik[65] | -3.05587 | 0.00137053 | 0.4334 | -4.06844 | -3.86188 | -3.3063 | -2.9924 | -2.73696 | -2.46901 | -2.40783 | 100000 | 0.999993 |
log_lik[66] | -2.5683 | 0.000819629 | 0.259189 | -3.21226 | -3.06855 | -2.70014 | -2.51147 | -2.38054 | -2.25607 | -2.22414 | 100000 | 1.00001 |
log_lik[67] | -2.49983 | 0.00075334 | 0.238227 | -3.10952 | -2.96599 | -2.60912 | -2.4412 | -2.33334 | -2.22697 | -2.19767 | 100000 | 1 |
log_lik[68] | -2.36999 | 0.00057377 | 0.149574 | -2.75292 | -2.65382 | -2.43233 | -2.34244 | -2.27251 | -2.18818 | -2.16267 | 67957 | 1.00001 |
log_lik[69] | -2.36999 | 0.00057377 | 0.149574 | -2.75292 | -2.65382 | -2.43233 | -2.34244 | -2.27251 | -2.18818 | -2.16267 | 67957 | 1.00001 |
log_lik[70] | -2.36999 | 0.00057377 | 0.149574 | -2.75292 | -2.65382 | -2.43233 | -2.34244 | -2.27251 | -2.18818 | -2.16267 | 67957 | 1.00001 |
log_lik[71] | -2.46125 | 0.000725389 | 0.215494 | -3.01841 | -2.88701 | -2.5549 | -2.4086 | -2.31443 | -2.21489 | -2.1868 | 88253 | 0.999968 |
log_lik[72] | -2.38928 | 0.000607463 | 0.165664 | -2.81825 | -2.70878 | -2.45625 | -2.35539 | -2.28072 | -2.19397 | -2.16786 | 74373 | 1.00004 |
log_lik[73] | -2.38928 | 0.000607463 | 0.165664 | -2.81825 | -2.70878 | -2.45625 | -2.35539 | -2.28072 | -2.19397 | -2.16786 | 74373 | 1.00004 |
log_lik[74] | -3.67851 | 0.0018051 | 0.570824 | -4.96139 | -4.70864 | -4.02303 | -3.61492 | -3.26414 | -2.86209 | -2.75324 | 100000 | 0.999999 |
log_lik[75] | -2.42412 | 0.000559944 | 0.17707 | -2.87225 | -2.76381 | -2.50505 | -2.38741 | -2.30404 | -2.20929 | -2.18255 | 100000 | 1.00006 |
log_lik[76] | -4.85816 | 0.00258955 | 0.818888 | -6.64133 | -6.31093 | -5.36921 | -4.78662 | -4.26925 | -3.64582 | -3.46943 | 100000 | 0.999995 |
log_lik[77] | -2.39 | 0.000493201 | 0.155964 | -2.78444 | -2.6866 | -2.4595 | -2.36016 | -2.28576 | -2.1976 | -2.17153 | 100000 | 1 |
log_lik[78] | -2.9613 | 0.00121791 | 0.385138 | -3.86024 | -3.67563 | -3.1825 | -2.90355 | -2.6775 | -2.44311 | -2.38634 | 100000 | 0.999986 |
log_lik[79] | -2.39 | 0.000493201 | 0.155964 | -2.78444 | -2.6866 | -2.4595 | -2.36016 | -2.28576 | -2.1976 | -2.17153 | 100000 | 1 |
log_lik[80] | -3.41139 | 0.00160314 | 0.506957 | -4.55559 | -4.33322 | -3.71439 | -3.35026 | -3.03988 | -2.69708 | -2.60639 | 100000 | 0.999997 |
log_lik[81] | -3.67851 | 0.0018051 | 0.570824 | -4.96139 | -4.70864 | -4.02303 | -3.61492 | -3.26414 | -2.86209 | -2.75324 | 100000 | 0.999999 |
log_lik[82] | -2.78322 | 0.00103984 | 0.328827 | -3.56402 | -3.39689 | -2.96697 | -2.72775 | -2.54094 | -2.35577 | -2.31086 | 100000 | 1.00001 |
log_lik[83] | -2.84325 | 0.0011636 | 0.367961 | -3.72012 | -3.53811 | -3.04991 | -2.78014 | -2.57055 | -2.36481 | -2.31905 | 100000 | 1 |
log_lik[84] | -2.36323 | 0.000525796 | 0.140284 | -2.71569 | -2.62599 | -2.42459 | -2.33932 | -2.27104 | -2.18727 | -2.16231 | 71184 | 0.999996 |
log_lik[85] | -2.66497 | 0.000953673 | 0.301578 | -3.40396 | -3.24526 | -2.82653 | -2.60503 | -2.44195 | -2.29206 | -2.25465 | 100000 | 1.00001 |
log_lik[86] | -2.37956 | 0.000542134 | 0.152967 | -2.77117 | -2.67185 | -2.44493 | -2.35106 | -2.27869 | -2.19225 | -2.16674 | 79613 | 1 |
log_lik[87] | -2.37956 | 0.000542134 | 0.152967 | -2.77117 | -2.67185 | -2.44493 | -2.35106 | -2.27869 | -2.19225 | -2.16674 | 79613 | 1 |
log_lik[88] | -2.37956 | 0.000542134 | 0.152967 | -2.77117 | -2.67185 | -2.44493 | -2.35106 | -2.27869 | -2.19225 | -2.16674 | 79613 | 1 |
log_lik[89] | -2.61597 | 0.000890215 | 0.281511 | -3.30934 | -3.15762 | -2.7638 | -2.55755 | -2.40974 | -2.27316 | -2.23913 | 100000 | 1 |
log_lik[90] | -2.37956 | 0.000542134 | 0.152967 | -2.77117 | -2.67185 | -2.44493 | -2.35106 | -2.27869 | -2.19225 | -2.16674 | 79613 | 1 |
log_lik[91] | -2.35169 | 0.000524891 | 0.13297 | -2.6827 | -2.595 | -2.41106 | -2.33035 | -2.26492 | -2.18286 | -2.1578 | 64175 | 0.999993 |
log_lik[92] | -3.39297 | 0.00177845 | 0.562394 | -4.69445 | -4.4338 | -3.71977 | -3.31813 | -2.98064 | -2.61791 | -2.52666 | 100000 | 0.999975 |
log_lik[93] | -2.56962 | 0.000873156 | 0.276116 | -3.26589 | -3.1083 | -2.7045 | -2.50481 | -2.37047 | -2.24856 | -2.21752 | 100000 | 0.999975 |
log_lik[94] | -2.73918 | 0.00110666 | 0.349957 | -3.59245 | -3.40636 | -2.92461 | -2.66867 | -2.48011 | -2.3068 | -2.26847 | 100000 | 0.999968 |
log_lik[95] | -2.35169 | 0.000524891 | 0.13297 | -2.6827 | -2.595 | -2.41106 | -2.33035 | -2.26492 | -2.18286 | -2.1578 | 64175 | 0.999993 |
log_lik[96] | -3.73209 | 0.00205986 | 0.651384 | -5.20444 | -4.91743 | -4.1244 | -3.65239 | -3.25571 | -2.80926 | -2.69637 | 100000 | 0.999975 |
log_lik[97] | -2.56962 | 0.000873156 | 0.276116 | -3.26589 | -3.1083 | -2.7045 | -2.50481 | -2.37047 | -2.24856 | -2.21752 | 100000 | 0.999975 |
水準の差
changeup | curve | cutball | fork | slider | straight | |
---|---|---|---|---|---|---|
changeup | 0 | 0.99474 | 0 | 1e-05 | 0.37879 | 0 |
curve | 0.00526 | 0 | 0 | 0 | 0.00256 | 0 |
cutball | 1 | 1 | 0 | 1 | 1 | 0 |
fork | 0.99999 | 1 | 0 | 0 | 1 | 0 |
slider | 0.62121 | 0.99744 | 0 | 0 | 0 | 0 |
straight | 1 | 1 | 1 | 1 | 1 | 0 |
考察
- 水準の効果の分布より、走者の状況も交絡因子も殆ど影響力が無いことが推し量れる
- 説明量や効果量についても、圧倒的に走者(b)についての値が大きく、モデルの殆どがこれで説明される事がわかる
- 交絡因子の大きさについても図れるのはメリットぽい。
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その5
その5です。今回は第4章の章末問題に取り組んでいきます。
4章 対応ある2群の差と相関の推測
内容
- 対応ある2群のt検定のオルタナとして機能します
- 対応ある、の意味とは?
- 同じ観察対象から2回測定しているもの
- beforeの体重とafterの体重のセット * n個など
- この解析をするときには実験デザインが大事になります
- どちらかの群にバイアスがかかることを避ける事が必要です
- 対応ある2群の実験デザインを行う
- マッチング : 施策実行前の状態が同じ2つを組にして、ランダムに2群に割り当てる
- プリテスト・ポストテスト : 施策の前後で同じ対象を観察する
- 相関関係を表現する要約統計量の例
- 共分散
- 平均偏差の積の平均値
- 正の相関がある時に正、負の相関がある時に負の値になる
- 相関の強さは表現できない
- 相関係数 or
- データの正規化を行った後に、積の平均値を計算
- -1 <= r <= 1
- 共分散
- 2変量正規分布の導入
- 共分散: の関係性を持つ
- 正規分布に従う2変量が観測される確率をモデル化する際に、あてはめが可能となる理論分布
- 2群の差異の考察のバリエーション
- 独立した2群の差の分析
- 対応ある2群の群間差の分析(Inter)
- 対応ある2群の個人内差の分析(Intra)
- 対応ある場合の生成量(2変量に相関が無い場合は、とすればokです
章末問題
データの取得
取り敢えずデータを準備します。
import numpy as np import os from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger logger = get_logger() # データの準備 a = np.array([62,54,19,54,47,22,35,77,64,60,27,41,41,44,57,16,42,89,40,67,69,46,74,62,60,87,32,42,73,25,42,57,31,35,33,38,43,53,55,62,67,56,76,5,31,70,66,65,34,48]) b = np.array([73,72,56,58,71,42,78,77,75,72,56,71,69,77,84,51,62,88,56,58,84,91,71,82,81,77,65,78,79,60,66,70,65,57,64,61,56,67,75,64,68,67,80,55,48,85,56,62,65,79]) x = np.stack((a, b), axis=1) # 出力用ディレクトリの用意 result_dir = os.path.join("result", "chapter4") if os.path.exists(result_dir) is False: os.makedirs(result_dir) logger.info("{0} is made".format(result_dir))
各種統計量の計算と可視化
モデルから推定をする前に、統計量を計算していきます。 対応関係を使って散布図を書けば相関が見えてきますね。
import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from tabulate import tabulate from IPython.core.display import display # データの要約統計量を計算 df = pd.DataFrame({"a":a, "b":b}) df_desc = df.describe() display(df_desc) df_desc_path = os.path.join(result_dir, "df_describe.md") with open(df_desc_path, "w") as f: f.write(tabulate(df_desc, df_desc.columns, tablefmt="pipe")) logger.info("sample data summary is saved at {0}".format(df_desc_path)) # 散布図を可視化 sns.jointplot(data=df, x="a", y="b") jointplot_path = os.path.join(result_dir, "jointplot.png") plt.show() plt.savefig(jointplot_path) logger.info("jointplot result is saved at {0}".format(jointplot_path)) r = np.corrcoef(x.T)[0][1] logger.info("ピアソン相関係数は{0}です".format(r)) s = np.cov(x.T)[0][1] logger.info("共分散は{0}です".format(s))
- 結果(表と分布)
a | b | |
---|---|---|
count | 50 | 50 |
mean | 49.9 | 68.48 |
std | 18.6703 | 10.9958 |
min | 5 | 42 |
25% | 35.75 | 60.25 |
50% | 50.5 | 68.5 |
75% | 63.5 | 77 |
max | 89 | 91 |
Stanによる統計モデルの構築
これもしょうもないテクニックなんですが、2次元配列 or 行列がパラメータとして得られている場合には1次元の配列に入れています。PyStanの可視化部分が多次元構造に対応しておらず、突然カーネルが死んでしまうためです(summaryメソッドは対応しています)。PyStanの可視化は昔のPyMC3のモジュールを利用しているのですが、完璧に対応できていないことが原因のようです。その内PyMC3のチームが作っているらしいmcmcplotlibというモジュールに移行する予定らしいですが、まだその雰囲気はありません…
また同順率を計算する時にsinの逆関数としてasinを利用しました。
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "g2_pair.stan") stan_file_c = os.path.join("stan", "g2_pair.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f) logger.info("Stan model is compiled to {0}".format(stan_file_c))
- Stan
data { int<lower=0> n ; vector[2] x[n] ; real c_mu_diff ; real c_es ; real c_cohenu ; real c_pod ; real c_pbt ; real<lower=0, upper=1> cdash_pbt ; real c_diff_sd ; real c_pair_es ; real<lower=0, upper=1> c_pair_pod ; real<lower=0, upper=1> cdash_pair_pbt ; real c_rho ; real<lower=0, upper=1> c_poc; } parameters { vector[2] mu ; real<lower=0> sigma_a ; real<lower=0> sigma_b ; real<lower=-1, upper=1> rho ; } transformed parameters { # ダイレクトに共分散行列を与えると # PyStanの可視化でエラー無しに落ちるので、Stan内部で作る real<lower=0> sigma_ab ; cov_matrix[2] Sigma ; sigma_ab = sigma_a * sigma_b * rho ; Sigma[1,1] = pow(sigma_a, 2) ; Sigma[1,2] = sigma_ab ; Sigma[2,1] = sigma_ab ; Sigma[2,2] = pow(sigma_b, 2) ; } model { for(i in 1:n){ x[i] ~ multi_normal(mu, Sigma) ; } } generated quantities { vector[n] log_lik ; real mu_diff ; real es_a ; real es_b ; real cohenu_a ; real cohenu_b ; real<lower=0, upper=1> pod ; real<lower=0, upper=1> pbt ; real diff_sd ; real pair_es ; real pair_pod ; real pair_pbt ; real<lower=0, upper=1> poc ; int<lower=0, upper=1> prob_mu_diff_upper_0 ; int<lower=0, upper=1> prob_mu_diff_upper_c ; int<lower=0, upper=1> prob_es_a_upper_c ; int<lower=0, upper=1> prob_es_b_upper_c ; int<lower=0, upper=1> prob_cohenu_a_upper_c ; int<lower=0, upper=1> prob_pod_upper_c ; int<lower=0, upper=1> prob_pbt_upper_cdash ; int<lower=0, upper=1> prob_diff_sd_upper_c ; int<lower=0, upper=1> prob_pair_es_upper_c ; int<lower=0, upper=1> prob_pair_pod_upper_c ; int<lower=0, upper=1> prob_pair_pbt_upper_cdash ; int<lower=0, upper=1> prob_rho_upper_c ; int<lower=0, upper=1> prob_poc_upper_c ; for(i in 1:n){ log_lik[i] = multi_normal_lpdf(x[i] | mu, Sigma) ; } mu_diff = mu[1] - mu[2] ; es_a = fabs(mu[1] - mu[2]) / sigma_a ; es_b = fabs(mu[1] - mu[2]) / sigma_b ; cohenu_a = normal_cdf(mu[2], mu[1], sigma_a) ; cohenu_b = normal_cdf(mu[1], mu[2], sigma_b) ; pod = normal_cdf(mu_diff / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; pbt = normal_cdf((mu_diff - c_pbt) / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; diff_sd = sqrt(pow(sigma_a, 2) + pow(sigma_b, 2) - 2 * rho * sigma_a * sigma_b) ; pair_es = fabs(mu_diff) / diff_sd; pair_pod = normal_cdf(pair_es, 0, 1) ; pair_pbt = normal_cdf((mu_diff - c_pbt) / diff_sd, 0, 1) ; poc = 0.5 + 1 / pi() * asin(rho) ; prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ; prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ; prob_es_a_upper_c = es_a > c_es ? 1 : 0 ; prob_es_b_upper_c = es_b > c_es ? 1 : 0 ; prob_cohenu_a_upper_c = cohenu_a > c_cohenu ? 1 : 0 ; prob_pod_upper_c = pod > c_pod ? 1 : 0 ; prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ; prob_diff_sd_upper_c = diff_sd > c_diff_sd ? 1 : 0 ; prob_pair_es_upper_c = pair_es > c_pair_es ? 1 : 0 ; prob_pair_pod_upper_c = pair_pod > c_pair_pod ? 1 : 0 ; prob_pair_pbt_upper_cdash = pair_pbt > cdash_pair_pbt ? 1 : 0 ; prob_rho_upper_c = rho > c_rho ? 1 : 0 ; prob_poc_upper_c = poc > c_poc ? 1 : 0 ; }
統計モデルによる事後分布のサンプリング
prob_mu_diff_upper_0については全てのサンプリング値で0となり、pystanのモジュールで可視化出来ない(カーネル密度推定出来ない)のでその直前で除きます。
import pandas as pd import pickle import pystan import matplotlib import matplotlib.pyplot as plt import os from IPython.core.display import display from tabulate import tabulate matplotlib.rcParams['figure.figsize'] = (10, 80) # Stanで使うデータの用意 stan_data = {"n": len(x), "x": x, "c_mu_diff": -15, "c_es": 1.0, "c_cohenu": 0.8, "c_pod": 0.1, "c_pbt": 15, "cdash_pbt": 0.05, "c_diff_sd": 20, "c_pair_es": 1.0, "c_pair_pod": 0.8, "cdash_pair_pbt": 0.05, "c_rho": 0.75, "c_poc": 0.75} # 興味のあるパラメータの設定 par = ["mu", "sigma_a", "sigma_b", "rho", "sigma_ab", "mu_diff", "es_a", "es_b", "cohenu_a", "cohenu_b", "pod", "pbt", "diff_sd", "pair_es", "pair_pod", "pair_pbt", "poc", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_a_upper_c", "prob_es_b_upper_c", "prob_cohenu_a_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash", "prob_diff_sd_upper_c", "prob_pair_es_upper_c", "prob_pair_pod_upper_c", "prob_pair_pbt_upper_cdash", "prob_rho_upper_c", "prob_poc_upper_c", "log_lik"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_pair.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング logger.info("Start MCMC sampling") fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) summary_df_path = os.path.join(result_dir, "df_summary.md") with open(summary_df_path, "w") as f: f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe")) logger.info("MCMC result summary is saved at {0}".format(summary_df_path)) # 事後分布の可視化 par.remove("prob_mu_diff_upper_0") fit.traceplot(par) traceplot_path = os.path.join(result_dir, "traceplot.png") plt.savefig(traceplot_path) plt.show() logger.info("traceplot result is saved at {0}".format(traceplot_path)) # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu[0] | 49.8926 | 0.00957801 | 2.73487 | 44.4876 | 45.3974 | 48.0818 | 49.8909 | 51.7234 | 54.3679 | 55.2544 | 81531 | 1.00003 |
mu[1] | 68.4765 | 0.00566431 | 1.60689 | 65.3396 | 65.8384 | 67.4029 | 68.4804 | 69.5427 | 71.1116 | 71.6325 | 80478 | 1.00001 |
sigma_a | 19.2229 | 0.00716584 | 2.00791 | 15.7698 | 16.2345 | 17.8062 | 19.061 | 20.45 | 22.7575 | 23.6287 | 78515 | 0.999995 |
sigma_b | 11.3205 | 0.00418226 | 1.18074 | 9.2944 | 9.56625 | 10.4895 | 11.2246 | 12.0418 | 13.4129 | 13.9161 | 79705 | 0.999964 |
rho | 0.594667 | 0.000335941 | 0.0935382 | 0.390446 | 0.428657 | 0.536483 | 0.602234 | 0.661714 | 0.733422 | 0.754116 | 77527 | 0.999996 |
sigma_ab | 131.615 | 0.153909 | 38.1313 | 69.8853 | 77.6958 | 104.862 | 127.231 | 153.284 | 200.727 | 219.867 | 61381 | 0.999985 |
mu_diff | -18.5839 | 0.00693064 | 2.19166 | -22.9048 | -22.1834 | -20.0325 | -18.5904 | -17.1309 | -14.9834 | -14.2762 | 100000 | 1.00003 |
es_a | 0.977031 | 0.000478766 | 0.151399 | 0.687624 | 0.732986 | 0.873885 | 0.97449 | 1.07755 | 1.22976 | 1.27962 | 100000 | 1.00005 |
es_b | 1.65908 | 0.000820216 | 0.259375 | 1.18273 | 1.25308 | 1.47896 | 1.64867 | 1.82674 | 2.10267 | 2.19663 | 100000 | 1 |
cohenu_a | 0.832983 | 0.000118209 | 0.0373811 | 0.754155 | 0.768216 | 0.808909 | 0.835093 | 0.859383 | 0.890606 | 0.89966 | 100000 | 1.00004 |
cohenu_b | 0.0540274 | 8.64856e-05 | 0.0273491 | 0.0140236 | 0.0177475 | 0.0338697 | 0.0496074 | 0.0695757 | 0.105088 | 0.118457 | 100000 | 0.999984 |
pod | 0.202597 | 0.000108955 | 0.0344545 | 0.139349 | 0.148333 | 0.178551 | 0.201236 | 0.225017 | 0.261581 | 0.273534 | 100000 | 1.00004 |
pbt | 0.0674462 | 7.64904e-05 | 0.0216195 | 0.0324478 | 0.0364988 | 0.0518413 | 0.065037 | 0.0803561 | 0.10656 | 0.116507 | 79887 | 1.00001 |
diff_sd | 15.4018 | 0.00515225 | 1.62928 | 12.6101 | 12.9776 | 14.2577 | 15.2707 | 16.3993 | 18.2913 | 18.9678 | 100000 | 0.999976 |
pair_es | 1.21977 | 0.000598837 | 0.189369 | 0.851131 | 0.910184 | 1.09139 | 1.21844 | 1.34708 | 1.5329 | 1.59392 | 100000 | 1.00002 |
pair_pod | 0.884636 | 0.00011549 | 0.0365212 | 0.802652 | 0.818637 | 0.86245 | 0.888472 | 0.911022 | 0.93735 | 0.944523 | 100000 | 1.00001 |
pair_pbt | 0.0165938 | 3.57684e-05 | 0.0113109 | 0.00309664 | 0.00399766 | 0.0085756 | 0.013873 | 0.0216362 | 0.0384226 | 0.0455154 | 100000 | 0.999988 |
poc | 0.704289 | 0.000133159 | 0.0369491 | 0.627679 | 0.641013 | 0.680247 | 0.705723 | 0.730171 | 0.762078 | 0.771934 | 76996 | 0.999993 |
prob_mu_diff_upper_0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 100000 | nan |
prob_mu_diff_upper_c | 0.05073 | 0.000763322 | 0.219447 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 82650 | 0.999997 |
prob_es_a_upper_c | 0.43221 | 0.00156655 | 0.495386 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 100000 | 1.00004 |
prob_es_b_upper_c | 0.99731 | 0.000176451 | 0.0517957 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 86167 | 1.00002 |
prob_cohenu_a_upper_c | 0.81293 | 0.00136484 | 0.38997 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 81639 | 0.999982 |
prob_pod_upper_c | 0.99972 | 5.65467e-05 | 0.0167309 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 87544 | 1.00006 |
prob_pbt_upper_cdash | 0.78259 | 0.0013449 | 0.412486 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 94067 | 1.00003 |
prob_diff_sd_upper_c | 0.00844 | 0.00033347 | 0.0914814 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 75258 | 1.00001 |
prob_pair_es_upper_c | 0.87754 | 0.00116844 | 0.327818 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 78714 | 0.999991 |
prob_pair_pod_upper_c | 0.97798 | 0.000521882 | 0.146749 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 79069 | 0.999996 |
prob_pair_pbt_upper_cdash | 0.01667 | 0.000463955 | 0.128032 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 76153 | 1.00003 |
prob_rho_upper_c | 0.0288 | 0.000634729 | 0.167245 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 69427 | 1.00002 |
prob_poc_upper_c | 0.10329 | 0.0011511 | 0.304339 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 69902 | 0.999987 |
log_lik[0] | -7.20677 | 0.000510806 | 0.153203 | -7.51933 | -7.46593 | -7.30743 | -7.20317 | -7.10048 | -6.96204 | -6.9202 | 89954 | 1.00002 |
log_lik[1] | -7.04915 | 0.000502755 | 0.148317 | -7.35091 | -7.29856 | -7.14718 | -7.04574 | -6.94693 | -6.81087 | -6.76865 | 87030 | 0.999998 |
log_lik[2] | -8.37423 | 0.0010673 | 0.313067 | -9.05396 | -8.9247 | -8.56952 | -8.3508 | -8.15363 | -7.90214 | -7.82931 | 86041 | 0.999984 |
log_lik[3] | -7.93282 | 0.000748297 | 0.236632 | -8.43927 | -8.34739 | -8.0812 | -7.91782 | -7.76714 | -7.57149 | -7.51186 | 100000 | 1.00001 |
log_lik[4] | -7.09005 | 0.000503034 | 0.148439 | -7.39144 | -7.3404 | -7.18795 | -7.08658 | -6.98816 | -6.85275 | -6.80871 | 87077 | 1 |
log_lik[5] | -9.88101 | 0.00200767 | 0.594664 | -11.1672 | -10.9272 | -10.2558 | -9.83603 | -9.45866 | -8.98453 | -8.85144 | 87732 | 0.999988 |
log_lik[6] | -8.71399 | 0.00119432 | 0.377678 | -9.53705 | -9.38041 | -8.9495 | -8.68412 | -8.44366 | -8.15071 | -8.06613 | 100000 | 1.00001 |
log_lik[7] | -8.05166 | 0.000808855 | 0.255782 | -8.60068 | -8.49872 | -8.2127 | -8.03498 | -7.87169 | -7.66212 | -7.60015 | 100000 | 1.00006 |
log_lik[8] | -7.29648 | 0.000531128 | 0.158762 | -7.62103 | -7.56603 | -7.40109 | -7.29165 | -7.18636 | -7.04438 | -7.00056 | 89350 | 1.00003 |
log_lik[9] | -7.14304 | 0.000503928 | 0.150479 | -7.44972 | -7.39716 | -7.24229 | -7.13974 | -7.03927 | -6.90195 | -6.85998 | 89169 | 1.00001 |
log_lik[10] | -7.86887 | 0.000793813 | 0.227803 | -8.35733 | -8.26465 | -8.01223 | -7.85569 | -7.70969 | -7.51842 | -7.4618 | 82354 | 0.999979 |
log_lik[11] | -7.31585 | 0.000503531 | 0.15923 | -7.6417 | -7.58604 | -7.41964 | -7.31028 | -7.20657 | -7.06351 | -7.01762 | 100000 | 1.00002 |
log_lik[12] | -7.19628 | 0.000480615 | 0.151984 | -7.50551 | -7.45265 | -7.29581 | -7.19238 | -7.09229 | -6.95323 | -6.90778 | 100000 | 1.00001 |
log_lik[13] | -7.76793 | 0.000671836 | 0.212453 | -8.21805 | -8.13685 | -7.90243 | -7.75448 | -7.61944 | -7.44257 | -7.38864 | 100000 | 1.00003 |
log_lik[14] | -8.15601 | 0.000870612 | 0.275312 | -8.75044 | -8.63974 | -8.32869 | -8.137 | -7.96147 | -7.7388 | -7.67263 | 100000 | 1 |
log_lik[15] | -8.83067 | 0.0014309 | 0.395669 | -9.69151 | -9.52681 | -9.07889 | -8.79907 | -8.55174 | -8.23532 | -8.14442 | 76462 | 0.999986 |
log_lik[16] | -7.17365 | 0.000504963 | 0.151681 | -7.4807 | -7.43007 | -7.27362 | -7.1703 | -7.06917 | -6.93037 | -6.88726 | 90228 | 0.999968 |
log_lik[17] | -9.38805 | 0.00180669 | 0.499185 | -10.4735 | -10.2703 | -9.70282 | -9.35085 | -9.03296 | -8.63428 | -8.52274 | 76341 | 1.00004 |
log_lik[18] | -7.6546 | 0.000618485 | 0.195582 | -8.06579 | -7.99312 | -7.77964 | -7.64522 | -7.51888 | -7.34957 | -7.29857 | 100000 | 0.999987 |
log_lik[19] | -9.16164 | 0.00145992 | 0.461667 | -10.1567 | -9.98055 | -9.45124 | -9.12617 | -8.83245 | -8.47093 | -8.36446 | 100000 | 1 |
log_lik[20] | -8.01269 | 0.00084993 | 0.250002 | -8.54689 | -8.45199 | -8.17131 | -7.99616 | -7.83703 | -7.6313 | -7.57052 | 86521 | 1.00001 |
log_lik[21] | -10.6807 | 0.00242399 | 0.766533 | -12.3495 | -12.0334 | -11.1564 | -10.6244 | -10.1335 | -9.53191 | -9.35498 | 100000 | 0.999991 |
log_lik[22] | -8.05483 | 0.000812016 | 0.256782 | -8.609 | -8.50441 | -8.21509 | -8.03812 | -7.87417 | -7.6654 | -7.60062 | 100000 | 1.00004 |
log_lik[23] | -7.7555 | 0.000666349 | 0.210718 | -8.20175 | -8.11848 | -7.89039 | -7.74328 | -7.6084 | -7.43055 | -7.37751 | 100000 | 1.00001 |
log_lik[24] | -7.65779 | 0.00062335 | 0.19712 | -8.07091 | -7.99542 | -7.78414 | -7.64716 | -7.52059 | -7.35168 | -7.30084 | 100000 | 1.00001 |
log_lik[25] | -9.09361 | 0.00140706 | 0.444952 | -10.0591 | -9.8797 | -9.37279 | -9.05784 | -8.77697 | -8.42635 | -8.32266 | 100000 | 1.00004 |
log_lik[26] | -7.50563 | 0.000563826 | 0.178297 | -7.87608 | -7.81074 | -7.6211 | -7.49828 | -7.38201 | -7.22696 | -7.17773 | 100000 | 0.999996 |
log_lik[27] | -8.05778 | 0.000818254 | 0.258755 | -8.61498 | -8.51212 | -8.22005 | -8.03921 | -7.87442 | -7.66727 | -7.60576 | 100000 | 1.00003 |
log_lik[28] | -7.79465 | 0.000734035 | 0.215728 | -8.25129 | -8.16903 | -7.93286 | -7.78268 | -7.64341 | -7.46225 | -7.40685 | 86373 | 1.00006 |
log_lik[29] | -7.88166 | 0.000749915 | 0.230262 | -8.37495 | -8.28305 | -8.02686 | -7.86715 | -7.72116 | -7.52979 | -7.47299 | 94280 | 0.999984 |
log_lik[30] | -7.08697 | 0.00050002 | 0.148871 | -7.38986 | -7.33785 | -7.18472 | -7.08402 | -6.98546 | -6.84782 | -6.80542 | 88643 | 0.999981 |
log_lik[31] | -7.0756 | 0.000500711 | 0.148709 | -7.37755 | -7.32674 | -7.17359 | -7.072 | -6.97293 | -6.83728 | -6.79434 | 88206 | 0.999993 |
log_lik[32] | -7.57171 | 0.000588542 | 0.186113 | -7.9608 | -7.89077 | -7.69187 | -7.56353 | -7.44214 | -7.28218 | -7.23107 | 100000 | 0.999997 |
log_lik[33] | -7.56191 | 0.000618204 | 0.184662 | -7.94683 | -7.87941 | -7.68109 | -7.55408 | -7.43403 | -7.27271 | -7.22362 | 89226 | 0.999976 |
log_lik[34] | -7.41856 | 0.000535287 | 0.169273 | -7.76974 | -7.70685 | -7.52824 | -7.41221 | -7.30165 | -7.15134 | -7.10458 | 100000 | 0.999989 |
log_lik[35] | -7.26815 | 0.000523735 | 0.157009 | -7.58864 | -7.53398 | -7.37085 | -7.26407 | -7.15992 | -7.01727 | -6.97304 | 89873 | 0.999969 |
log_lik[36] | -7.71186 | 0.00064242 | 0.203151 | -8.14296 | -8.06432 | -7.84074 | -7.70171 | -7.57035 | -7.39709 | -7.34306 | 100000 | 0.999994 |
log_lik[37] | -7.05357 | 0.000500421 | 0.148306 | -7.35433 | -7.30301 | -7.15194 | -7.05027 | -6.95145 | -6.81583 | -6.77401 | 87831 | 0.999979 |
log_lik[38] | -7.17771 | 0.000505919 | 0.151559 | -7.48541 | -7.43305 | -7.278 | -7.17341 | -7.07353 | -6.93566 | -6.89261 | 89743 | 1.00001 |
log_lik[39] | -7.69658 | 0.000637628 | 0.201636 | -8.12377 | -8.04458 | -7.82498 | -7.68571 | -7.55632 | -7.38327 | -7.33149 | 100000 | 1.00001 |
log_lik[40] | -7.6861 | 0.000633419 | 0.200305 | -8.11 | -8.03307 | -7.81326 | -7.67564 | -7.54736 | -7.37526 | -7.32353 | 100000 | 1.00002 |
log_lik[41] | -7.135 | 0.000501372 | 0.14986 | -7.43817 | -7.3869 | -7.2345 | -7.13109 | -7.03222 | -6.89567 | -6.85188 | 89341 | 0.999986 |
log_lik[42] | -8.00336 | 0.000847046 | 0.247836 | -8.53242 | -8.435 | -8.16072 | -7.98708 | -7.82881 | -7.6244 | -7.56385 | 85608 | 1.00006 |
log_lik[43] | -9.90667 | 0.00191457 | 0.605442 | -11.2268 | -10.9773 | -10.2848 | -9.8613 | -9.47535 | -8.99769 | -8.86278 | 100000 | 0.999986 |
log_lik[44] | -8.72917 | 0.00118485 | 0.374681 | -9.5436 | -9.38778 | -8.96398 | -8.70135 | -8.46409 | -8.16497 | -8.07738 | 100000 | 0.999993 |
log_lik[45] | -8.14419 | 0.000925485 | 0.271828 | -8.72857 | -8.62312 | -8.31587 | -8.12515 | -7.95293 | -7.73213 | -7.66608 | 86268 | 1 |
log_lik[46] | -9.48036 | 0.00165699 | 0.523985 | -10.6101 | -10.409 | -9.80909 | -9.43988 | -9.10587 | -8.69649 | -8.57541 | 100000 | 1 |
log_lik[47] | -8.21578 | 0.000900439 | 0.284744 | -8.82678 | -8.71809 | -8.39352 | -8.19606 | -8.01586 | -7.78551 | -7.7145 | 100000 | 1.00001 |
log_lik[48] | -7.38682 | 0.000525221 | 0.16609 | -7.7291 | -7.66924 | -7.49486 | -7.38077 | -7.27242 | -7.12421 | -7.07814 | 100000 | 0.999993 |
log_lik[49] | -7.80583 | 0.000690243 | 0.218274 | -8.26971 | -8.18366 | -7.94377 | -7.79197 | -7.65319 | -7.47189 | -7.41726 | 100000 | 1.00003 |
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その4.1
対応のあるt検定を行う前に、3章の内容を等分散モデルで書くとどうなるか、をちゃんと検証します。
番外編 3章の等分散モデルと変分ベイズによる推定
等分散モデルによる差の推定
といっても、Stanのモデルをちょっと変えるだけなのでした。
まずコンパイル
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする # 等分散モデル stan_file = os.path.join("stan", "g2_ind_equ.stan") stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f)
Stanのモデルはこんな感じです。
data { int<lower=0> n_a ; int<lower=0> n_b ; real<lower=0> a[n_a] ; real<lower=0> b[n_b] ; real<lower=0> c_mu_diff ; real<lower=0> c_es ; real<lower=0> c_cohenu ; real<lower=0> c_pod ; real<lower=0> c_pbt ; real<lower=0> cdash_pbt ; } parameters { real<lower=0> mu_a ; real<lower=0> sigma ; real<lower=0> mu_b ; } model { a ~ normal(mu_a, sigma) ; b ~ normal(mu_b, sigma) ; } generated quantities { vector[n_a] log_lik ; real mu_diff ; real es ; real cohenu ; real pod ; real pbt ; int<lower=0, upper=1> prob_mu_diff_upper_0 ; int<lower=0, upper=1> prob_mu_diff_upper_c ; int<lower=0, upper=1> prob_es_upper_c ; int<lower=0, upper=1> prob_cohenu_upper_c ; int<lower=0, upper=1> prob_pod_upper_c ; int<lower=0, upper=1> prob_pbt_upper_cdash ; for(i in 1:n_a){ log_lik[i] = normal_lpdf(a[i] | mu_a, sigma) + normal_lpdf(b[i] | mu_b, sigma) ; } mu_diff = mu_a - mu_b ; es = mu_diff / sigma ; cohenu = normal_cdf(mu_a, mu_b, sigma) ; pod = normal_cdf(es / sqrt(2), 0, 1) ; pbt = normal_cdf((mu_diff - c_pbt) / ( sqrt(2) * sigma ), 0, 1) ; prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ; prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ; prob_es_upper_c = es > c_es ? 1 : 0 ; prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ; prob_pod_upper_c = pod > c_pod ? 1 : 0 ; prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ; }
後はサンプリングを同様に行うだけ。
import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # 等分散モデル # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 pars = ["mu_a", "sigma", "mu_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング fit = model.sampling(data=stan_data, pars=pars, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=pars, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) # 事後分布の可視化 for par in summary_df.index[summary_df["sd"] == 0]: pars.remove(par) fit.traceplot(pars) plt.show() # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
これを実行すると、s.d.についてはEAPが8.63と推定されます。a群とb群のちょうど中間ぐらいですね。 両方の分布に同じ分散を使っているので、推定される値としては納得できるものです。 ただ尤度については716と不等分散モデルより高くなってしまいました。よってデータにはフィッティングしていないと捉えます。
変分ベイズによる推定
これまでは事後分布の推定を、NUTSアルゴリズムによるサンプリングで行ってきました。 ところでStanには、MCMCサンプリング以外にも変分ベイズによる分布の推定も一応実装されています。ついでにこれも試してみましょう。モデルのコンパイル部分までは、上述のものと共通です。
# サンプリングをADVIでやってみる import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display from collections import OrderedDict from pystan.external.pymc import plots %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # 等分散モデル # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 pars = ["mu_a", "sigma", "mu_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでADVIサンプリング fit = model.vb(data=stan_data, pars=pars, iter=100000, seed=1234) # 事後分布の表を取得 # ADVIはAPIがまだ完成されていないので、summaryの表を作る方式が違う # chainを利用したサンプリングもしていないため、Rhatも計算できない vb_sample = pd.read_csv(fit["args"]["sample_file"].decode("utf-8"), comment="#") vb_sample = vb_sample.drop("lp__", 1) summary_df = vb_sample.describe(percentiles=prob).T display(summary_df) # カーネル密度推定出来ないパラメータの削除 for par in summary_df.index[summary_df["std"] == 0]: pars.remove(par) ## 事後分布の可視化 od = OrderedDict() for i, par in enumerate(fit["sampler_param_names"]): par_s = par.split(".") if len(par_s) == 1: od[par] = np.array(fit["sampler_params"][i]) else: par = par_s[0] if par in od.keys(): od[par] = np.vstack([od[par], np.array(fit["sampler_params"][i])]) else: od[par] = np.array(fit["sampler_params"][i]) plots.traceplot(od, pars) plt.show() # WAICの計算 log_lik = od["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
注意すべきところは、PyStanの変分ベイズのAPIが完成されていないためsamplingとやり方がぜんぜん違う点です。 サンプリング結果の集計にはpandasのdescribeメソッドを使っています。
結果は次のとおりです。
count | mean | std | min | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | max | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu_a | 1001 | 56.4035 | 1.14076 | 52.6137 | 54.1746 | 54.6317 | 55.6578 | 56.3594 | 57.2304 | 58.3143 | 58.6232 | 59.8892 |
sigma | 1001 | 7.8148 | 0.593153 | 5.79781 | 6.71698 | 6.90312 | 7.43062 | 7.78376 | 8.18896 | 8.82866 | 9.04493 | 10.1979 |
mu_b | 1001 | 40.3831 | 1.19016 | 36.5435 | 38.1027 | 38.3704 | 39.5665 | 40.4121 | 41.1783 | 42.2879 | 42.7135 | 44.4097 |
log_lik.1 | 1001 | -6.42276 | 0.16823 | -7.06181 | -6.76429 | -6.6899 | -6.53466 | -6.42269 | -6.31037 | -6.15472 | -6.10705 | -5.86712 |
log_lik.2 | 1001 | -6.07869 | 0.150718 | -6.67179 | -6.3649 | -6.31169 | -6.17576 | -6.08422 | -5.98092 | -5.81884 | -5.78335 | -5.58263 |
log_lik.3 | 1001 | -8.86784 | 0.466876 | -10.957 | -9.91909 | -9.69209 | -9.1477 | -8.81625 | -8.52498 | -8.19867 | -8.08506 | -7.54151 |
log_lik.4 | 1001 | -6.23111 | 0.153535 | -6.66425 | -6.54282 | -6.48627 | -6.33461 | -6.2304 | -6.12679 | -5.96821 | -5.93276 | -5.71942 |
log_lik.5 | 1001 | -6.0363 | 0.152214 | -6.58138 | -6.32926 | -6.28372 | -6.13988 | -6.03141 | -5.93468 | -5.77738 | -5.73873 | -5.49402 |
log_lik.6 | 1001 | -9.86775 | 0.610265 | -12.5208 | -11.2156 | -11.0038 | -10.2397 | -9.81467 | -9.41468 | -9.01542 | -8.86263 | -8.08728 |
log_lik.7 | 1001 | -6.96474 | 0.214733 | -7.70444 | -7.41328 | -7.32943 | -7.10373 | -6.95857 | -6.81509 | -6.63545 | -6.57965 | -6.26646 |
log_lik.8 | 1001 | -7.96041 | 0.353478 | -10.3051 | -8.72111 | -8.56551 | -8.16825 | -7.93375 | -7.71524 | -7.44256 | -7.35968 | -7.04291 |
log_lik.9 | 1001 | -6.55416 | 0.179014 | -7.26925 | -6.9163 | -6.838 | -6.66726 | -6.55088 | -6.42996 | -6.27155 | -6.21939 | -6.00683 |
log_lik.10 | 1001 | -6.20992 | 0.155187 | -6.83945 | -6.50715 | -6.45275 | -6.31807 | -6.21339 | -6.10742 | -5.94728 | -5.89824 | -5.662 |
log_lik.11 | 1001 | -7.74777 | 0.313827 | -9.05947 | -8.43614 | -8.28447 | -7.94685 | -7.71843 | -7.52812 | -7.30075 | -7.21452 | -6.86988 |
log_lik.12 | 1001 | -6.23652 | 0.159109 | -6.74122 | -6.5479 | -6.49699 | -6.34863 | -6.22906 | -6.12679 | -5.9821 | -5.94032 | -5.73857 |
log_lik.13 | 1001 | -6.24006 | 0.160162 | -6.76222 | -6.54558 | -6.50324 | -6.35045 | -6.23568 | -6.134 | -5.98247 | -5.93949 | -5.70956 |
log_lik.14 | 1001 | -6.3109 | 0.159837 | -6.75209 | -6.61471 | -6.57248 | -6.42249 | -6.30973 | -6.19668 | -6.06357 | -6.02272 | -5.77449 |
log_lik.15 | 1001 | -6.76248 | 0.194958 | -7.7655 | -7.14874 | -7.08354 | -6.89143 | -6.75921 | -6.62693 | -6.45886 | -6.40091 | -6.21029 |
log_lik.16 | 1001 | -9.50787 | 0.558317 | -12.0117 | -10.7875 | -10.4896 | -9.84995 | -9.46429 | -9.09205 | -8.71395 | -8.58059 | -7.86643 |
log_lik.17 | 1001 | -6.23197 | 0.159713 | -6.77176 | -6.55355 | -6.51252 | -6.33496 | -6.22916 | -6.12582 | -5.97284 | -5.91782 | -5.47086 |
log_lik.18 | 1001 | -10.937 | 0.792785 | -16.4895 | -12.6988 | -12.3467 | -11.3937 | -10.8907 | -10.4062 | -9.72509 | -9.5612 | -8.87464 |
log_lik.19 | 1001 | -6.55719 | 0.178717 | -7.10137 | -6.92829 | -6.86605 | -6.66928 | -6.54786 | -6.43631 | -6.2802 | -6.22669 | -5.77744 |
log_lik.20 | 1001 | -6.92278 | 0.21248 | -7.83707 | -7.3648 | -7.28035 | -7.0519 | -6.91023 | -6.78217 | -6.60749 | -6.55017 | -6.33407 |
log_lik.21 | 1001 | -7.55595 | 0.295133 | -9.58726 | -8.17106 | -8.06674 | -7.72853 | -7.5413 | -7.3452 | -7.10152 | -7.04768 | -6.85252 |
log_lik.22 | 1001 | -7.55767 | 0.281643 | -9.09015 | -8.19547 | -8.05543 | -7.72132 | -7.53802 | -7.36901 | -7.14595 | -7.075 | -6.81468 |
log_lik.23 | 1001 | -7.35302 | 0.270502 | -8.89604 | -7.91456 | -7.82611 | -7.51521 | -7.32474 | -7.17505 | -6.95945 | -6.88067 | -6.60574 |
log_lik.24 | 1001 | -6.94368 | 0.215861 | -8.22401 | -7.36714 | -7.31273 | -7.08547 | -6.93466 | -6.79251 | -6.616 | -6.54842 | -6.35364 |
log_lik.25 | 1001 | -6.71402 | 0.191982 | -7.69412 | -7.09788 | -7.03066 | -6.84066 | -6.71253 | -6.58261 | -6.40948 | -6.34514 | -6.1755 |
log_lik.26 | 1001 | -9.53694 | 0.584227 | -13.4754 | -10.843 | -10.567 | -9.86342 | -9.49629 | -9.13708 | -8.68395 | -8.54657 | -8.1407 |
log_lik.27 | 1001 | -6.78402 | 0.201555 | -7.52225 | -7.20148 | -7.13041 | -6.91297 | -6.78015 | -6.64539 | -6.46899 | -6.42453 | -6.18105 |
log_lik.28 | 1001 | -6.44432 | 0.167778 | -6.96219 | -6.78804 | -6.72768 | -6.55666 | -6.44535 | -6.32604 | -6.18597 | -6.13832 | -5.85676 |
log_lik.29 | 1001 | -7.69239 | 0.315538 | -9.81127 | -8.36852 | -8.23455 | -7.87711 | -7.67732 | -7.4738 | -7.2188 | -7.15845 | -6.87702 |
log_lik.30 | 1001 | -7.7927 | 0.319736 | -9.13938 | -8.50667 | -8.34904 | -7.99599 | -7.76069 | -7.5676 | -7.33473 | -7.24417 | -6.90979 |
log_lik.31 | 1001 | -6.16013 | 0.157118 | -6.71291 | -6.45853 | -6.41783 | -6.26426 | -6.15737 | -6.05708 | -5.90701 | -5.86343 | -5.49712 |
log_lik.32 | 1001 | -6.13017 | 0.151897 | -6.74355 | -6.41578 | -6.37076 | -6.23383 | -6.13207 | -6.0305 | -5.87105 | -5.83521 | -5.59475 |
log_lik.33 | 1001 | -7.00073 | 0.223585 | -7.81207 | -7.48043 | -7.40014 | -7.14758 | -6.99534 | -6.84432 | -6.65199 | -6.60361 | -6.3401 |
log_lik.34 | 1001 | -6.92585 | 0.215031 | -7.65724 | -7.38092 | -7.30914 | -7.06934 | -6.91029 | -6.78167 | -6.59183 | -6.54201 | -6.25093 |
log_lik.35 | 1001 | -6.78402 | 0.201555 | -7.52225 | -7.20148 | -7.13041 | -6.91297 | -6.78015 | -6.64539 | -6.46899 | -6.42453 | -6.18105 |
log_lik.36 | 1001 | -6.45221 | 0.172966 | -7.0252 | -6.79691 | -6.752 | -6.56508 | -6.448 | -6.33613 | -6.17511 | -6.12069 | -5.69448 |
log_lik.37 | 1001 | -6.43874 | 0.169585 | -6.94675 | -6.7849 | -6.73926 | -6.54679 | -6.43412 | -6.32283 | -6.16711 | -6.11841 | -5.65075 |
log_lik.38 | 1001 | -5.98581 | 0.149329 | -6.52304 | -6.27182 | -6.22623 | -6.0843 | -5.98738 | -5.88368 | -5.74347 | -5.69836 | -5.48248 |
log_lik.39 | 1001 | -6.22339 | 0.155309 | -6.81162 | -6.52915 | -6.47324 | -6.33043 | -6.22686 | -6.12094 | -5.9492 | -5.9203 | -5.72975 |
log_lik.40 | 1001 | -6.38489 | 0.163779 | -6.94275 | -6.72587 | -6.64711 | -6.49179 | -6.38321 | -6.27196 | -6.11557 | -6.0736 | -5.90043 |
log_lik.41 | 1001 | -6.66083 | 0.18857 | -7.36905 | -7.064 | -6.97551 | -6.77774 | -6.65256 | -6.5262 | -6.36894 | -6.30808 | -6.14803 |
log_lik.42 | 1001 | -6.03561 | 0.149328 | -6.59567 | -6.32253 | -6.27471 | -6.13574 | -6.03972 | -5.93328 | -5.79028 | -5.74892 | -5.56119 |
log_lik.43 | 1001 | -7.9404 | 0.350003 | -10.3391 | -8.68739 | -8.53891 | -8.1428 | -7.92091 | -7.69262 | -7.41976 | -7.33845 | -7.05618 |
log_lik.44 | 1001 | -11.5909 | 0.862142 | -15.4103 | -13.5076 | -13.1583 | -12.1299 | -11.5137 | -10.9734 | -10.3584 | -10.1675 | -9.0104 |
log_lik.45 | 1001 | -7.77097 | 0.315238 | -9.01645 | -8.47879 | -8.3403 | -7.97236 | -7.74171 | -7.5439 | -7.3051 | -7.22634 | -6.98188 |
log_lik.46 | 1001 | -7.6209 | 0.303068 | -9.7381 | -8.25764 | -8.14321 | -7.80149 | -7.6036 | -7.40669 | -7.16314 | -7.0904 | -6.8879 |
log_lik.47 | 1001 | -6.97629 | 0.217257 | -7.94213 | -7.42008 | -7.33636 | -7.10676 | -6.9649 | -6.82924 | -6.64425 | -6.58657 | -6.38334 |
log_lik.48 | 1001 | -6.66957 | 0.186988 | -7.39538 | -7.06688 | -6.98032 | -6.78594 | -6.65865 | -6.54712 | -6.38723 | -6.3303 | -6.12713 |
log_lik.49 | 1001 | -6.63225 | 0.187684 | -7.31207 | -7.0097 | -6.94733 | -6.75122 | -6.63039 | -6.50275 | -6.33472 | -6.2794 | -6.0357 |
log_lik.50 | 1001 | -6.32385 | 0.159621 | -6.76867 | -6.63689 | -6.58101 | -6.43788 | -6.32294 | -6.2061 | -6.07091 | -6.03085 | -5.80358 |
mu_diff | 1001 | 16.0204 | 1.66137 | 10.9645 | 12.6538 | 13.2746 | 14.9254 | 16.0157 | 17.0902 | 18.7446 | 19.2714 | 20.8654 |
es | 1001 | 2.06108 | 0.259806 | 1.27133 | 1.58195 | 1.65403 | 1.87671 | 2.05535 | 2.23764 | 2.50222 | 2.59257 | 2.90178 |
cohenu | 1001 | 0.977041 | 0.0139309 | 0.898194 | 0.94317 | 0.95094 | 0.969721 | 0.980078 | 0.987378 | 0.993829 | 0.995237 | 0.998145 |
pod | 1001 | 0.924156 | 0.0256932 | 0.815664 | 0.868346 | 0.878915 | 0.907751 | 0.926937 | 0.943204 | 0.961581 | 0.966615 | 0.979909 |
pbt | 1001 | 0.705715 | 0.0529771 | 0.538187 | 0.593225 | 0.61623 | 0.671596 | 0.708586 | 0.741193 | 0.786155 | 0.803224 | 0.849742 |
prob_mu_diff_upper_0 | 1001 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_mu_diff_upper_c | 1001 | 0.879121 | 0.32615 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_es_upper_c | 1001 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
prob_cohenu_upper_c | 1001 | 0.951049 | 0.215874 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_pod_upper_c | 1001 | 0.152847 | 0.36002 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
prob_pbt_upper_cdash | 1001 | 0.97003 | 0.17059 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その4
その4です。今回は第3章の話をしていきます。
3章 独立した2群の差の推測
内容
- 2群に分けた群比較をする際には、ランダム化による交絡因子の排除が重要!(ランダマイゼーション)
- 2群の、どんな差を見たいかで仮説は変わる
- 群1の平均が群2の平均を上回る確率
- 群1, 2の平均値の差のt年推定、区間推定
- 群1, 2の平均値の差が基準点cより大きい確率
- 効果量が基準点cより大きい確率
- 非重複度、優越率、閾上率…etc
- 2群でモデルを作る時に、標準偏差は共通させるか?分けるか?
- 共通させる場合
- 独立に定義する場合
- かならずしも独立に定義したほうが良いわけではない。モデルの複雑性は増してしまうので。
- 共通させる場合
- 生成量の例
- モデルの選択
章末問題
データの取得
とりあえずデータを読み込みます。 今回は2群の対応のない場合の比較です。
import numpy as np from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger # データの準備 a = np.array([56,55,55,62,54,63,47,58,56,56,57,52,53,50,50,57,57,55,60,65,53,43,60,51,52, 60,54,49,56,54,55,57,53,58,54,57,60,57,53,61,60,58,56,52,62,52,66,63,54,50]) b = np.array([33,37,59,41,42,61,46,25,32,35,55,44,45,41,33,61,46,16,48,34,27,37,28,31,32, 20,50,42,26,55,45,36,51,51,50,48,47,39,36,35,32,38,25,66,54,27,35,34,49,39]) logger = get_logger()
各種統計量の計算と可視化
- 各種統計量はpandas.describeを使えば大体出せますね。
- 分布はdistplot以外にも色々ありますが、matplotlibよりseabornの方が簡単に色々でるので使いました。
import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from IPython.core.display import display %matplotlib inline df = pd.DataFrame({"a":a, "b":b}) display(df.describe()) # 分散は出ない sns.distplot(a) plt.show() sns.distplot(b) plt.show()
- 結果
a | b | |
---|---|---|
count | 50 | 50 |
mean | 55.76 | 40.38 |
std | 4.57839 | 11.1427 |
min | 43 | 16 |
25% | 53 | 33 |
50% | 56 | 39 |
75% | 58 | 48 |
max | 66 | 66 |
さてこの結果から何が言えるかというと…
- a群の方が平均値が15ぐらい高くなっている
- 中央値となる50%点もa群の方が高いが、なんといっても最小値が大幅に違う
- この2つの群はともに正規分布っぽい感じにはなっている
- 分散は2群で結構違いそう。
みたいなものがわかります。よくある統計でしたら、t検定で2群の間には有意差があるということを出して終わりですね。 今回は最初に書いた内容を使って、
- 差の大きさがどれぐらい幅があるのか(ちょっとだけの差でしたら、殆ど意味がないですよね)。
- 差は14ぐらいありそうである。それはどれぐらい確かなのか。
- 2つの分布はどれぐらい離れているのか。重複してる部分はどれぐらいあるのか。
- 群としてはaの方が大きくても、この群から適当に1つずつ選んだ時に、どれぐらいの確立でaの方の個体が大きな値になるのか
と言ったことを調べていくわけです。
Stanによる統計モデルの構築
Stanでモデルを作って、2群を比較します。 教科書に従って、とりあえず不偏分散のモデルを作っています。aとbの分散値が先程の要約テーブルで結構違いましたので、等分散よりモデルとしては適当になりそうです(後で検証します)
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "g2_ind.stan") stan_file_c = os.path.join("stan", "g2_ind.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f)
# Stanのモデル data { int<lower=0> n_a ; int<lower=0> n_b ; real<lower=0> a[n_a] ; real<lower=0> b[n_b] ; real<lower=0> c_mu_diff ; real<lower=0> c_es ; real<lower=0> c_cohenu ; real<lower=0> c_pod ; real<lower=0> c_pbt ; real<lower=0> cdash_pbt ; } parameters { real<lower=0> mu_a ; real<lower=0> sigma_a ; real<lower=0> mu_b ; real<lower=0> sigma_b ; } model { a ~ normal(mu_a, sigma_a) ; b ~ normal(mu_b, sigma_b) ; } generated quantities { vector[n_a] log_lik ; real mu_diff ; real es ; real cohenu ; real pod ; real pbt ; int<lower=0, upper=1> prob_mu_diff_upper_0 ; int<lower=0, upper=1> prob_mu_diff_upper_c ; int<lower=0, upper=1> prob_es_upper_c ; int<lower=0, upper=1> prob_cohenu_upper_c ; int<lower=0, upper=1> prob_pod_upper_c ; int<lower=0, upper=1> prob_pbt_upper_cdash ; for(i in 1:n_a){ log_lik[i] = normal_lpdf(a[i] | mu_a, sigma_a) + normal_lpdf(b[i] | mu_b, sigma_b) ; } mu_diff = mu_a - mu_b ; es = mu_diff / sigma_a ; cohenu = normal_cdf(mu_a, mu_b, sigma_b) ; pod = normal_cdf(mu_diff / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; pbt = normal_cdf((mu_diff - c_pbt) / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ; prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ; prob_es_upper_c = es > c_es ? 1 : 0 ; prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ; prob_pod_upper_c = pod > c_pod ? 1 : 0 ; prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ; }
統計モデルによる事後分布のサンプリング
注意点は2つあります。
まず1つは可視化をする部分です。
PyStanのモジュールをそのまま使ったサンプリング値の可視化では、サンプリング値が全て同値である場合(=分散が0の時)はtraceplotでエラーが発生してしまいます。
これは可視化するために行っている、カーネル密度推定の行列計算で失敗するためです。
今回、僕は自分でPyStanの一部を改造したものを利用しました。
が、これは試してみたかっただけなので、こんな面倒なことをせず、summary_df
のsd columnを見て、0だったらtrace plotを書かなければ良いと思います。
もう1つはWAICの計算の部分です。 計算をする際に、データ数分だけfor文を回して対数尤度を計算しています この実装をしているものが検索してみると見つかるし、豊田本上級でもこの紹介がされています(定義的にも正しい?)。 ただ今回の場合は、モデルが単純で50個のデータが共通の母数を持つ正規分布から生成されていので、対数尤度は1次元配列で計算し、waicの計算でsumをする部分が無くても大丈夫です。 本の付録のモデルでは、こちらの簡潔版で計算されています。
import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 par = ["mu_a", "sigma_a", "mu_b", "sigma_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) # 事後分布の可視化 fit.traceplot(par, {"prob_mu_diff_upper_0":np.int}) plt.show() # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu_a | 55.7604 | 0.00210572 | 0.665888 | 54.4483 | 54.6698 | 55.3161 | 55.7601 | 56.2025 | 56.8535 | 57.0771 | 100000 | 0.99999 |
sigma_a | 4.69936 | 0.00163133 | 0.492013 | 3.85845 | 3.97235 | 4.35244 | 4.65696 | 5.00083 | 5.56823 | 5.78077 | 90964 | 0.999993 |
mu_b | 40.3768 | 0.00513303 | 1.62321 | 37.1566 | 37.7085 | 39.2942 | 40.3763 | 41.4601 | 43.0385 | 43.5743 | 100000 | 0.999979 |
sigma_b | 11.4351 | 0.00377584 | 1.19403 | 9.39226 | 9.66662 | 10.5891 | 11.3316 | 12.1638 | 13.5548 | 14.0673 | 100000 | 1.00003 |
log_lik[0] | -6.04735 | 0.000490513 | 0.152051 | -6.35677 | -6.30426 | -6.14733 | -6.04291 | -5.94225 | -5.80452 | -5.76237 | 96090 | 0.999988 |
log_lik[1] | -5.88977 | 0.000482951 | 0.146824 | -6.18808 | -6.13823 | -5.98677 | -5.88593 | -5.78838 | -5.65612 | -5.61159 | 92424 | 1.00001 |
log_lik[2] | -7.21332 | 0.000987215 | 0.312185 | -7.8894 | -7.76195 | -7.40848 | -7.19029 | -6.99224 | -6.74368 | -6.67145 | 100000 | 0.999982 |
log_lik[3] | -6.74241 | 0.000737068 | 0.233082 | -7.24164 | -7.15023 | -6.89021 | -6.72715 | -6.57867 | -6.38654 | -6.32838 | 100000 | 0.999992 |
log_lik[4] | -5.91407 | 0.000486231 | 0.147301 | -6.21372 | -6.16309 | -6.01166 | -5.90981 | -5.8126 | -5.67868 | -5.63577 | 91775 | 1.00002 |
log_lik[5] | -8.73421 | 0.00139258 | 0.440372 | -9.66647 | -9.49826 | -9.01462 | -8.7067 | -8.42342 | -8.06105 | -7.95471 | 100000 | 0.999977 |
log_lik[6] | -7.74988 | 0.00123701 | 0.391176 | -8.59949 | -8.44025 | -7.99428 | -7.7188 | -7.47338 | -7.16121 | -7.07202 | 100000 | 1.00001 |
log_lik[7] | -6.88161 | 0.000752406 | 0.237932 | -7.39441 | -7.29783 | -7.03121 | -6.86651 | -6.71561 | -6.51828 | -6.45829 | 100000 | 0.999984 |
log_lik[8] | -6.10953 | 0.000492635 | 0.155785 | -6.42726 | -6.37336 | -6.21194 | -6.10435 | -6.00181 | -5.86187 | -5.81865 | 100000 | 0.999984 |
log_lik[9] | -5.94667 | 0.00048357 | 0.148166 | -6.24786 | -6.19765 | -6.04461 | -5.9428 | -5.8443 | -5.70916 | -5.66716 | 93881 | 0.999998 |
log_lik[10] | -6.71085 | 0.000707805 | 0.223828 | -7.19011 | -7.1025 | -6.85252 | -6.6971 | -6.55379 | -6.36867 | -6.30991 | 100000 | 0.999989 |
log_lik[11] | -6.21359 | 0.000504989 | 0.159691 | -6.54217 | -6.48589 | -6.31792 | -6.20807 | -6.10355 | -5.96017 | -5.91534 | 100000 | 1.00002 |
log_lik[12] | -6.09369 | 0.00049108 | 0.150984 | -6.40111 | -6.34936 | -6.19326 | -6.08984 | -5.98977 | -5.85288 | -5.80928 | 94527 | 1.00002 |
log_lik[13] | -6.60844 | 0.000670976 | 0.212181 | -7.05939 | -6.97869 | -6.7431 | -6.59501 | -6.461 | -6.28327 | -6.22806 | 100000 | 1.00001 |
log_lik[14] | -6.82174 | 0.000683128 | 0.216024 | -7.28002 | -7.19688 | -6.95918 | -6.80882 | -6.67223 | -6.48941 | -6.43331 | 100000 | 1.00001 |
log_lik[15] | -7.5454 | 0.0011654 | 0.368531 | -8.34713 | -8.19487 | -7.77589 | -7.51674 | -7.28485 | -6.99185 | -6.90637 | 100000 | 0.999982 |
log_lik[16] | -5.99181 | 0.000501445 | 0.14957 | -6.29642 | -6.24462 | -6.09066 | -5.9877 | -5.8892 | -5.7523 | -5.71028 | 88970 | 1.00003 |
log_lik[17] | -8.19007 | 0.00156933 | 0.496266 | -9.27191 | -9.06673 | -8.49926 | -8.15016 | -7.83863 | -7.44239 | -7.33214 | 100000 | 1 |
log_lik[18] | -6.48041 | 0.000548936 | 0.173589 | -6.83737 | -6.77563 | -6.59425 | -6.47343 | -6.35979 | -6.20713 | -6.15985 | 100000 | 1.00001 |
log_lik[19] | -7.98662 | 0.00135638 | 0.428926 | -8.91464 | -8.74592 | -8.25546 | -7.95257 | -7.68182 | -7.34138 | -7.24086 | 100000 | 0.999979 |
log_lik[20] | -6.71568 | 0.000647957 | 0.204902 | -7.14689 | -7.07082 | -6.84698 | -6.70495 | -6.57247 | -6.39928 | -6.34525 | 100000 | 0.999987 |
log_lik[21] | -9.68219 | 0.00249132 | 0.787825 | -11.3858 | -11.0712 | -10.1776 | -9.6205 | -9.12183 | -8.49679 | -8.32186 | 100000 | 1.00001 |
log_lik[22] | -6.85578 | 0.000652113 | 0.206216 | -7.2915 | -7.21196 | -6.98743 | -6.84574 | -6.71229 | -6.53602 | -6.47974 | 100000 | 0.999978 |
log_lik[23] | -6.70807 | 0.000604833 | 0.191265 | -7.10608 | -7.03765 | -6.83147 | -6.69943 | -6.57574 | -6.40914 | -6.35915 | 100000 | 1 |
log_lik[24] | -6.43882 | 0.000529426 | 0.167419 | -6.78252 | -6.724 | -6.54813 | -6.43286 | -6.32297 | -6.17362 | -6.12905 | 100000 | 0.999994 |
log_lik[25] | -7.88994 | 0.0011714 | 0.37043 | -8.69871 | -8.54522 | -8.12159 | -7.86367 | -7.6282 | -7.33072 | -7.24402 | 100000 | 0.999994 |
log_lik[26] | -6.26907 | 0.000533125 | 0.164099 | -6.60602 | -6.54939 | -6.37576 | -6.26344 | -6.15559 | -6.00966 | -5.96361 | 94744 | 1.00001 |
log_lik[27] | -6.90996 | 0.000818656 | 0.258882 | -7.4684 | -7.36431 | -7.0731 | -6.8915 | -6.72761 | -6.51943 | -6.45387 | 100000 | 1.00001 |
log_lik[28] | -6.64837 | 0.000691255 | 0.218594 | -7.11348 | -7.02718 | -6.7864 | -6.63571 | -6.49597 | -6.31147 | -6.25681 | 100000 | 0.999982 |
log_lik[29] | -6.74747 | 0.000709018 | 0.224211 | -7.22664 | -7.13774 | -6.88949 | -6.73451 | -6.59022 | -6.40498 | -6.34879 | 100000 | 0.999986 |
log_lik[30] | -5.92905 | 0.000495427 | 0.148163 | -6.23007 | -6.17892 | -6.02743 | -5.92463 | -5.82661 | -5.69268 | -5.64909 | 89438 | 1.00002 |
log_lik[31] | -5.9427 | 0.000481788 | 0.147672 | -6.24341 | -6.19239 | -6.04019 | -5.93912 | -5.84076 | -5.70632 | -5.66325 | 93948 | 1 |
log_lik[32] | -6.45466 | 0.000551094 | 0.174271 | -6.81692 | -6.7533 | -6.56731 | -6.4479 | -6.33399 | -6.1815 | -6.13349 | 100000 | 1 |
log_lik[33] | -6.39362 | 0.00054583 | 0.172607 | -6.75128 | -6.68942 | -6.50549 | -6.38655 | -6.27466 | -6.12209 | -6.07451 | 100000 | 1.00001 |
log_lik[34] | -6.26907 | 0.000533125 | 0.164099 | -6.60602 | -6.54939 | -6.37576 | -6.26344 | -6.15559 | -6.00966 | -5.96361 | 94744 | 1.00001 |
log_lik[35] | -6.09634 | 0.000515248 | 0.154159 | -6.40996 | -6.35791 | -6.19751 | -6.09192 | -5.99059 | -5.85023 | -5.80673 | 89517 | 1.00002 |
log_lik[36] | -6.42419 | 0.000541343 | 0.171188 | -6.77728 | -6.71547 | -6.53626 | -6.41756 | -6.30555 | -6.15428 | -6.1082 | 100000 | 1.00001 |
log_lik[37] | -5.87456 | 0.000486567 | 0.147509 | -6.1755 | -6.12359 | -5.97152 | -5.87056 | -5.77306 | -5.63864 | -5.59665 | 91908 | 1.00002 |
log_lik[38] | -6.08502 | 0.000474243 | 0.149969 | -6.38969 | -6.33788 | -6.18427 | -6.08081 | -5.98184 | -5.84684 | -5.80406 | 100000 | 1 |
log_lik[39] | -6.58677 | 0.000616932 | 0.195091 | -6.9986 | -6.92179 | -6.71234 | -6.57743 | -6.45114 | -6.2843 | -6.23202 | 100000 | 0.999988 |
log_lik[40] | -6.52813 | 0.000557474 | 0.176289 | -6.89347 | -6.82835 | -6.64266 | -6.52131 | -6.40592 | -6.25027 | -6.20306 | 100000 | 0.999983 |
log_lik[41] | -5.97066 | 0.00048391 | 0.148823 | -6.27342 | -6.2218 | -6.06889 | -5.9667 | -5.86849 | -5.73262 | -5.69075 | 94583 | 1.00001 |
log_lik[42] | -6.7658 | 0.000749145 | 0.2369 | -7.27327 | -7.17892 | -6.91408 | -6.75068 | -6.60047 | -6.40354 | -6.34414 | 100000 | 0.999985 |
log_lik[43] | -8.75266 | 0.00173329 | 0.548115 | -9.93939 | -9.71752 | -9.09476 | -8.70939 | -8.36307 | -7.93373 | -7.80486 | 100000 | 0.999986 |
log_lik[44] | -7.4732 | 0.000864603 | 0.273411 | -8.05226 | -7.94825 | -7.64826 | -7.45718 | -7.28092 | -7.05404 | -6.98571 | 100000 | 0.99998 |
log_lik[45] | -6.86811 | 0.000670866 | 0.212146 | -7.31424 | -7.2334 | -7.00434 | -6.85678 | -6.71946 | -6.53968 | -6.48454 | 100000 | 0.999993 |
log_lik[46] | -8.39545 | 0.00163259 | 0.516269 | -9.51436 | -9.31025 | -8.71891 | -8.35481 | -8.02785 | -7.61944 | -7.50049 | 100000 | 0.999979 |
log_lik[47] | -7.21641 | 0.00090867 | 0.287347 | -7.83597 | -7.72231 | -7.398 | -7.19558 | -7.01405 | -6.78238 | -6.71172 | 100000 | 0.99998 |
log_lik[48] | -6.19707 | 0.000517114 | 0.15826 | -6.52008 | -6.46738 | -6.30066 | -6.19209 | -6.08782 | -5.94499 | -5.90124 | 93663 | 1.00001 |
log_lik[49] | -6.61441 | 0.000670758 | 0.212112 | -7.0674 | -6.98435 | -6.74872 | -6.60168 | -6.46637 | -6.28944 | -6.23353 | 100000 | 1.00001 |
mu_diff | 15.3836 | 0.00555407 | 1.75635 | 11.9336 | 12.5127 | 14.207 | 15.3843 | 16.5573 | 18.2745 | 18.853 | 100000 | 0.999983 |
es | 3.30844 | 0.00160306 | 0.506931 | 2.36452 | 2.50846 | 2.95786 | 3.28843 | 3.63776 | 4.17218 | 4.35433 | 100000 | 0.99998 |
cohenu | 0.90849 | 0.000107709 | 0.0335828 | 0.831913 | 0.847353 | 0.888497 | 0.912552 | 0.932937 | 0.955767 | 0.961685 | 97214 | 0.999988 |
pod | 0.891195 | 0.00010542 | 0.0333368 | 0.816538 | 0.831072 | 0.87091 | 0.89468 | 0.915269 | 0.939324 | 0.945675 | 100000 | 0.999983 |
pbt | 0.667754 | 0.000166733 | 0.0527255 | 0.560294 | 0.578451 | 0.632602 | 0.669254 | 0.704582 | 0.751648 | 0.766491 | 100000 | 0.999974 |
prob_mu_diff_upper_0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_mu_diff_upper_c | 0.78579 | 0.00138494 | 0.410275 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 87758 | 0.999979 |
prob_es_upper_c | 0.72119 | 0.00151356 | 0.448416 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 87773 | 0.999991 |
prob_cohenu_upper_c | 0.0856 | 0.000957681 | 0.279774 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 85344 | 0.999981 |
prob_pod_upper_c | 0.01434 | 0.00040641 | 0.118889 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 85576 | 0.999976 |
prob_pbt_upper_cdash | 0.89648 | 0.00108438 | 0.304638 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 78924 | 1.00001 |
問題なく、パラメータ値が収束して計算されました(Rhat<1.1で判定しています)
WAIC = 683.147として計算されます。 なお等分散モデルだとWAIC=716.444なので、実際に不偏分散モデルのほうが感覚としても、情報量としてもデータとモデルがフィットしていることがわかります。
仮説と検証
- A群の平均値がB群より高い確率は?
- prob_mu_diff_upper_0のサンプリング結果より、100%の確率で、A群の方がB群より高いことがわかります。
- AぐんとB群の平均値の差の点推定、95%信頼区間推定
- 15.384(1.756)[11.934, 18.853]となりました。
- 平均値の差の片側95%信頼区間推定の上限、下限はどのあたりか
- mu_diff > 12.513、あるいはmu_diff < 18.275の区間に95%の確率で入る
- 平均値の差が14より大きい確率は?
- prob_mu_diff_upper_cより、78.6%の確率で、14より大きいことがわかります。
- A群から見たB群の効果量の点推定、95%信頼区間推定 両側/片側
- A群から見たB群の効果量が3より大きい確率は?
- prob_es_upper_cより、72.1%の確率で、3より大きいことがわかります。
- A群の分布とB群の分布の重複する割合(非重複度)の点推定、95%信頼区間推定 両側/片側
- 非重複度が0.95より大きい確率
- prob_cohenu_upper_cより、8.6% の確率で、0.95より
- 無作為に選んだ際に、A群の値がB群より大きくなる確率(優越率)はどのぐらいか
- 0.891(0.033)[0.817, 0.946]の確率で高くなることが推定されます。
- 優越率が0.95より大きい確率
- prob_pod_upper_cより 1.4%の確率で0.95より高くなることが推定されます。
- 無作為に選んだそれぞれの測定値の差が10より大きい確率(閾上率)
- 0.668(0.053)[0.560, 0.766]の確率で、大きいことが推定されます。
- 閾上率が0.60(60%)より大きい確率
- prob_pbt_upper_cdashより89.6%の確率で60%より高くなることが推定されます。
以上のように、MCMCによるサンプリングの試行結果を使うことで、様々な情報を抜き出すことが出来ました。 次回は第4章の、対応のあるt検定です。