はじめての 統計データ分析 ―ベイズ的〈ポスト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検定です。
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その3
その3です。今回は第2章の話をしていきます。
2章 MCMCと正規分布の推測
内容
- 解析的に事後分布を推定するのは難しい!
事後分布に従う乱数を発生させる、という発想転換を行う(マルコフ連鎖モンテカルロ法; MCMC)
- 乱数群の結果を比較して、解析する
- HM法とGS法が具体的なアルゴリズムとしては代表的な存在
- GS法は特定の事前分布であることが必要(やや恣意的)
- HM法は自由な事前分布を使うことが出来る
- 事前分布の設定 : 常識的な知識から決める
- 乱数の発生 : 初期乱数は使用しない(burn in、あるいはwarm upと呼ばれる)
- 値が登ったり下りたり(ドリフトした)形状が観測される場合には、好ましくない乱数生成となる
- 複数chainを発生させて、全部で統一的な傾向が見られてない場合、好ましくない乱数生成が行われていると考えられる
- 収束判定指標(Rcap) : 複数chainの分散, 1.1以下が理想
- 有効標本数(neff) : 乱数の内、いくつが無関係な乱数として扱うことが出来るかを示す指標
- (以上2つの指標の導出が気になった場合、基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門に結構丁寧に書いてあるので、そちらをお薦めします。)
事後分布の解釈をどう行うか?
- 点推定量
- 確信区間(信用区間)
- 予測分布
- 新たなデータが得られうる分布のこと
- 事後予測分布
- 推定された母数の事後分布を用いて、未知データを予測する統計モデルの平均を計算
- 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
- 条件付き予測分布
- 事後分布の特定の点推定量を利用する
- 簡単だが、せっかく推定したモデルを有効活用出来ない
- (これについては、その1で紹介した「はじめての統計データ分析」 豊田秀樹のメモ - StatModeling Memorandumさんの方で問題点が指摘されています)
- (が、豊田先生もこのような指摘を結構受けたのか、本に付録でついてくる追加Q&Aに解説が追記されています(後述))
- ベイズ的推測
- 何らかの仮説を定めて、データから推定すること
- 統計的仮説検定では帰無仮説に基づいたが、そんなことはしなくてよいのがベイズ統計
- 生成量
- 推定した母数を利用して、各種統計量の分布も計算することが出来る
- 母数より生成される値なので、生成量という
- 推定した母数を利用して、各種統計量の分布も計算することが出来る
- 研究仮説が正しい確率
- 真偽を判定する二値変数を利用して、研究仮説が正しい確率を評価することが出来る
- 0 or 1しか使わないので、分散や信頼区間は使えない
条件付き予測分布の取扱について(2017-05-28追記)
ここまでで書いた概要部分でも説明しましたが、議論のされている部分なので補足説明です。 このテキストではデータが得られた状態で、次に得られるであろう未知データについての予測方法を2つ紹介しています。
1つは事後予測分布です。事後予測分布の数式はと記述されます。これを踏まえたやり方です。 数式だけ見ると混乱しがちなのですが、MCMCで計算する時にはそんなに難しくありません。 MCMCの各サンプリングにおいて得られるがを示していますし、このを構築したモデルに代入して得られる値が分布からサンプリングされたうちの1つということになります。 MCMCが積分の代わりをしている、という認識ですね。
もう1つは条件付き予測分布です。これはMCMCの結果として得られるのEAPやMAPといった点推定値をモデルに代入して、予測分布を得るというものです。 これは理解もし易いし、簡単です。分布の中で最もそれっぽいと理解出来た値をモデルに代入して予測をすることができます。 著者の豊田先生としては、
- 理解が簡単
- 速度が速い(生成量や予測量として計算しなくて済む)
- オンライン性(事後予測分布は過去データが必要だが、条件付き予測分布は点推定値だけでよい)
という利点を挙げています。
しかしながらリンク先にもあるように、条件付き予測分布には色々問題点もありますし、何より個人的には「もったいない」という点が目立つと考えています。 せっかくMCMCまでして分布を計算したんですから、それを有効に活用しましょう。 それに「仮説が正しい確率」を計算するには、どっちにしろ事後予想分布が必要になります。
以上の理由で、以後の説明では条件付き予測分布は使用せず、事後予想分布で推定を行っています。
章末問題
データの準備
- 1章のデータをそのまま利用します
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 # データの準備 x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33]) logger = get_logger()
Stanによる統計モデルの構築
- MCMCをStanで行います
- Stan言語で書かれたファイルを読み込んで、C++にコンパイルすることで、高速にサンプリングを行うことが可能となります。
- 標本データや問題文を見た時に、値の範囲について情報が全く得られなかったので、パラメータ分布には0以上という制限しか利用していません。
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "g1_mean.stan") stan_file_c = os.path.join("stan", "g1_mean.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 ; real<lower=0> x[n] ; int<lower=0> c ; real<lower=0, upper=1> percent ; real ces ; } parameters { real<lower=0> mu ; real<lower=0> sigma ; } model { x ~ normal(mu,sigma); } generated quantities{ real xaste ; real sigma2 ; real cv ; real es; real percent_point ; real<lower=0, upper=1> prob_under_c_cdf ; real rate_c ; real<lower=0, upper=1> prob_mu_under_c ; real<lower=0, upper=1> prob_under_c_bool ; real<lower=0, upper=1> prob_es_under_ces ; real<lower=0, upper=1> prob_prob_under_c_bool_over_percent ; xaste = normal_rng(mu,sigma) ; sigma2 = sigma * sigma ; cv = sigma / mu ; es = (mu - c) / sigma ; percent_point = mu + normal_cdf(percent, 0, 1) * sigma ; rate_c = xaste / c ; prob_under_c_cdf = normal_cdf(c, mu, sigma) ; prob_under_c_bool = xaste < c ? 1 : 0 ; prob_mu_under_c = mu < c ? 1 : 0 ; prob_es_under_ces = es < ces ? 1 : 0 ; prob_prob_under_c_bool_over_percent = prob_under_c_bool > percent ? 1 : 0 ; }
統計モデルによる事後分布のサンプリング
- コンパイルしたモデルを使って、正規分布からデータが得られた際に、正規分布の母数がどのように分布するか調べます
- 母数の分布から、各種効果量を推定しています。
- 仮説が得られる確率については、各サンプリング値において判定を行い、その結果のbool値の頻度で確認します。
- 付録でついてきたスクリプトにはここらへんの記述がありませんので、何か問題が有った場合には指摘頂けると助かります。
import pandas as pd import pickle import pystan import matplotlib import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 40) # Stanで使うデータの用意 stan_data = {"n": x.size, "x": x, "c": 41, "percent": 0.70, "ces": -0.05} # 興味のあるパラメータの設定 par = ["mu", "sigma", "xaste", "sigma2", "cv", "es", "percent_point", "rate_c", "prob_under_c_cdf", "prob_under_c_bool", "prob_mu_under_c", "prob_es_under_ces", "prob_prob_under_c_bool_over_percent"] # 興味のある母数の%点を書く prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # コンパイルしたモデルの読み込み stan_file_c = os.path.join("stan", "g1_mean.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="HMC") # 事後分布の表を取得、可視化 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.plot(pars=par) plt.show()
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu | 40.660 | 0.010 | 0.649 | 39.383 | 39.603 | 40.224 | 40.659 | 41.091 | 41.735 | 41.947 | 4120.000 | 1.004 |
sigma | 6.523 | 0.014 | 0.468 | 5.660 | 5.797 | 6.200 | 6.499 | 6.828 | 7.321 | 7.488 | 1154.000 | 1.008 |
xaste | 40.693 | 0.021 | 6.560 | 27.784 | 29.927 | 36.314 | 40.679 | 45.108 | 51.424 | 53.495 | 96434.000 | 1.000 |
sigma2 | 42.772 | 0.182 | 6.178 | 32.033 | 33.604 | 38.445 | 42.237 | 46.623 | 53.602 | 56.063 | 1153.000 | 1.008 |
cv | 0.160 | 0.000 | 0.012 | 0.139 | 0.142 | 0.152 | 0.160 | 0.168 | 0.181 | 0.185 | 1372.000 | 1.006 |
es | -0.052 | 0.002 | 0.099 | -0.248 | -0.215 | -0.120 | -0.053 | 0.014 | 0.112 | 0.144 | 4065.000 | 1.004 |
percent_point | 45.605 | 0.019 | 0.744 | 44.204 | 44.418 | 45.099 | 45.583 | 46.092 | 46.862 | 47.132 | 1509.000 | 1.009 |
rate_c | 0.993 | 0.001 | 0.160 | 0.678 | 0.730 | 0.886 | 0.992 | 1.100 | 1.254 | 1.305 | 96434.000 | 1.000 |
prob_under_c_cdf | 0.521 | 0.001 | 0.039 | 0.443 | 0.455 | 0.494 | 0.521 | 0.548 | 0.585 | 0.598 | 4055.000 | 1.004 |
prob_under_c_bool | 0.520 | 0.002 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 97815.000 | 1.000 |
prob_mu_under_c | 0.703 | 0.007 | 0.457 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 4861.000 | 1.001 |
prob_es_under_ces | 0.510 | 0.007 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 5481.000 | 1.002 |
prob_prob_under_c_bool_over_percent | 0.520 | 0.002 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 97815.000 | 1.000 |
Rhat値が全パラメータで1.1以下になっていることや、Traceplotの様子からMCMCが収束していることが確認できると思います。 というわけで推定値は一応問題がありません。実際に推定値を用いて推論を行うことが出来ます。
仮説と検証
平均値
- このデータにおける平均値はいくつか
- 事後分布のEAPが40.660(0.650)であるため、40.660と推定される
- このデータにおいて、平均値が95%の確率で分布する区間はどこか
- このデータにおいて平均値が片側で95%の信頼区間を持つ区間はどこか
- 事後分布の平均値の95%信頼区間より、41.735以下になる確率が95%である
- もしくは39.603以上になる確率が95%である
標準偏差
予測分布
- 新しい人からデータをとったときに、95%の確率で何点を取るでしょうか
- 27.769から53.520点の間のいずれかの点数をとる
生成量
- 分散の点推定と95%信頼区間推定
- 42.772(6.178)[32.033, 56.063]
- 基準値を41とした際の、変動係数の点推定と95%信頼区間推定
- 0.160(0.012)[0.139, 0.185]
- 基準値を41とした際の、効果量の点推定
- -0.524(0.099)
- 基準値を41とした際の、効果量の95%区間推定(両側/片側)
- 両側[-0.248, 0.144]、片側[-0.215, ] or [, 0.112]
- 70%点の点推定と区間推定
- 45.605(0.744)[44.204, 47.132]
- 41以下の値が観測される確率と95%信頼区間
- 分布関数を利用 : 0.521(0.039)[0.443, 0.598]
- 確率的命題を利用 : 0520
- 41と測定値の比の、点推定値と95%信頼区間推定値
- 0.993
研究仮説が正しい確率
- 平均値が41より小さい確率
- 0.703
- 基準値を41とした際の効果量が、-0.05より小さい確率
- 0.510
- 41未満の測定値が観測される確率が、70%より大きい確率
- 0.512
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その2
前置き
環境構築
まずPythonを実行する環境を作ります。 一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。 以下のものはOSXを想定しています。
自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。 導入は各リポジトリのページに書いてあるとおりです。
pythonは3系を使います。今回利用するライブラリはそんなに変なものではないので、2系である必要はありません。 解析はスクリプトでやるよりも、jupyterを使う方が僕は好きです。Stan言語をいじれない問題点はありますが、パラメータを変えるのが楽です。
pyenv install 3.6.1 pyenv virtualenv 3.6.1 3.6.1_postp mkdir -p ~/Desktop/postp # 作業ディレクトリは適当に変えてください。 cd ~/Desktop postp pyenv local 3.6.1_postp pip install numpy scipy pandas matplotlib cython seaborn jupyter pip install pystan
1章 データの整理とベイズの定理
内容
一応箇条書きで書いていますが、メモです。 ちゃんと勉強したいときは、本文を参考にしてくださいね。
- 分布の紹介
- 正規分布
- 一様分布
- データの解釈
章末問題
- pythonにて解析をしていきます。
データの準備
僕の好みですが、printではなくloggingモジュールを利用して途中経過を書いていきます。計算時間なども把握できるのが便利です。
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 # データの準備 x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33]) logger = get_logger()
度数分布表
いきなりですが、pandasでは階級幅を決めるメソッドが無いので超めんどくさい… 手動でデータを入れるbinsを作成した後、pd.cutを使って仕分けます。
import pandas as pd df = pd.DataFrame(x, columns=["weight"]) width = 5 minim_range = int((df.min() / width).apply(np.floor) * width) maxim_range = int((df.max() / width).apply(np.ceil) * width) bins = range(minim_range, maxim_range + 1, width) c = pd.cut(df["weight"], bins) d_df = df.groupby(c).count() d_df.columns = ["count"] d_df.columns.name = "range" d_df.index.name = None d_df
ヒストグラムの作成
seabornを使って可視化します。matplotlibではないのは、やはり好みです pandas.plotと違ってseaborn / matplotlibのhistメソッドはbinsに配列を受け取ってくれます。
import seaborn as sns import pandas as pd %matplotlib inline df = pd.DataFrame(x, columns=["weight"]) width = 5 minim_range = int((df.min() / width).apply(np.floor) * width) maxim_range = int((df.max() / width).apply(np.ceil) * width) bins = range(minim_range, maxim_range + 1, width) sns.distplot(df["weight"], kde=False, bins=bins)
統計量
最頻値だけはメソッドが無いので、countをとって最大の値を持つものを返しています。 pandasのvar, stdメソッドはデフォルトパラメータだと母分散、母標準偏差を返す事に注意してください。 numpyは標本分散、標本標準偏差を返す。謎の挙動。
import pandas as pd df = pd.DataFrame(x, columns=["weight"]) # 標本平均 df_mean = df["weight"].mean() logger.info("平均値は{0}".format(df_mean)) # 標本分散 df_var = df["weight"].var(ddof=0) logger.info("標本分散は{0}".format(df_var)) # 標本標準偏差 df_std = df["weight"].std(ddof=1) logger.info("標本標準偏差は{0}".format(df_std)) # データのソート df["weight"].sort_values() # 最大値 df_max = df["weight"].max() logger.info("最大値は{0}".format(df_max)) # 最小値 df_min = df["weight"].min() logger.info("最小値は{0}".format(df_min)) # 中央値 df_medium = df["weight"].median() logger.info("中央値は{0}".format(df_min)) # 最頻値 df_most_freq = df["weight"].value_counts().iloc[0] logger.info("最頻値は{0}".format(df_most_freq))
正規分部の扱い
scipy.statsのnormオブジェクトを利用します。 片側だけから特定の確率になる値を求めるメソッドが無いですが、intervalメソッドで両側の場合の値を求めて、2で割って考えれば良いです。 これは正規分布の対称性を利用したものですね。
import pandas as pd from scipy.stats import norm df = pd.DataFrame(x, columns=["weight"]) # パラメータを設定 df_mean = df["weight"].mean() df_std = df["weight"].std(ddof=1) logger.info("正規分布のパラメータ: 平均は{0}、標準偏差は{1}".format(df_mean, df_std)) # 分布の平均が40.64なので、40付近の値が観測されやすい # 45以上の値が観測される確率 # 1から45までの値が観測される確率を引く(1-cdf) # scipyでは、sfを使えば良い from scipy.stats import norm prob_upper_45 = norm.sf([45], loc=df_mean, scale=df_std)[0] logger.info("45以上の値が観測される確率は{0}".format(prob_upper_45)) # 35以上45未満の値が観測される確率 # 45以下が観測される確率から、35以下が観測される確率を引く prob_upper_35 = norm.cdf([35], loc=df_mean, scale=df_std)[0] prob_upper_45 = norm.cdf([45], loc=df_mean, scale=df_std)[0] prob_between_35_45 = prob_upper_45 - prob_upper_35 logger.info("35以上45未満の値が観測される確率は{0}".format(prob_between_35_45)) # 95%信頼区間 # 正規分布において、平均 ± 1.96 * 標準偏差の範囲が該当することを利用する #min_edge = df_mean - df_std * 1.96 #max_edge = df_mean + df_std * 1.96 #prob_min = norm.cdf([min_edge], loc=df_mean, scale=df_std)[0] #prob_max = norm.cdf([max_edge], loc=df_mean, scale=df_std)[0] # prob_sig_range = prob_max - prob_min # が、scipyだとintervalを利用すれば良い min_edge, max_edge = norm.interval(0.95, loc=df_mean, scale=df_std) logger.info("95%信頼区間は{0}から{1}".format(min_edge, max_edge)) # p(x>a) = 0.05であるような点a # F(平均 + 1.64 * 標準偏差)が大体0.95になる #prob_less_095 = norm.sf([df_mean + df_std * 1.64], loc=df_mean, scale=df_std)[0] # 何も情報が無いとしたら、正規分布性を元に、90%信頼区間を出して、上端を利用すれば良い _, less_095 = norm.interval(0.9, loc=df_mean, scale=df_std) prob_less_095 = norm.cdf([less_095], loc=df_mean, scale=df_std)[0] logger.info("a = {0} であれば、p(x>a)となる確率は{1}である".format(less_095, prob_less_095)) # 3つの4分位点 # 50%信頼区間と中央値を出せば良い quad_first, quad_third = norm.interval(0.5, loc=df_mean, scale=df_std) quad_second = df_medium logger.info("第1四分位点は{0}, 第2四分位点は{1}, 第3四分位点は{2}".format(quad_first, quad_second, quad_third))
Rで出来ることがPythonだとめんどくさかったりしますね。 特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります!