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