StanとPythonでベイズ統計モデリング その1

StanとRでベイズ統計モデリング(通称アヒル本)をだいたい読みました。

StanとRでベイズ統計モデリング (Wonderful R)

本の紹介

既に様々な書評もありますし、方々から賛辞の声を挙げられている本です。僕としても非常に分かりやすく、使える本だと感じました。著者の松浦さんがウリを書いてくださっているので、まずそれを読むのが良いと思います。様々な方の書評も纏められています。

statmodeling.hatenablog.com

読んでみての感想を、良い点と改訂版に期待する点(笑)で書いてみたいと思います。

良い点

使えそうな分布が結構紹介されている

これは僕の専門分野が生命情報解析、著者の松浦さんの専門が医療統計で近いということもあるのですが、使えそうでありつつ統計入門レベルでは出てこない分布が実例を交えて紹介されています。例えば

• コーシー分布
• ゼロ過剰ポアソン分布
• ディリクレ分布

などです。ここらへんは数式だけではよく分かりませんし、とっつきづらいなところです。Stanと交えながら活用例を紹介してくださっているのは、イメージがつかみやすかったです。

階層モデル以外のベイジアンモデリングがふんだんに登場する

ベイジアンモデリングに挑戦してみよう、となると久保先生の通称緑本が第一の選択肢にあがります。僕も1年前に読みました。緑本はベイズ化する前の線形モデルから始まり、ベイズまで展開していくことで個人差が生じる過程をどう解析するのか非常に分かりやすく理解することが出来ます。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

ただ、場所差、個体差などのカテゴリーの違いで本が終わってしまうので「ベイズ化すればパラメータ推定が楽になるんだな」ぐらいの感覚におちついてしまいました。これは緑本が「ベイジアンモデル」の本ではなく、「モデリング」の本ということです。

アヒル本はベイズモデルに特化しており、階層モデル以外にも

• 状態空間モデル
• トピックモデル
• 打ち切りや外れ値への対応

が紹介されています。そのためベイジアンモデルの「あいまいさ自体をパラメータ化してしまう」という特徴が他にどのように活かせるのか理解しやすく、応用範囲の広さにも気が付くことが出来ました。

ベイズモデルの難しいところは、「パラメータを生み出す分布がある」という感覚の理解かと思っています。トピックモデルなんかは、これがないと数式がごっついので難解ですし。アヒル本はその感覚の養成に役立ちます。

有用なMCMCを使ったベイジアンモデルの参考書

Stanを使ったMCMCでのモデルフィッティングにおいて、想定される問題点の解決策が述べられています。

• モデルが収束しない
• そもそもStanがコンパイル出来ない
• 離散値パラメータの処理が出来ない

ここらへんはStanのマニュアルを読めば書いてある事象ではあります。ただ僕は読んでみても行間が抜けすぎていて、わからない部分が結構出てきました。 なんたって、600ページありますから…。

アヒル本はマニュアルの時には行間として省略されていたところを初学者向けに補足してくれているので、緑本を読んだ後ぐらいだったら殆ど詰まること無く読めました。現状だと確率的プログラミング言語としてEdwardがアツい印象がありますが、アヒル本はベイジアンモデルの内容を包括的に書いてくれているのでStanが使えなくなったから無駄になる知識ではなさそうです。

文法が新しい

Stanは2.10.0辺りで新しい文法が提案されたようです。そのためマニュアルを読んでも、本中で紹介されているtarget記法が出てこなかったりします。これは2.9.0を参考に和訳されている日本語マニュアルでも同様です。

アヒル本の文法は、2017-06-14時点では問題なく新しいバージョンで書かれています。なのでエラーメッセージにうんざりすることもありません。

改訂版に期待する点

紹介した分布を全部使って欲しい

6章で分布を幾つか紹介しています。先述したコーシー分布やディリクレ分布などはここで登場します。

ただ、ここでラプラス分布は紹介されるものの、本中ではそれっきりです。 ぜひこれらの分布についても実例を挙げて解説していただき、活かせる点を見せてくれたらいいなと思いました。きっと何かしらの有用性があると思いますので。

練習問題を面白くして欲しい

4章以降は各章の最後に練習問題が載っています。学習した内容を復習・確認することが出来ますし、答えについてもレポジトリに R + Stanの実装で配布されています。最高です。

github.com

ただあくまで確認が目的なのか、他の参考書のデータを引用していたり、本文中の内容の焼き直しだったりしたので退屈に感じました。練習問題に+alphaの内容を盛り込んで、更なる学習が出来ればよりモチベーションアップに繋がるのではないでしょうか。

---

書評はこんな感じです。 「はじめての統計データ分析」みたいに全部の問題をやったりはしませんが、Python + PyStanで幾つか解いてみたいと思います。タイトルはそういう意味です笑。

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

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

6章 比率とクロス表の推測

内容

• 離散的な値をとるカウントデータの解析
  • 比率の差のz検定、クロス表のカイ二乗検定のオルタナ

• カテゴリカル分布 : カテゴリカルな分類においてカウントとして得られるデータに対応する分布

  • ベルヌイ試行
  • 2項分布
  • 多項分布
    • 各事象が起こる確率が独立であることに注意

• 1つの2項分布の比率の推測:「きのこの山とたけのこの山の好きな人の比率」など

  • 確率の事前分布は無情報であるので、を仮定してMCMCを行う
  • (二項分布を仮定した時の)事象が生起する確率を推測することが可能
  • 生成量
    • オッズ :
      • 2つの選択肢のどちらの確率が大きいかを測ることが出来る
    • 仮説(基準とした確率値に対して、pが下回る確率など; 即ち)の検証

• 1つの多項分布の比率の推測

  • 同様に事前分布は無情報であるので、番目のカテゴリの事前分布にを仮定してMCMCを行う
    • 全ての確率の和が1となるので、カテゴリが全部でn個なら n個目は として決定される制限がある
  • 0からn番目の事象が独立してそれぞれ起こるのであれば、母数が発生する確率は1つ1つの事象をとして である
  • 生成量
    • カテゴリ間の比較( など)
    • 連言命題が正しい確率

• 複数の2項分布の推測(独立したクロス表) : カテゴリが個あって、それぞれで2つの観測結果が有った際にカテゴリ毎に or 全体でどのような差があるかを調べたい場合

  • (独立した)2 × 2のクロス表
    • 興味があるベルヌイ試行の確率だけ見る
    • 男と女でベルヌイ試行の確率が異なる(独立した2項分布が2つある)と考える
    • 
    • 生成量
      • 比率の差: , 性質の比:
        • 性質の違いを考えることが出来る
      • オッズ比 :
        • 正反応は他方の反応の何倍生じやすいか、を表す
  • g × 2 のクロス表
    • 
    • 各f(theta)が得られる確率については何もわからないので、の無情報事前分布を採用

• 対応あるクロス表の推測(構造のある多項分布): 1つの標本から2回測定が有った場合

  • 2 × 2のクロス表: 2つの事象の発生にどれだけ関係があるか？
    • 2つの分布が独立ではなく、対応がある場合
    • あるいはベルヌイ試行の興味が2つの分布の同時性にある場合
    • ;
    • 独立と連関
      • 変数Aのカテゴリ, 変数Bのカテゴリが独立である場合、
      • ピアソン残差 :
        • セルごとに計算され、独立性/連関性を示す
        • 独立である場合0となる
        • 正の値は独立の場合より観測されやすく、負の場合は独立の場合より観測されづらいことを示す
      • クラメルの連関係数 : ;はカテゴリの総数
        • クロス表全体での連関の大きさを示す
        • 
        • 小さいほど独立、大きいほど連関である

章末問題

前準備

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要因モデル

  • モデル
    • 
    • 
    • 要因の番目の水準における番目の測定値が、番目の水準がもつ平均的な効果と正規分布から生じる乱数による現れる
  • 以上を纏めて、事後確率はと表される
  • 水準内の測定値は独立に観測されてると考える　
    • 同時確率は
    • 即ち、尤度は

• 生成量

  • 全平均 :
    • 各水準の平均の単純な平均
  • 水準の効果(eol):
    • 各平均から全平均を引いたもの
    • 水準の効果の和が0になることに注意
  • 説明率:
    • 測定値の分散に占める要因の分散の割合で、要因が観測値をどれほど説明するかを示す
    • は効果と誤差が互いに独立ならば、である
      • 測定値の分散は要因の分散 + 誤差の分散になるということ
    • で定義される
  • 効果量 :
    • 水準間の標準偏差が、水準内の標準偏差の何倍に相当するか示すもの
    • は要因の影響を受けない、水準内の平均的バラつきと等価であると考える事ができる

• モデルを使った要因の効果の方法例

  • 水準の効果はどのぐらいか、有るのか無いのか
    • 水準の効果が基準値cより大きい確率をEAPで評価する
  • その要素はどれぐらいの効果の大きさを持つのか
    • 効果量、説明量を利用
  • 水準間で比較をするとどうなるのか
    • 水準ごとの平均の差のEAPを調べる
  • 連言命題が正しい確率はどれぐらいか
    • 連言命題 : 複数の演算子により表される複数の命題
    • 各命題が成り立つ確率の同時確率のEAPを調べる
  • 興味のある2水準の推測
    • 基本的に、今までやってきた2群比較と同様
    • しかし2水準以外の水準から得られるデータも利用するため、全体のデータを使って事後分布を推測する。この方が全体で誤差値が共通であるため、正確な予測が可能となる。

• 独立した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です
  • 対応が無い場合に使っていた生成量ももちろん使えます(省略)
  • 対応のある値の差(pair diff)
    • のこと
      • 単に群同士を比較するのではなく、群間である対応関係を利用する
    • に従う
      • この導出についてはこれを参考に。
  • 対応ある値の差の効果量(pair es)
    • 
      • のこと
    • 差の平均値に対するばらつきの平均値を評価できる
  • 対応ある値の差の優越率(pair pod)
    • 
      • MCMCでは、を評価
      • なおここでは平均0、標準偏差1とした時の累積分布関数としています
    • 個人の変動が0より大きい確率を評価できる
  • 対応ある値の差の閾上率(pair pbt)
    • 
      • MCMCでは、を評価
    • 個人の変動がより大きい確率を評価できる
  • 同順率(pbc)
    • 
    • 2つのを選んだ時に、2変数の大小が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()

# データの準備
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群でモデルを作る時に、標準偏差は共通させるか？分けるか？
  • 共通させる場合
    • 
    • 
  • 独立に定義する場合
    • 
    • 
  • かならずしも独立に定義したほうが良いわけではない。モデルの複雑性は増してしまうので。
• 生成量の例
  • 平均の差
  • 平均の差の効果量(estimated size)
    • 
  • 非重複度(Cohen’s U) : u1が群2で何%点に相当するか
    • 
    • 0.5 で2群が完全に重複していることを示す
  • 優越率(pid) : 無作為に選んだ1方の群の測定値が、無作為に選んだもう一方の測定値を上回る確率
    • 
      • (zは標準正規分布に従う変数)
    • 評価をするには、標準正規分布におけるの分布関数を調べれば良い
  • 閾上率(pic) : 無作為に選んだ一方の群の測定値と、もういっぽうの測定値の差がcより大きくなる確率
    • 
    • 評価をするには、標準正規分布におけるの分布関数を調べればよい
• モデルの選択
  • 標準偏差は共通させたほうがよい？独立のほうがよい？
    • 基本的に独立の方が制約がなく、子細に利用可能
    • しかし、制約が正しい場合には、安定しやすい
  • WAICの導入
    • 真の分布が得られたときの、未知データに対する平均対数尤度の推定値(豊田実践ベイズ本より)
      • 未知データに対する予測力
      • 小さい方が良い
    • 
      • 詳細はhttp://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html
      • 実践 ベイズモデリング -解析技法と認知モデル-も詳しかったです。
    • WAICはMCMCのChain中にある全サンプル値を使って未知データの予測を行い、過学習している分を補正する感じ？

章末問題

データの取得

とりあえずデータを読み込みます。 今回は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%信頼区間推定 両側/片側
  • 3.308(0.507)[2.365, 4.354]となりました。
  • 片側区間なら、es > 2.508、あるいは es < 4.172の区間に95%の確率で入りました。
• A群から見たB群の効果量が3より大きい確率は？
  • prob_es_upper_cより、72.1%の確率で、3より大きいことがわかります。
• A群の分布とB群の分布の重複する割合(非重複度)の点推定、95%信頼区間推定 両側/片側
  • 0.908(0.034)[0.832, 0.962]となりました。
  • 片側区間なら、cohenu > 2.508、あるいはcohenu < 0.956の区間に95%の確率で入ります。
  • これはB群から見た時に、典型的なA群の測定値は40.8%情報に存在することを示す
• 非重複度が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) : 乱数の内、いくつが無関係な乱数として扱うことが出来るかを示す指標
      • (以上２つの指標の導出が気になった場合、基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門に結構丁寧に書いてあるので、そちらをお薦めします。)

• 事後分布の解釈をどう行うか？

  • 点推定量
    • 事後期待値(EAP) = 平均値として扱うことが出来る
    • 事後中央値(MED) = 中央値として(ry
    • 事後確率最大値(MAP) = 最頻値(ry
    • 事後分散
    • 事後標準偏差
      • 与えられたデータを元に推定した分布自体の広がりがどれぐらいあるかを示したもの
    • 事後期待値の標準偏差(標準誤差; se)
      • EAPがどれぐらいの広がりを持つか示したもの
  • 確信区間(信用区間)
    • 母数が存在する確率がx % である区間を確信区間とする
  • 予測分布
    • 新たなデータが得られうる分布のこと
    • 事後予測分布
      • 
      • 推定された母数の事後分布を用いて、未知データを予測する統計モデルの平均を計算
      • 
      • 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
    • 条件付き予測分布
      • 事後分布の特定の点推定量を利用する
      • 簡単だが、せっかく推定したモデルを有効活用出来ない
      • (これについては、その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%の確率で分布する区間はどこか
  • 事後分布の平均値の2.5%区間と97.5%信頼区間は39.383から41.947である
• このデータにおいて平均値が片側で95%の信頼区間を持つ区間はどこか
  • 事後分布の平均値の95%信頼区間より、41.735以下になる確率が95%である
  • もしくは39.603以上になる確率が95%である

標準偏差

• 標準偏差の平均値、両側/片側95%信頼区間はどこか
  • 6.523(0.468)、5.660から7.488、5.797以上あるいは7.321以下である

予測分布

• 新しい人からデータをとったときに、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