はじめての統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―　その6

その6です。今回は第5章の章末問題に取り組んでいきます。

5章実験計画による多群の差の推測

内容

分散分析によるF検定へのオルタナ
要因から生じる影響の調査
- 要因 : 離散的なカテゴリーによる影響のこと
- 水準: 要因が取りうる様々な状態のこと
- 水準数: 水準のパターン数。以下要因 $A$ の水準数を $a$ で表す。
独立した1要因モデル
- モデル
  - $y_{ij} = \mu_j + e$
  - $e \sim Normal(0, \sigma_e)$
  - 要因の $j$ 番目の水準における $i$ 番目の測定値が、 $j$ 番目の水準がもつ平均的な効果と正規分布から生じる乱数による現れる
- 以上を纏めて、事後確率は $f(y_{ij} | \mu_j, \sigma_e)$ と表される
- 水準内の測定値は独立に観測されてると考える　
  - 同時確率は $f(y_{ij} | \mu_j, \sigma_e) = \Pi_{i = 1}^n f(y_{ij}|\mu_j, \sigma_e)$
  - 即ち、尤度は $f(y|\theta) =\Pi_{i = 1}^a f(y_j|\mu_j, \sigma_e)$
生成量
- 全平均 :
  - 各水準の平均の単純な平均
- 水準の効果(eol):
  - 各平均から全平均を引いたもの
  - 水準の効果の和が0になることに注意
- 説明率:
  - 測定値の分散に占める要因の分散の割合で、要因が観測値をどれほど説明するかを示す
  - は効果と誤差が互いに独立ならば、である
    - 測定値の分散は要因の分散 + 誤差の分散になるということ
  - $0 \leq \eta \leq 1$ で定義される
- 効果量 :
  - 水準間の標準偏差が、水準内の標準偏差の何倍に相当するか示すもの
  - $\sigma_e$ は要因の影響を受けない、水準内の平均的バラつきと等価であると考える事ができる
モデルを使った要因の効果の方法例
- 水準の効果はどのぐらいか、有るのか無いのか
  - 水準の効果が基準値cより大きい確率をEAPで評価する
- その要素はどれぐらいの効果の大きさを持つのか
  - 効果量、説明量を利用
- 水準間で比較をするとどうなるのか
  - 水準ごとの平均の差のEAPを調べる
- 連言命題が正しい確率はどれぐらいか
  - 連言命題 : 複数の演算子により表される複数の命題
  - 各命題が成り立つ確率の同時確率のEAPを調べる
- 興味のある2水準の推測
  - 基本的に、今までやってきた2群比較と同様
  - しかし2水準以外の水準から得られるデータも利用するため、全体のデータを使って事後分布を推測する。この方が全体で誤差値が共通であるため、正確な予測が可能となる。
独立した2要因モデル
- 要因を複数導入
  - 水準の組み合わせ1単位のことをセルと以下では呼称する
- 交互作用: 一方の要因の水準毎に、他方の要因の平均の高低パターンが異なる現象
- モデル
  - $y_{ijk} = \mu + a_{j} + b_{k} + ab_{jk} + e$
  - $e \sim Normal(0, \sigma_e)$
  - 水準 $j, k$ の $i$ 番目の測定値が、全体の平均に要因 $a$ の $j$ 番目の水準の効果 + 要因 $b$ の $k$ 番目の水準の効果 + $a$ と $b$ の交互作用の効果 + 乱数から構成されることを示す
  - セルの平均 + 乱数
- それぞれ全平均(水準の効果)の性質から次の制約を持つ(水準の効果を全水準で足すと0になる性質を定式化)
  - $\Sigma_{i=1}^{a} a_{i} = 0$
  - $\Sigma_{i=1}^{b} b_{i} = 0$
  - $\Sigma_{j=1}^{a} ab_{jk} = 0$
  - $\Sigma_{k=1}^{a} ab_{jk} = 0$
- 尤度 $f(y|\theta) = f(y|\mu, \sigma_e = \Pi_{j=1}^{a}\Pi_{k=1}^{b} f(y_{jk} | \mu_{jk}, \sigma_e)$

章末問題

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

f:id:ajhjhaf:20170531233354p:plain

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
changeup	0	0.99474	0	1e-05	0.37879
curve	0.00526	0	0	0	0.00256
cutball	1	1	0	1	1
fork	0.99999	1	0	0	1
slider	0.62121	0.99744	0	0	0
straight	1	1	1	1	1

f:id:ajhjhaf:20170531233333p:plain

考察

水準の効果の分布より、走者の状況も交絡因子も殆ど影響力が無いことが推し量れる
- 説明量や効果量についても、圧倒的に走者(b)についての値が大きく、モデルの殆どがこれで説明される事がわかる
- 交絡因子の大きさについても図れるのはメリットぽい。

Easy to type

個人的な勉強の記録です。データ分析、可視化などをメイントピックとしています。

はじめての統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―　その6

5章実験計画による多群の差の推測

内容

章末問題

1要因の推測

データの取得

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

2要因の推測

データの取得

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

考察

5章 実験計画による多群の差の推測

内容

章末問題

1要因の推測

データの取得

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

2要因の推測

データの取得

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

考察

5章実験計画による多群の差の推測