Easy to type

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

トップ > ベイジアンモデリング > はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―　その7

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

ベイジアンモデリング

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

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

内容

  • 離散的な値をとるカウントデータの解析

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

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

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

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

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

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

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

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

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

    • 2 × 2のクロス表: 2つの事象の発生にどれだけ関係があるか？
      • 2つの分布が独立ではなく、対応がある場合
      • あるいはベルヌイ試行の興味が2つの分布の同時性にある場合
      • f(x|\theta) = f(x_{11}, x_{12}, x_{21}, x_{22}| p_{11}, p_{12}, p_{21}, p_{22}) ; fは多項分布
      • 独立と連関
        • 変数AのカテゴリA_i, 変数BのカテゴリB_jが独立である場合、f(B_j) = f(B_j|A_i)
        • ピアソン残差 : e_{ij} = \frac{p_{ij} - p_{i.}p_{.j}}{\sqrt{p_{i.}p_{.j}}}
          • セルごとに計算され、独立性/連関性を示す
          • 独立である場合0となる
          • 正の値は独立の場合より観測されやすく、負の場合は独立の場合より観測されづらいことを示す
        • クラメルの連関係数 : V = \sqrt{\frac{\Sigma_{i=1}^n \Sigma_{j=1}^m e_{ij}^2}{min(a, b) - 1}};n, mはカテゴリの総数
          • クロス表全体での連関の大きさを示す
          • 0 \lt V \lt 1
          • 小さいほど独立、大きいほど連関である

章末問題

前準備

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

f:id:ajhjhaf:20170607005421p:plain

比率の差の推測

相談相手問題(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

f:id:ajhjhaf:20170607005519p:plain

治療法問題(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

f:id:ajhjhaf:20170607005541p:plain

連関の推測

広告効果(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

f:id:ajhjhaf:20170607005555p:plain

ワイン(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

f:id:ajhjhaf:20170607005605p:plain

備考

  • 個人的には、クロス集計はよく使いますのでここで扱ってくれたのはいいものだと思いました。

    • ただクロス集計では次のような問題には対応できない
      • A群とB群の人間が居て、それぞれa人とb人存在する
      • 各人は好きな色(全部でn種類)の飴を好きなだけピックアップする
      • 各人の飴の本数、m本は任意の値をとる
      • この時、A群とB群で特定の色の飴が取得されやすいか、されづらいかはどう判定すればよい？
    • こういう問題ですと、より階層的なモデルが必要になってきます。ぱっと思いつくのはディリクレ分布を使って、トピックの混合率をモデル化するトピックモデルを使う手法でしょうか。

  • 「はじめての統計データ分析」はこれで最終回ですね。大分Stanへの恐怖心もなくなってきましたので、次は「実践ベイズモデリング」とか「StanとRでベイズ統計モデリング」あたりを流していこうかと思います。