Easy to type

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

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

アヒル本(StanとRでベイズ統計モデリング)のChapter5にPythonで取り組んでいきます。 練習問題を解いて、本文中に書かれてるグラフをPythonで描いてみます。

なおChapter1~3は導入だったのと、Chapter4は練習問題の内容が「はじめての統計データ分析」と被っていたのでパスします。


Chapter5 基礎的な回帰とモデルのチェック

  • 重回帰
    • 複数の説明変数を用いた回帰のこと
    • 目的
      • 説明変数からの応答変数の予想、及び説明変数の寄与率
    • 分布
      • 複数の説明変数ならScatterplot matrixを利用すると良い
    • MCMCの設定について
      • スケーリング: MCMCを行う前に、各データのオーダーを大体(1前後に)そろえること。
    • 結果の解釈 : モデルの改善に活かす!
      • 因果関係の推定
        • どちらが先に起こりうるか、といった事前知識が必要
      • 外挿
        • 利用したデータが  a \lt x \lt bの範囲であるにも関わらず、 c \lt\lt aのようなデータ cをモデルに当てはめて予測することを外挿と呼ぶ。大体は予想外のデータなので、上手く行かない。
      • パラメータの幅
        • EAPなどの点推定量より、幅がどうなっているかが大事(0を跨いだりしていないか?)
      • 可視化
        • 実測値と予測値のプロット(データと予測分布)
          • X軸にデータ、Y軸にEAPから得られる予測値をプロット
          • y = xの直線から大きく外れることがなければOKとする
        • 推定ノイズの分布
          • ノイズの分布を仮定した時に、それがモデルに従っているか確認する
      • MCMCサンプルの散布図行列
        • パラメータの関係について読み取れる
  • (二項)ロジスティック回帰
    • ロジスティック関数\frac{1}{1+\exp(-x)}をつかった方法
      • 応答変数が0から1の間にある場合、説明変数に適用することで値の範囲を変換することが出来る
    • 2項分布モデル
      • 正規分布が当てはまらない場合には特に有効
      • 値が整数で幅が小さく、連続値とみなせないような場合に使う
      • Nが十分に大きい場合は正規分布に近似できる(正規分布の方が演算が高速なのでよい)
    • ベルヌーイ分布モデル
      • 応答変数がbool値であるような時にはベルヌーイ分布を利用すれば良い
    • 生成量
      • オッズ
        • 二項ロジスティック回帰モデルにおいて、Y \sim Binomial(M, q)であり、q = \frac{1}{1+\exp(-x)}であるならば\frac{q}{1-q} = exp(x)と変換することが出来る
        • \frac{q}{1-q}はオッズと呼ばれる
        • その事象が起こるサンプリング確率の比率を示す
    • 可視化
      • ベルヌーイ分布を利用した場合
        • Y \sim Bernoulli(q) より推定されたqを利用する
        • 判断基準の閾値cを決定して、qとcからTP, FP, FN, TNの計算を行う
        • FP, TPを各qについてplotして曲線(ROC)を書く
        • ROCを下回った部分がAUCで、0.8以上だと良いものとされる
  • ポアソン回帰
    • カウントデータに対して予測を行う時に用いる
      • Y \sim poisson(\lambda)でカウントデータYが決定されるものとして、\lambdaは説明変数の線形結合を0以上に変換したものを使う(exp(x)がよく使われる)
      • パフォーマンスの理由より、\lambdaが十分に大きいなら、正規分布で近似して良い

練習問題

Scatterplot matrixについて

Pythonの実装はpandasマスターのid:sinhrksさんの記事を見るのが良いと思います。

sinhrks.hatenablog.com

個人的な使った感想は次の通りです。

  • 離散値の可視化が上手く行かないことがある(6離散値指定したのに、4つしか表示されないことがある)?
    • たぶんこれはseabornのPairPlotのバグです……
    • 離散値/連続値判定の初期値が大きすぎると上手くいかないので、離散値の判定閾値を4, 5辺りにいじって使いました。
  • 離散値の場合は対角のカウントはcountplotを使ったほうが良い?
    • 他の下半分、上半分のグラフと軸を合わせやすくするためです。
    • どうせseabornはPairGridで使うことになりますし…

そういうわけで改造させて頂いたバージョンがこちらです。残念ながら1つ目のバグは取れていません。そもそも僕の環境に特有なのかもしれませんし、諦めました。 上三角行列の部分はid:statmodelingさんの記事のまんまです。

statmodeling.hatenablog.com

#!/usr/bin/env python
# vim:fileencoding=utf-8
# Created: 2017-06-06
# Original scripts: http://sinhrks.hatenablog.com/entry/2016/11/01/075527

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib.patches import Ellipse


class Dispatcher(object):

    def __init__(self, fontsize=20, alpha=0.6, cmap='RdBu', threshold=5):
        self.fontsize = fontsize
        self.alpha = alpha
        self.cmap = plt.get_cmap(cmap)

        # 離散値 / 連続値とみなす閾値
        self.threshold = threshold

    def comb(self, x_series, y_series, label=None, color=None):
        """ 下三角部分のプロット """

        x_nunique = x_series.nunique()
        y_nunique = y_series.nunique()

        if x_nunique < self.threshold and y_nunique < self.threshold:
            # 離散値 x 離散値のプロット
            return self._dd_plot(x_series, y_series, label=label, color=color)

        elif x_nunique < self.threshold or y_nunique < self.threshold:
            # 離散値 x 連続値のプロット
            return self._dc_plot(x_series, y_series, label=label, color=color)

        else:
            # 連続値 x 連続値のプロット
            return plt.scatter(x_series, y_series, label=label, color=color)

    def _dd_plot(self, x_series, y_series, label=None, color=None):
        """ 離散値 x 離散値のプロット """

        # x, y 各組み合わせの個数を集計
        total = y_series.groupby([x_series, y_series]).count()

        # x, y軸をプロットする位置を取得
        xloc = total.index.labels[0]
        yloc = total.index.labels[1]
        values = total.values

        ax = plt.gca()
        for xp, yp, vp in zip(xloc, yloc, values):
            ax.annotate(vp,
                        (xp, yp),
                        fontsize=self.fontsize,
                        ha='center',
                        va='center')

        # 組み合わせの個数を散布図としてプロット
        size = values / (values.max() * 1.1) * 100 * 20
        ax.scatter(xloc, yloc, s=size, label=label, color=color)
        ax.set_ylim(yloc[0] - 0.5, yloc[-1] + 0.5)

    def _dc_plot(self, x_series, y_series, label=None, color=None):
        """ 離散値 x 連続値のプロット """

        if y_series.nunique() < x_series.nunique():
            # y軸が離散値の場合は、x, yを入替
            # 水平方向に箱ひげ図をプロット
            orient = "h" 
        else:
            orient = "v" 

        ax = sns.boxplot(x_series, y_series, orient=orient, color=color)
        ax = sns.swarmplot(x_series, y_series, orient=orient, color=color) 

    def diag(self, x_series, label=None, color=None):
        """ 対角部分のプロット """

        x_nunique = x_series.nunique()
        if x_nunique < self.threshold:
            ax = sns.countplot(x_series, color=color)
        else:
            ax = sns.distplot(x_series, kde=False, color=color)
        ax.yaxis.set_visible(False)

    def ellipse(self, x_series, y_series, label=None, color=None):
        """ 上三角部分のプロット """

        # 相関係数を楕円としてプロット
        r = x_series.corr(y_series)
        c = self.cmap(0.5 * (r + 1))
        ax = plt.gca()
        ax.axis('off')
        artist = Ellipse(xy=[.5, .5],
                         width=np.sqrt(1+r),
                         height=np.sqrt(1-r),
                         angle=45,
                         facecolor=c,
                         edgecolor='none',
                         transform=ax.transAxes)
        ax.add_artist(artist)
        ax.text(.5,
                .5,
                '{:.0f}'.format(r*100),
                fontsize=28,
                horizontalalignment='center',
                verticalalignment='center',
                transform=ax.transAxes)

       
d = Dispatcher()

# データを作る
x = np.random.normal(100, 5, 30)
y = np.random.choice([1,2,3], 30, p=[0.1,0.2,0.7])
z = np.random.choice([1,2,3,4], 30, p=[0.2,0.2,0.2,0.4])
w = np.random.normal(4, 0.1, 30)
df = pd.DataFrame(np.array([x,y,z,w]).T, columns=["x", "y", "z", "w"])

# 可視化する
g = sns.PairGrid(df, diag_sharey=False)
g.map_diag(d.diag)
g.map_lower(d.comb)
g.map_upper(d.ellipse)
plt.show()

f:id:ajhjhaf:20170606191215p:plain

これとは別に、勉強として自分でも実装してみましたが、離散×離散の可視化が難しいところです。僕はヒートマップを使おうかと思ったのですが、軸を表示しようとする関係でズレます。つらい。バブルプロットを使うのがテキスト通りの適切なアイディアです。 一応載せておきます。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib.patches import Ellipse
from scipy import stats

def value_type(v):
    if v.unique().size < 5:
        v_type = "discrete"
    else:
        v_type = "continuous"
    return v_type


def correlation(x, y, **kws):
    r, _ = stats.spearmanr(x, y)
    ax = plt.gca()
    ax.axis('off')
    ellcolor = plt.cm.RdBu(0.5*(r+1))
    txtcolor = 'black' if np.absolute(r) < 0.5 else 'white'
    artist = Ellipse(xy=[.5, .5],
                     width=np.sqrt(1+r),
                     height=np.sqrt(1-r),
                     angle=45,
                     facecolor=ellcolor,
                     edgecolor='none',
                     transform=ax.transAxes)
    ax.add_artist(artist)
    ax.text(.5,
            .5,
            '{:.0f}'.format(r*100),
            color=txtcolor,
            fontsize=28,
            horizontalalignment='center',
            verticalalignment='center',
            transform=ax.transAxes)


def diagonal(x, **kws):
    ax = plt.gca()
    x_type = value_type(x)
    if x_type == "discrete":
        sns.countplot(x)
    else:
        sns.distplot(x, kde=False)


def distribution(x, y, **kws):
    ax = plt.gca()
    x_type = value_type(x)
    y_type = value_type(y)
    if x_type == "discrete":
        if y_type == "discrete":
            df = pd.DataFrame([x, y]).T
            df.columns=[x.name, y.name]
            pdf = df.pivot_table(index=x.name, columns=y.name, aggfunc=lambda x: len(x)).fillna(0).T
            plt.pcolor(pdf)
        else:
            sns.boxplot(x, y)
            sns.swarmplot(x, y)
    else:
        if y_type == "discrete":
            sns.boxplot(x, y, orient="h")
            sns.swarmplot(x, y, orient="h")
        else:
            plt.scatter(x, y)


def scatterplot_matrix(df):
    g = sns.PairGrid(df)
    g.map_upper(correlation)
    g.map_diag(diagonal)
    g.map_lower(distribution)
    return g

x = np.random.normal(100, 5, 30)
y = np.random.choice([1,2,3], 30, p=[0.1,0.2,0.7])
z = np.random.choice([1,2,3,4], 30, p=[0.4,0.3,0.2,0.1])
w = np.random.normal(4, 0.1, 30)
df = pd.DataFrame(np.array([x,y,z,w]).T, columns=["x", "y", "z", "w"])
g = scatterplot_matrix(df)
plt.show(g)

f:id:ajhjhaf:20170530195929p:plain

練習問題(1)

まぁほとんど弄る部分もありません。 いつも通りPyStanを回すだけです。

import pystan
import pandas as pd
import os
import pickle
from IPython.core.display import display

stan_file = os.path.join("RStanBook", "chap05", "model", "model5-3.stan")
model = pystan.StanModel(stan_file)
    
df = pd.read_csv("./RStanBook/chap05/input/data-attendance-1.txt")
N = len(df)
A = df["A"].values
Score = df["Score"].values / 200
Y = df["Y"].values
stan_data = {
    "N": N,
    "A": A,
    "Score": Score,
    "Y": Y 
}
fit = model.sampling(
    data=stan_data,
    iter=2100,
    chains=4,
    warmup=100,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

mu = fit.extract()["mu"]
diff = np.apply_along_axis(lambda x: Y - x, axis=1, arr=mu)

モデルはgithubで公開されているものをそのまま使いました。

github.com

練習問題(2)

Stanの内部の生成量を出すところで残差を計算します。 個人的には汚い部分を全部Stanの中に隠せるのでこっちのほうが好きです。

import pystan
import pandas as pd
import os
import pickle
from IPython.core.display import display

stan_file = os.path.join("RStanBook", "chap05", "model", "model5-3_kai.stan")
stan_file_c = os.path.join("RStanBook", "chap05", "model", "model5-3_kai.pkl")

# コンパイルチェック
if os.path.exists(stan_file_c) == False:
    # コンパイル済みのバイナリが無かったらstanファイルからコンパイル
    model = pystan.StanModel(stan_file)    
    with open(stan_file_c, "wb") as f:
        pickle.dump(model, f)
else:
    # バイナリがあったらそれを利用
    with open(stan_file_c, "rb") as f:
        model = pickle.load(f)
    
df = pd.read_csv("./RStanBook/chap05/input/data-attendance-1.txt")
N = len(df)
A = df["A"].values
Score = df["Score"].values / 200
Y = df["Y"].values
stan_data = {
    "N": N,
    "A": A,
    "Score": Score,
    "Y": Y 
}
fit = model.sampling(
    data=stan_data,
    iter=2100,
    chains=4,
    warmup=100,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

モデルは次の通りです。

data {
  int N;
  int<lower=0, upper=1> A[N];
  real<lower=0, upper=1> Score[N];
  real<lower=0, upper=1> Y[N];
}

parameters {
  real b1;
  real b2;
  real b3;
  real<lower=0> sigma;
}

transformed parameters {
  real mu[N];
  for (n in 1:N)
    mu[n] = b1 + b2*A[n] + b3*Score[n];
}

model {
  for (n in 1:N)
    Y[n] ~ normal(mu[n], sigma);
}

generated quantities {
  real y_pred[N];
  real diff[N] ;
  for (n in 1:N){
    y_pred[n] = normal_rng(mu[n], sigma) ;
    diff[n] = Y[n] - mu[n] ;
  }
}

つづいてpythonでグラフを描いて、目視で比較をします。

まず実測値とモデルによる予測値を比較してみます。 テキストでは予測分布を使っていますが、試しにMCMCで推定されたパラメータ分布の中央値を使って値を予測してみます。この目的は、後述する事後予測分布との比較です。これは「はじめての統計データ分析」で言うところの、条件付き予測分布を使った手法です。こちらの手法だとStanファイル中で予測をしなくていい代わりに、エラーバーをつけるのがめんどくさいです。何故ならYを導出するための\muにも \sigmaにも分布があるので、その両方を考慮する必要があります。今回はパスさせて下さい。

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# Aの値でindexを分類する
A0_samples = df[df["A"] == 0].index
A1_samples = df[df["A"] == 1].index

# summaryのtableからmuの中央値とsigmaを取得する
summary_df_A0 = summary_df.loc[["mu[{0}]".format(i) for i in A0_samples], :]
summary_df_A1 = summary_df.loc[["mu[{0}]".format(i) for i in A1_samples], :]
map_A0 = summary_df_A0["50%"]
map_A1 = summary_df_A1["50%"]
sigma = summary_df.loc["sigma", "50%"]

# パラメータを使って条件付き予測分布を推定する
Y_A0 = Y[A0_samples]
Y_A1 = Y[A1_samples]
Y_A0A1 = np.hstack([Y_A0, Y_A1])
Y_ast_A0 = np.random.normal(map_A0, sigma)
Y_ast_A1 = np.random.normal(map_A1, sigma)
Y_ast_A0A1 = np.hstack([Y_ast_A0, Y_ast_A1])
# グラフのために相関係数を導出する
r = np.corrcoef(Y_A0A1, Y_ast_A0A1)[0][1]

# 可視化
fig, ax = plt.subplots()
ax.plot(Y_ast_A0, Y_A0, "o", label="A=0")
ax.plot(Y_ast_A1, Y_A1, "^", label="A=1")
ax.plot([0, 0.5], [0, 0.5], "k-")
ax.legend()

# 体裁を整える
ax.set_aspect('equal')
ax.set_xlabel("Predicted value(by EAP)")
ax.set_ylabel("Observed value")
ax.set_title("r={0}".format(r))
plt.show()

f:id:ajhjhaf:20170615161247p:plain

Aについてバラけているわけではないですが、あまり精度が良く無さそうです。モデルの精度が良くないのでしょうか(後述しますが、これは違います)?

続いて、MCMCサンプリングの中でyの予測値を推定します。これは「はじめての統計データ分析」で言うところの、事後予測分布を使った手法です。こちらのほうが分布全体を考慮して行っているので、正確な予測値を計算することが出来ます。エラーバーについてもMCMCの分布を使えばいいので簡単です。今回は教科書の80%区間ではなく、デフォルトで計算される75%区間を利用しています。

# Aの値でindexを分類する
A0_samples = df[df["A"] == 0].index
A1_samples = df[df["A"] == 1].index

Y_A0 = Y[A0_samples]
Y_A1 = Y[A1_samples]
Y_A0A1 = np.hstack([Y_A0, Y_A1])

# プロット点の計算
summary_df_A0 = summary_df.loc[["y_pred[{0}]".format(i) for i in A0_samples], :]
summary_df_A1 = summary_df.loc[["y_pred[{0}]".format(i) for i in A1_samples], :]
Y_pred_A0 = summary_df_A0["50%"]
Y_pred_A1 = summary_df_A1["50%"]
Y_pred = np.hstack([Y_pred_A0, Y_pred_A1])
# 誤差の計算
Y_pred_A0_err_down = (Y_pred_A0 - Y_pred_summary_A0["25%"]).values
Y_pred_A1_err_down = (Y_pred_A1 - Y_pred_summary_A1["25%"]).values
Y_pred_A0_err_up = (Y_pred_summary_A0["75%"] - Y_pred_A0).values
Y_pred_A1_err_up = (Y_pred_summary_A1["75%"] - Y_pred_A1).values
Y_pred_A0_err = np.vstack([Y_pred_A0_err_down, Y_pred_A0_err_up])
Y_pred_A1_err = np.vstack([Y_pred_A1_err_down, Y_pred_A1_err_up])
# 相関係数の計算
r = np.corrcoef(Y_A0A1, Y_pred)[0][1]

# 可視化
fig, ax = plt.subplots()
ax.errorbar(Y_pred_A0, Y_A0, yerr=Y_pred_A0_err, fmt="o", label="A=0")
ax.errorbar(Y_pred_A1, Y_A1, yerr=Y_pred_A1_err, fmt="^", label="A=1")
# y=xの線を書く
ax.plot([0, 0.5], [0, 0.5], "k-")
# 凡例をつける
ax.legend()

# 体裁を整える
ax.set_aspect('equal')
ax.set_xlabel("Predicted value")
ax.set_ylabel("Observed value")
ax.set_title("r={0}".format(r))
plt.show()

f:id:ajhjhaf:20170615161515p:plain

結果がわりと異なります。相関係数に如実に現れました。条件付き予測分布ではダメかと思ったモデルですが、こちらでは問題がなさげなことが推定されます。以上のことから、条件付き予測分布はオンラインで利用できるなどの利点がありますが、その分結果の浮動性があることを頭に入れたほうが良いと思います。とはいえ、どちらの方法でもAに関する偏ったバラつきが無いことはわかりますので、何を見たいかの使い分けが結局のところ肝心です。

最後に実測値とパラメータ\muの誤差を調べます。

summary_df_diff = summary_df.loc[["diff[{0}]".format(i) for i in range(0, 50)], :]

fig, ax = plt.subplots()
sns.distplot(summary_df_diff["50%"], bins=np.arange(-0.15, 0.15, 0.03), ax=ax)
ax.set_xlabel("Epsilon EAP")
ax.set_ylabel("Count")
plt.show()

f:id:ajhjhaf:20170615164959p:plain

教科書の結果から分かっていましたが、きちんとあてはめた正規分布に従っています。

練習問題(3)

離散値と離散値の集計です。pandasのpivot_tableメソッドを使えばいいのですが、芸がないのでヒートマップで可視化しておきます。AとYの各組み合わせで足していけばいいので、ラムダ式を読み出してlenでカウントしました。

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.read_csv("./RStanBook/chap05/input/data-attendance-3.txt")
pdf = df[["A", "Y"]].pivot_table(index="A", columns="Y", aggfunc=lambda x:len(x))
g = sns.heatmap(pdf, annot=True)
plt.show(g)

f:id:ajhjhaf:20170615170145p:plain

練習問題(4)

晴れがもたらす影響を基準値として、曇と雨が持つ影響についてもMCMCから推定してみます。これだけだと芸が無いので、WAICによるモデルのフィッティングについても評価してみます。

import pystan
import pandas as pd
import os
import pickle
import numpy as np
from IPython.core.display import display
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()

stan_file = os.path.join("RStanBook", "chap05", "model", "model5-5_kai.stan")
model = pystan.StanModel(stan_file)    
    
trans = {"A": 1, "B": 2, "C": 3}
df = pd.read_csv("./RStanBook/chap05/input/data-attendance-3.txt")
I = len(df)
A = df["A"].values
Score = df["Score"].values / 200
W = df["Weather"].values
W = np.array([trans[w] for w in W])
Y = df["Y"].values
stan_data = {
    "I": I,
    "A": A,
    "Score": Score,
    "W": W,
    "Y": Y 
}
fit = model.sampling(
    data=stan_data,
    iter=1200,
    chains=4,
    warmup=200,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

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))

モデルは以下の通りです。

data {
    int I;
    int<lower=0, upper=1> A[I];
    real<lower=0, upper=1> Score[I];
    int<lower=1, upper=3> W[I];
    int<lower=0, upper=1> Y[I];
}

parameters {
    real b[3];
    real b_w[2] ;
}

transformed parameters {
    real q[I];
    real beta_w[3] ;
    beta_w[1] = 0 ;
    beta_w[2] = b_w[1] ;
    beta_w[3] = b_w[2] ;

    for (i in 1:I)
        q[i] = inv_logit(b[1] + b[2]*A[i] + b[3]*Score[i] + beta_w[W[i]]);
}

model {
    for (i in 1:I)
        Y[i] ~ bernoulli(q[i]);
}

generated quantities {
    real log_lik[I] ;

    for (i in 1:I){
        log_lik[i] = bernoulli_lpmf(Y[i]|q[i]) ;
    }
}

さて経験的にパラメータを決定した時と比べて、WAICはどうなったでしょうか。

  • 経験的に決めた場合 : WAIC=2769
  • MCMCからサンプリングした場合: WAIC=2764

ちょっとだけ、MCMCからサンプリングした場合のほうが値が良くなります。EAPを見てモデルを書いてみると

  • 経験的に決めた場合

 Y \sim Bernoulli(q)

q = \frac{1}{1 + \frac{1}{exp(0.17 - 0.62 * A + 1.95 * Score - 0.45 * Weather)}}

f:id:ajhjhaf:20170606192239p:plain:w200

  • MCMCからサンプリングした場合

Y \sim Bernoulli(q)

q = \frac{1}{1 + \frac{1}{exp(0.26 - 0.62 * A + 1.98 * Score +  Weather)}}

f:id:ajhjhaf:20170606192250p:plain:w200

となります。Weatherの部分も係数の-0.45を考慮すると、そこまで変わりません。僕は経験的に得られた0:0.2:1が結構いいものであったと解釈しました。

つづいてpythonでグラフを描いて、目視で結果を比較します。 テキストに描いてあるように、ロジスティック回帰の観測値は0/1になるので、可視化が難しいところです。qの値が、データからどのように変化するかを調べます。テキストの例では連続値のScoreと離散値のAが可視化されています。

まずScoreの影響を可視化します。

  • 推測したqの中央値をデータと対応付けるために、indexを振り直します。
  • カラムを複数使ったデータの選択なので、queryメソッドで取得します。
  • qの中央値の線をきれいに書くためにソートを挟んでいます
# qの値を取得
summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :]
q_range = summary_df_q[["25%", "50%", "75%"]]
q_range.index = range(len(q_eap))
# データと結合(順番はこの時点では保存されている)
df_concat = pd.concat([df, q_range], axis=1)
# 一部のデータだけ選択
df_filtered = df_concat.query("A == 0 & Weather == 'A'")
df_filtered = df_filtered.sort_values("Score")

# Yの点の位置をそれぞれちょっとずらす
y_randomized = [y + np.random.normal(0, 0.05) for y in df_filtered["Y"]]

# 可視化
fig, ax = plt.subplots()
ax.plot(df_concat_filtered["Score"], df_concat_filtered["50%"], "k-")
ax.fill_between(df_concat_filtered["Score"], df_concat_filtered["25%"], df_concat_filtered["75%"], facecolor="pink")
ax.plot(df_concat_filtered["Score"], y_randomized, "o", label="Y")
ax.legend()
ax.set_xlabel("Score")
ax.set_ylabel("q")
plt.show()

f:id:ajhjhaf:20170615181717p:plain

次にAに対するプロットを書きましょう。matplotlibで頑張るのが面倒だったのでseabornで書きました(なんかごめんなさい)。適材適所です!

# qの値を取得
summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :]
q_range = summary_df_q[["25%", "50%", "75%"]]
q_range.index = range(len(q_eap))
# データと結合(順番はこの時点では保存されている)
df_concat = pd.concat([df, q_range], axis=1)

# 可視化
fig, ax = plt.subplots()
sns.violinplot(data=df_concat, x="50%", y="Y", ax=ax, orient="h")
sns.swarmplot(data=df_concat, x="50%", y="Y", hue="A", ax=ax, orient="h", palette=sns.color_palette("Set2", n_colors=2))
ax.set_xlabel("q")
ax.set_ylabel("Y")
plt.show()

f:id:ajhjhaf:20170615183431p:plain

最後にROCカーブを書いて、AUCを計算します。これについてはScikit-learnを使うのが簡単でしょう。

from sklearn.metrics import roc_curve, roc_auc_score

summary_df_q = summary_df.loc[["q[{0}]".format(i) for i in range(0, I)], :]
q = summary_df_q["50%"]

# AUCの計算
auc = roc_auc_score(Y, q)

# ROCカーブの座標を計算
false_positive_rate, true_positive_rate, _ = roc_curve(Y, q)
fig, ax = plt.subplots()
ax.plot(false_positive_rate, true_positive_rate)
# y=xを描く
ax.plot([0, 1], [0, 1], "k-")

# 体裁を整える
ax.set_aspect('equal')
ax.set_xlabel("False Positive Rate")
ax.set_ylabel("True Positive Rate")
ax.set_title("AUC = {0}".format(auc))
plt.show()

f:id:ajhjhaf:20170615184747p:plain

練習問題(5)

import matplotlib.pyplot as plt
import os
import pandas as pd
import pickle
import pystan
from IPython.core.display import display

stan_file = os.path.join("RStanBook", "chap05", "model", "model5-6b.stan")
model = pystan.StanModel(stan_file)    
    
trans = {"A": 1, "B": 2, "C": 3}
df = pd.read_csv("./RStanBook/chap05/input/data-attendance-2.txt")
N = len(df)
A = df["A"].values
Score = df["Score"].values / 200
M = df["M"].values
stan_data = {
    "N": N,
    "A": A,
    "Score": Score,
    "M": M 
}
fit = model.sampling(
    data=stan_data,
    iter=2100,
    chains=4,
    warmup=100,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

Mpred = np.median(fit.extract()["m_pred"], axis=0)
fig = plt.figure(figsize=(5,5))
ax = fig.gca()
ax.plot(Mpred, M, "o")
ax.plot(range(0, 100), range(0, 100))
ax.set_xlim(0, 100)
ax.set_ylim(0, 100)
ax.set_xlabel("Predicted")
ax.set_ylabel("Measured")

plt.show(ax)

f:id:ajhjhaf:20170615184732p:plain

全然y=xに乗っていません…(合っているのかめちゃめちゃ不安です)。

この結果からは、予測が出来ていないことが分かります。説明変数が足りてないということですね。よくあることですが、計測が終わった後にこうなった場合は、新しくメタデータを足して計測し直しになるわけです。つらい。

練習問題(6)

練習問題(5)と(6)はちょっと気をつける必要があります。アヒル本だけ読んだ時には「処理」の違いというものが

  • 両方に別々の処理を施した
  • 片方には処理を施して、片方には何もしなかった

という2つのパターンのどちらなのか分かりません。前者の場合にはベースラインとなるものに別々の効果が生じることを予測できますが、後者の場合には一方だけはベースラインのままで、もう一方に何か効果が現れることになります。

データの元となった緑本を読んでみると、後者であることが分かります。前者を踏まえてモデルを作ってみると、MCMCが収束せずに頭を悩まるハメになります(僕は緑本を読んだはずなのに忘れていて喘ぎました)。

前者の場合でも相対的な効果を踏まえると同じじゃないの?と思っていたのですが収束しないですし、そうではないんですね。ちょっとした発見でした。これはChapter 10の一番初めに出てくる識別可能性の問題に引っかかるからですね。

import pystan
import pandas as pd
import os
import pickle
from IPython.core.display import display

stan_file = os.path.join("RStanBook", "chap05", "model", "practice6.stan")
model = pystan.StanModel(stan_file)    
    
trans = {"C": 0, "T": 1}
df = pd.read_csv("./RStanBook/chap05/input/data3a.csv")
N = len(df)
X = df["x"].values
F = df["f"].values
F = np.array([trans[x] for x in F])
Y = df["y"].values
stan_data = {
    "N": N,
    "X": X,
    "F": F,
    "Y": Y 
}
fit = model.sampling(
    data=stan_data,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

モデルは次です。

data {
    int N;
    real<lower=0> X[N];
    int<lower=0, upper=1> F[N];
    int<lower=0> Y[N];
}

parameters {
    real b[3] ;
}

transformed parameters {
    real<lower=0> lambda[N] ;

    for (n in 1:N){
        lambda[n] = exp(b[1] + b[2]*X[n] + b[3] * F[n]) ;
    }
}

model {
    for (n in 1:N){
        Y[n] ~ poisson(lambda[n]);
    }
}

練習問題(7)

(6)と同様に処理の効果にだけ気をつけて下さい。

import pystan
import pandas as pd
import os
import pickle
from IPython.core.display import display

stan_file = os.path.join("RStanBook", "chap05", "model", "practice7.stan")
model = pystan.StanModel(stan_file)    
    
trans = {"C": 0, "T": 1}
df = pd.read_csv("./RStanBook/chap05/input/data4a.csv")
I = len(df)
N = df["N"].values
X = df["x"].values
F = df["f"].values
F = np.array([trans[x] for x in F])
Y = df["y"].values
stan_data = {
    "I": I,
    "N": N,
    "X": X,
    "F": F,
    "Y": Y 
}
fit = model.sampling(
    data=stan_data,
    iter=1500,
    chains=4,
    warmup=500,
    seed=1234,
)
summary = fit.summary()
summary_df = pd.DataFrame(
    summary["summary"],
    index =summary["summary_rownames"],
    columns=summary["summary_colnames"]
)
display(summary_df)

モデルは以下の通りです。

data {
    int I;
    int<lower=0> N[I];
    real<lower=0> X[I];
    int<lower=0, upper=1> F[I];
    int<lower=0> Y[I];
}

parameters {
    real b[3] ;
}

transformed parameters {
    real<lower=0> q[I] ;

    for (i in 1:I){
        q[i] = inv_logit(b[1] + b[2]*X[i] + b[3]*F[i]) ;
    }
}

model {
    for (i in 1:I){
        Y[i] ~ binomial(N[i], q[i]) ;
    }
}

練習問題(4), (6)と(7)は初めに取り組んだ時に、モデルが全然収束しなくて足踏みしました。そういう時にはMCMCサンプリングが無理をして変な部分を探索しようとするのでめちゃめちゃ時間がかかります。なので待機時間で「失敗したな」と悟れたりするのはちょっと面白いところです。

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

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

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

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

本の紹介

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

statmodeling.hatenablog.com

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

良い点

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

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

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

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

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

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

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

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

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

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

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

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

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

ここらへんは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章 比率とクロス表の推測

内容

  • 離散的な値をとるカウントデータの解析
  • カテゴリカル分布 : カテゴリカルな分類においてカウントとして得られるデータに対応する分布

    • ベルヌイ試行
    • 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でベイズ統計モデリング」あたりを流していこうかと思います。

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

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

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

内容

  • 分散分析によるF検定へのオルタナ
  • 要因から生じる影響の調査
    • 要因 : 離散的なカテゴリーによる影響のこと
    • 水準: 要因が取りうる様々な状態のこと
    • 水準数: 水準のパターン数。以下要因Aの水準数をaで表す。
  • 独立した1要因モデル

    • モデル
      • y_{ij} = \mu_j + e
      •  e \sim Normal(0, \sigma_e)
      • 要因のj番目の水準におけるi番目の測定値が、j番目の水準がもつ平均的な効果と正規分布から生じる乱数による現れる
    • 以上を纏めて、事後確率はf(y_{ij} | \mu_j, \sigma_e)と表される
    • 水準内の測定値は独立に観測されてると考える 
      • 同時確率は f(y_{ij} | \mu_j, \sigma_e) = \Pi_{i = 1}^n f(y_{ij}|\mu_j, \sigma_e)
      • 即ち、尤度は f(y|\theta) =\Pi_{i = 1}^a f(y_j|\mu_j, \sigma_e)
  • 生成量

    • 全平均 :  \mu = \Sigma_{i=1}^a \frac{\mu}{a}
      • 各水準の平均の単純な平均
    • 水準の効果(eol): a_j = \mu_j - \mu
      • 各平均から全平均を引いたもの
      • 水準の効果の和が0になることに注意
    • 説明率:  \eta^{2} = \frac{\sigma_{a}^{2}}{\sigma_{y}^{2}}
      • 測定値の分散に占める要因の分散の割合で、要因が観測値をどれほど説明するかを示す
      •  \sigma_{y}^{2}は効果と誤差が互いに独立ならば、 \sigma_{y}^{2} = \sigma_{a}^{2} + \sigma_{e}^{2}である
        • 測定値の分散は要因の分散 + 誤差の分散になるということ
      •  0 \leq \eta \leq 1で定義される
    • 効果量 :  \delta = \frac{\sigma_a}{\sigma_e}
      • 水準間の標準偏差が、水準内の標準偏差の何倍に相当するか示すもの
      •  \sigma_eは要因の影響を受けない、水準内の平均的バラつきと等価であると考える事ができる
  • モデルを使った要因の効果の方法例

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

    • 要因を複数導入
      • 水準の組み合わせ1単位のことをセルと以下では呼称する
    • 交互作用: 一方の要因の水準毎に、他方の要因の平均の高低パターンが異なる現象
    • モデル
      •  y_{ijk} = \mu  + a_{j} + b_{k} + ab_{jk} + e
      •  e \sim Normal(0, \sigma_e)
      • 水準 j, k i番目の測定値が、全体の平均に要因 a j番目の水準の効果 + 要因 b k番目の水準の効果 +  a bの交互作用の効果 + 乱数から構成されることを示す
      • セルの平均 + 乱数
    • それぞれ全平均(水準の効果)の性質から次の制約を持つ(水準の効果を全水準で足すと0になる性質を定式化)
      •  \Sigma_{i=1}^{a} a_{i} = 0
      •  \Sigma_{i=1}^{b} b_{i} = 0
      •  \Sigma_{j=1}^{a} ab_{jk} = 0
      •  \Sigma_{k=1}^{a} ab_{jk} = 0
    • 尤度 f(y|\theta) = f(y|\mu, \sigma_e = \Pi_{j=1}^{a}\Pi_{k=1}^{b} f(y_{jk} | \mu_{jk}, \sigma_e)

章末問題

1要因の推測

データの取得
import numpy as np
import os
from logging import getLogger, Formatter, StreamHandler, DEBUG

# printではなくloggerを使う
def get_logger():
    logger = getLogger(__name__)
    logger.setLevel(DEBUG)
    log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s'
    formatter = Formatter(log_fmt)
    stream_handler = StreamHandler()
    stream_handler.setLevel(DEBUG)
    stream_handler.setFormatter(formatter)
    logger.addHandler(stream_handler)
    return logger
logger = get_logger()

# データの準備
ld = np.array([5.02, 6.67, 8.17, 2.79, 8.13, 6.34, 6.32, 3.97])
ll = np.array([9.89, 9.58, 11.20, 9.05, 12.33, 9.39, 10.88, 9.37, 17.40])
dm = np.array([10.20, 7.29, 7.57, 3.42, 5.82, 10.92, 5.21, 13.47, 8.64, 6.05])
y = np.concatenate([a, b, c])
A_ld = np.full(ld.size, 1, np.int)
A_ll = np.full(ll.size, 2, np.int)
A_dm = np.full(dm.size, 3, np.int)
A = np.concatenate([A_ld, A_ll, A_dm])

# 出力用ディレクトリの用意
result_dir = os.path.join("result", "chapter5")
if os.path.exists(result_dir) is False:
    os.makedirs(result_dir)
    logger.info("{0} is made".format(result_dir))
Stanによる統計モデルの構築

特に難しいところは無いような気がします。説明量を計算するために、分散値を計算する手間があるぐらいでしょうか?

import os
import pystan
import pickle

# Stanのモデルを読み込んでコンパイルする
stan_file = os.path.join("stan", "e1_ind.stan")
stan_file_c = os.path.join("stan", "e1_ind.pkl")
model = pystan.StanModel(file=stan_file)
with open(stan_file_c, "wb") as f:
    pickle.dump(model, f)
logger.info("Stan model is compiled to {0}".format(stan_file_c))

e1_ind.stan

data {
    int<lower=0> n ;
    int<lower=0> a ;
    real<lower=0> y[n] ;
    int<lower=1, upper=3> A[n] ;
}

parameters {
    vector[a] mu_a ;
    real<lower=0> sigma_e ;
}

model {
    for(i in 1:n){
        y[i] ~ normal(mu_a[A[i]], sigma_e) ;
    }
}

generated quantities {
    real mu ;
    vector[a] eol ;
    vector[a] eol2 ;
    vector[a] prob_eol_upper_0 ;
    real var_e ;
    real var_a ;
    real eta2 ;
    real delta ;
    int<lower=0, upper=1> prob_ll_upper_dm ;
    int<lower=0, upper=1> prob_ll_upper_ld ;
    int<lower=0, upper=1> prob_dm_upper_ld ;
    vector[n] log_lik ;

    mu = mean(mu_a) ;
    for(i in 1:a){
        eol[i] = mu_a[i] - mu ;
        eol2[i] = pow(eol[i], 2) ;
        prob_eol_upper_0[i] = eol[i] > 0 ? 1 : 0 ;
    }
    var_a = mean(eol2) ;
    var_e = pow(sigma_e, 2) ;
    eta2 = var_a / (var_a + var_e) ;
    delta = sqrt(var_a) / sigma_e ;
    prob_ll_upper_dm = mu_a[2] - mu_a[3] > 0 ? 1 : 0 ;
    prob_ll_upper_ld = mu_a[2] - mu_a[1] > 0 ? 1 : 0 ;
    prob_dm_upper_ld = mu_a[3] - mu_a[1] > 0 ? 1 : 0 ;
    for(i in 1:n){
        log_lik[i] = normal_lpdf(y[i] | mu_a[A[i]], sigma_e) ;
    }
}
統計モデルによる事後分布のサンプリング
import pandas as pd
import pickle
import pystan
import matplotlib
import matplotlib.pyplot as plt
import os
from IPython.core.display import display
from tabulate import tabulate
matplotlib.rcParams['figure.figsize'] = (10, 50)

# Stanで使うデータの用意
stan_data = {"n": y.size,
             "a": 3,
             "y": y,
             "A": A}
# 興味のあるパラメータの設定
par = ["mu_a",
       "sigma_e",
       "eol",
       "eta2",
       "delta",
       "prob_eol_upper_0",
       "prob_ll_upper_dm",
       "prob_ll_upper_ld",
       "prob_dm_upper_ld",
       "log_lik"]
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

# モデルの読み込み
stan_file_c = os.path.join("stan", "e1_ind.pkl")
with open(stan_file_c, "rb") as f:
    model = pickle.load(f)

# MCMCでサンプリング
logger.info("Start MCMC sampling")
fit = model.sampling(data=stan_data,
                     pars=par,
                     iter=21000,
                     chains=5,
                     warmup=1000,
                     seed=1234,
                     algorithm="NUTS")


# 事後分布の表を取得
summary = fit.summary(pars=par, probs=prob)
summary_df = pd.DataFrame(summary["summary"],
                          index=summary["summary_rownames"],
                          columns=summary["summary_colnames"])
display(summary_df)
summary_df_path = os.path.join(result_dir, "df_summary_1.md")
with open(summary_df_path, "w") as f:
    f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe"))
logger.info("MCMC result summary is saved at {0}".format(summary_df_path))

# 事後分布の可視化
fit.traceplot(par)
traceplot_path = os.path.join(result_dir, "traceplot_1.png")
plt.savefig(traceplot_path)
plt.show()
logger.info("traceplot result is saved at {0}".format(traceplot_path))

# WAICの計算
log_lik = fit.extract("log_lik")["log_lik"]
waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))
logger.info("WAICの値は{0}です".format(waic))
  • 結果(表と分布)
mean se_mean sd 2.5% 5% 25% 50% 75% 95% 97.5% n_eff Rhat
mu_a[0] 5.92465 0.0030795 0.973823 3.99253 4.32772 5.28835 5.92998 6.56208 7.51919 7.84787 100000 1.00007
mu_a[1] 11.0072 0.00292346 0.924478 9.17396 9.49254 10.4026 11.0076 11.6092 12.5191 12.8357 100000 0.999985
mu_a[2] 7.85955 0.00276635 0.874797 6.12594 6.42644 7.28474 7.86113 8.43353 9.29189 9.58305 100000 1
sigma_e 2.74451 0.00151073 0.42629 2.06146 2.14526 2.44184 2.69385 2.9907 3.51953 3.72462 79623 0.999968
eol[0] -2.33916 0.00245045 0.7749 -3.88038 -3.61524 -2.8447 -2.33923 -1.83011 -1.07629 -0.817259 100000 1.00008
eol[1] 2.74342 0.00239336 0.756846 1.24128 1.50518 2.24837 2.74077 3.23758 3.98729 4.24112 100000 0.999999
eol[2] -0.404261 0.00233393 0.738053 -1.86101 -1.60891 -0.887243 -0.405952 0.0805197 0.804455 1.05156 100000 1.00003
eta2 0.384362 0.000399327 0.126278 0.122444 0.162176 0.299146 0.391642 0.476528 0.579842 0.608479 100000 0.999981
delta 0.804616 0.000704594 0.222812 0.373535 0.439964 0.653323 0.802353 0.954108 1.17476 1.24665 100000 0.999984
prob_eol_upper_0[0] 0.00215 0.000165199 0.0463185 0 0 0 0 0 0 0 78613 1.00006
prob_eol_upper_0[1] 0.99963 6.2415e-05 0.0192319 1 1 1 1 1 1 1 94944 1.00003
prob_eol_upper_0[2] 0.28713 0.00152318 0.452425 0 0 0 0 1 1 1 88225 1.00002
prob_ll_upper_dm 0.99168 0.000333967 0.0908342 1 1 1 1 1 1 1 73976 0.999978
prob_ll_upper_ld 0.99975 5.19414e-05 0.0158095 1 1 1 1 1 1 1 92642 0.999991
prob_dm_upper_ld 0.93277 0.000930336 0.250421 0 0 1 1 1 1 1 72454 1.00007
log_lik[0] -2.03678 0.000758605 0.199372 -2.49373 -2.39575 -2.14793 -2.01548 -1.89931 -1.75258 -1.7086 69071 0.999994
log_lik[1] -2.01783 0.000754794 0.191828 -2.4532 -2.35815 -2.12696 -1.99886 -1.88578 -1.74113 -1.69764 64590 1.00002
log_lik[2] -2.33612 0.000992587 0.313883 -3.08375 -2.92534 -2.50565 -2.2841 -2.11391 -1.92039 -1.86752 100000 1.00004
log_lik[3] -2.67679 0.00135735 0.429231 -3.69657 -3.47937 -2.91461 -2.60959 -2.36613 -2.10742 -2.04192 100000 1.00004
log_lik[4] -2.32348 0.000978475 0.309421 -3.06092 -2.90399 -2.49026 -2.27211 -2.10467 -1.91389 -1.86127 100000 1.00004
log_lik[5] -1.99069 0.000735103 0.179753 -2.39033 -2.3073 -2.0972 -1.97534 -1.86648 -1.72545 -1.6824 59794 1.00001
log_lik[6] -1.98954 0.00073424 0.17923 -2.38764 -2.30472 -2.09579 -1.9742 -1.86567 -1.72457 -1.6816 59586 1.00001
log_lik[7] -2.25015 0.000893585 0.282577 -2.92974 -2.78158 -2.40036 -2.20673 -2.05154 -1.87364 -1.82393 100000 1.00003
log_lik[8] -2.06094 0.000756641 0.203617 -2.52871 -2.42581 -2.17552 -2.0385 -1.91954 -1.77099 -1.72765 72418 1
log_lik[9] -2.11696 0.000790546 0.223723 -2.63793 -2.52526 -2.2399 -2.08855 -1.96137 -1.80487 -1.76051 80088 0.999999
log_lik[10] -1.9749 0.000705525 0.171411 -2.3497 -2.27598 -2.07886 -1.96221 -1.85611 -1.71601 -1.67376 59027 1.00001
log_lik[11] -2.24435 0.000846443 0.267669 -2.87634 -2.74126 -2.38938 -2.20473 -2.05568 -1.88209 -1.83201 100000 0.999994
log_lik[12] -2.09647 0.000768147 0.216744 -2.60167 -2.49056 -2.21654 -2.06987 -1.94642 -1.79251 -1.74675 79617 1
log_lik[13] -2.15804 0.000816763 0.238127 -2.71684 -2.59706 -2.2877 -2.12586 -1.9917 -1.83017 -1.78421 85001 0.999997
log_lik[14] -1.97342 0.000706486 0.170781 -2.34448 -2.27222 -2.07722 -1.96171 -1.85505 -1.71519 -1.67338 58435 1.00001
log_lik[15] -2.16266 0.000819785 0.239732 -2.72544 -2.60451 -2.29303 -2.13016 -1.99503 -1.83272 -1.78664 85517 0.999997
log_lik[16] -4.87403 0.00339268 1.07286 -7.32982 -6.8243 -5.5138 -4.74535 -4.09518 -3.35255 -3.15931 100000 0.999977
log_lik[17] -2.35558 0.000919091 0.290642 -3.03811 -2.89536 -2.51793 -2.31435 -2.14877 -1.95975 -1.90793 100000 1
log_lik[18] -1.98969 0.000680232 0.173859 -2.37117 -2.2949 -2.09513 -1.97689 -1.86882 -1.72784 -1.68611 65325 0.999965
log_lik[19] -1.97261 0.000674084 0.168196 -2.33564 -2.26748 -2.07642 -1.96134 -1.85564 -1.7164 -1.67546 62259 0.999965
log_lik[20] -3.36627 0.00188208 0.595165 -4.73815 -4.46265 -3.7144 -3.2838 -2.93321 -2.54231 -2.44062 100000 0.999985
log_lik[21] -2.26206 0.000825123 0.260927 -2.87464 -2.74387 -2.40867 -2.22729 -2.07788 -1.90249 -1.85212 100000 0.999979
log_lik[22] -2.6317 0.00118696 0.37535 -3.50938 -3.32991 -2.84462 -2.57612 -2.36118 -2.12327 -2.05933 100000 1
log_lik[23] -2.46518 0.00102346 0.323646 -3.2244 -3.0644 -2.64582 -2.41881 -2.23303 -2.02502 -1.96827 100000 0.999983
log_lik[24] -4.20173 0.00267471 0.845818 -6.14049 -5.74926 -4.70266 -4.09605 -3.58666 -3.0153 -2.86741 100000 0.99999
log_lik[25] -2.00989 0.000689887 0.180727 -2.40828 -2.32687 -2.11763 -1.99477 -1.88481 -1.74166 -1.69933 68626 0.999981
log_lik[26] -2.19919 0.000762851 0.241235 -2.76133 -2.64207 -2.33492 -2.16884 -2.03002 -1.86336 -1.81479 100000 0.999977

f:id:ajhjhaf:20170531233354p:plain

2要因の推測

データの取得

データの取得がめちゃめちゃ汚い気がしますが、許してください…

import numpy as np
import os
from logging import getLogger, Formatter, StreamHandler, DEBUG

# printではなくloggerを使う
def get_logger():
    logger = getLogger(__name__)
    logger.setLevel(DEBUG)
    log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s'
    formatter = Formatter(log_fmt)
    stream_handler = StreamHandler()
    stream_handler.setLevel(DEBUG)
    stream_handler.setFormatter(formatter)
    logger.addHandler(stream_handler)
    return logger
logger = get_logger()

# データの準備
y = np.array([
140,146,149,136,147,147,143,143,143,141,
139,136,136,140,135,132,140,134,
123,127,131,130,138,128,129,
115,120,118,118,121,124,129,119,128,
128,124,123,121,122,126,131,122,
121,121,120,116,117,113,118,
143,141,142,145,149,145,143,141,142,155,
138,134,142,136,135,136,131,133,
131,128,128,128,127,130,130,
117,125,132,122,119,122,129,117,127,
117,120,124,122,122,122,118,122,
119,125,122,116,119,113,122])
A_loaded = np.full(49, 1, np.int)
A_empty = np.full(y.size - A_loaded.size, 2, np.int)
A = np.concatenate([A_loaded, A_empty])
B_loaded_st = np.full(10, 1, np.int)
B_loaded_cut = np.full(8, 2, np.int)
B_loaded_f = np.full(7, 3, np.int)
B_loaded_ch = np.full(9, 4, np.int)
B_loaded_sl = np.full(8, 5, np.int)
B_loaded_cur = np.full(7, 6, np.int)
B_empty_st = np.full(10, 1, np.int)
B_empty_cut = np.full(8, 2, np.int)
B_empty_f = np.full(7, 3, np.int)
B_empty_ch = np.full(9, 4, np.int)
B_empty_sl = np.full(8, 5, np.int)
B_empty_cur = np.full(7, 6, np.int)
B = np.concatenate([B_loaded_st,
                    B_loaded_cut,
                    B_loaded_f,
                    B_loaded_ch,
                    B_loaded_sl,
                    B_loaded_cur,
                    B_empty_st,
                    B_empty_cut,
                    B_empty_f,
                    B_empty_ch,
                    B_empty_sl,
                    B_empty_cur])


# 出力用ディレクトリの用意
result_dir = os.path.join("result", "chapter5")
if os.path.exists(result_dir) is False:
    os.makedirs(result_dir)
    logger.info("{0} is made".format(result_dir))
Stanによる統計モデルの構築
import os
import pystan
import pickle

# Stanのモデルを読み込んでコンパイルする
stan_file = os.path.join("stan", "e2_ind.stan")
stan_file_c = os.path.join("stan", "e2_ind.pkl")
model = pystan.StanModel(file=stan_file)
with open(stan_file_c, "wb") as f:
    pickle.dump(model, f)
logger.info("Stan model is compiled to {0}".format(stan_file_c))

e2_ind.stan

  • PyStanがmatrixを可視化する際にうまくうけとってくれないので、ベクトルに変換してからPythonへ渡しています
  • 本の付録のStanコードではtransformed parametersの部分が冗長だったので、forとif文で対応しています。
data {
    int<lower=1> n ;
    int<lower=1> n_a ;
    int<lower=1> n_b ;
    int<lower=0> y[n] ;
    int<lower=1, upper=n_a> A[n] ;
    int<lower=1, upper=n_b> B[n] ;
}

parameters {
    real<lower=0> mu ;
    vector[n_a-1] a ;
    vector[n_b-1] b ;
    matrix[n_a-1, n_b-1] ab ;
    real<lower=0> sigma_e ;
}

transformed parameters {
    # 水準の効果の和が0になる制約条件を使う
    vector[n_a] a_r ;
    vector[n_b] b_r ;
    matrix[n_a, n_b] ab_r ;
    vector[n_a*n_b] ab_v ;

    for(i in 1:n_a-1){
        a_r[i] = a[i] ;
    }
    a_r[n_a] = - sum(a_r[1:n_a - 1]) ;

    for(i in 1:(n_b-1)){
        b_r[i] = b[i] ;
    }
    b_r[n_b] = - sum(b_r[1:n_b - 1]) ;

    for (j in 1:n_a) {
        for (k in 1:n_b) {
            if (j == n_a && k == n_b) {
                ab_r[j, k] = sum(ab[1:j-1, 1:k-1]) ;
            } else if (j == n_a){
                ab_r[j, k] = -sum(ab[1:j-1, k]) ;
            } else if (k == n_b){
                ab_r[j, k] = -sum(ab[j, 1:k-1]) ;
            } else {
                ab_r[j, k] = ab[j, k] ;
            }
            ab_v[(j-1) * n_b + k] = ab_r[j, k] ;
        }
    }
}

model {
    for (i in 1:n) {
        y[i] ~ normal(mu + a_r[A[i]] + b_r[B[i]] + ab_r[A[i], B[i]], sigma_e) ;
    }
}

generated quantities {
    vector[n] log_lik ;
    vector[n_b] prob_b_upper_0 ;
    vector[n_a] a_r2 ;
    vector[n_b] b_r2 ;
    vector[n_a*n_b] ab_r2 ;
    real var_a ;
    real var_b ;
    real var_ab ;
    real var_e ;
    real var_y ;
    real eta_a2 ;
    real eta_b2 ;
    real eta_ab2 ;
    real eta_t2 ;
    real delta_a ;
    real delta_b ;
    real delta_ab ;
    matrix[n_b, n_b] prob_b_comparison ;

    for (i in 1:n) {
        log_lik[i] = normal_lpdf(y[i] | mu + a_r[A[i]] + b_r[B[i]] + ab_r[A[i], B[i]], sigma_e) ;
    }
    for (i in 1:n_b) {
        prob_b_upper_0[i] = b_r[i] > 0 ? 1 : 0 ;
    }
    for (i in 1:n_a) {
        a_r2[i] = pow(a_r[i], 2) ;
    }
    for (i in 1:n_b) {
        b_r2[i] = pow(b_r[i], 2) ;
    }
    for (i in 1:n_a*n_b) {
        ab_r2[i] = pow(ab_v[i], 2) ;
    }
    var_a = mean(a_r2) ;
    var_b = mean(b_r2) ;
    var_ab = mean(ab_r2) ;
    var_e = pow(sigma_e, 2) ;
    var_y = var_a + var_b + var_ab + var_e ;
    eta_a2 = var_a / var_y ;
    eta_b2 = var_b / var_y ;
    eta_ab2 = var_ab / var_y ;
    eta_t2 = (var_a + var_b + var_ab) / var_y ;
    delta_a = sqrt(var_a) / sigma_e ;
    delta_b = sqrt(var_b) / sigma_e ;
    delta_ab = sqrt(var_ab) / sigma_e ;

    for (j in 1:n_b){
        for (k in 1:n_b){
            prob_b_comparison[j, k] = b_r[j] > b_r[k] ? 1 : 0 ;
        }
    }
}
統計モデルによる事後分布のサンプリング
  • データ型が行列のパラメータは、可視化する前に除いています
  • 水準の効果の差については、別途ピボットテーブルを作成して可視化しました。
import pandas as pd
import pickle
import pystan
import matplotlib
import matplotlib.pyplot as plt
import os
from IPython.core.display import display
from tabulate import tabulate
matplotlib.rcParams['figure.figsize'] = (10, 50)

# Stanで使うデータの用意
stan_data = {"n": y.size,
             "n_a": np.unique(A).size,
             "n_b": np.unique(B).size,
             "y": y,
             "A": A,
             "B": B}
# 興味のあるパラメータの設定
pars = ["mu",
       "sigma_e",
       "a_r",
       "b_r",
       "ab_v",
       "prob_b_upper_0",
       "eta_a2",
       "eta_b2",
       "eta_ab2",
       "eta_t2",
       "delta_a",
       "delta_b",
       "delta_ab",
       "prob_b_comparison",
       "log_lik"]
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

# モデルの読み込み
stan_file_c = os.path.join("stan", "e2_ind.pkl")
with open(stan_file_c, "rb") as f:
    model = pickle.load(f)

# MCMCでサンプリング
logger.info("Start MCMC sampling")
fit = model.sampling(data=stan_data,
                     pars=pars,
                     iter=21000,
                     chains=5,
                     warmup=1000,
                     seed=1234,
                     algorithm="NUTS")


# 事後分布の表を取得
summary = fit.summary(pars=pars, probs=prob)
summary_df = pd.DataFrame(summary["summary"],
                          index=summary["summary_rownames"],
                          columns=summary["summary_colnames"])
display(summary_df)
summary_df_path = os.path.join(result_dir, "df_summary_2.md")
with open(summary_df_path, "w") as f:
    f.write(tabulate(summary_df, summary_df.columns, tablefmt="pipe"))
logger.info("MCMC result summary is saved at {0}".format(summary_df_path))

# 水準毎に上回る確率を調べる
comparison = fit.summary(pars=["prob_b_comparison"], probs=prob)
comparison_df = pd.DataFrame(comparison["summary"],
                             index=comparison["summary_rownames"],
                             columns=comparison["summary_colnames"])
comparison_index = pd.Series(comparison_df.index.str.strip("prob_b_comparison[").str.strip("]").str.split(","))
b_dict = {"0":"straight", "1":"cutball", "2":"fork", "3":"changeup", "4":"slider", "5":"curve"}
upper_index = comparison_index.apply(lambda x:b_dict[x[0]])
lower_index = comparison_index.apply(lambda x:b_dict[x[1]])
comparison_triple = pd.DataFrame([upper_index, lower_index, comparison_df["mean"].values]).T
comparison_triple.columns = ["upper", "lower", "possibility"]
comparison_pivot = comparison_triple.pivot(index="upper", columns="lower", values="possibility")
display(comparison_pivot)
b_df_path = os.path.join(result_dir, "b_comparison.md")
with open(b_df_path, "w") as f:
    f.write(tabulate(comparison_pivot, comparison_pivot.columns, tablefmt="pipe"))
logger.info("comparison table is saved at {0}".format(b_df_path))

# 事後分布の可視化
pars.remove("prob_b_upper_0")
pars.remove("prob_b_comparison")
fit.traceplot(pars)
traceplot_path = os.path.join(result_dir, "traceplot_2.png")
plt.savefig(traceplot_path)
plt.show()
logger.info("traceplot result is saved at {0}".format(traceplot_path))

# WAICの計算
log_lik = fit.extract("log_lik")["log_lik"]
waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))
logger.info("WAICの値は{0}です".format(waic))
  • 結果(表と分布)
mean se_mean sd 2.5% 5% 25% 50% 75% 95% 97.5% n_eff Rhat
mu 128.842 0.00125739 0.397622 128.059 128.188 128.577 128.842 129.107 129.495 129.622 100000 0.999982
sigma_e 3.89006 0.000955725 0.302227 3.35328 3.42584 3.67831 3.87124 4.08093 4.41686 4.53412 100000 1.00008
a_r[0] 0.0566725 0.00125659 0.39737 -0.722588 -0.595841 -0.208101 0.0562403 0.321384 0.710725 0.840036 100000 0.999989
a_r[1] -0.0566725 0.00125659 0.39737 -0.840036 -0.710725 -0.321384 -0.0562403 0.208101 0.595841 0.722588 100000 0.999989
b_r[0] 15.2089 0.00256537 0.811242 13.6165 13.8763 14.6635 15.2074 15.7524 16.5397 16.8083 100000 0.999994
b_r[1] 7.22147 0.00280979 0.888533 5.47369 5.76311 6.62748 7.22389 7.81345 8.68103 8.97568 100000 1
b_r[2] 0.303162 0.00295644 0.93491 -1.53609 -1.2345 -0.321006 0.298546 0.928776 1.83405 2.15155 100000 0.999977
b_r[3] -6.50771 0.00270552 0.85556 -8.17839 -7.90278 -7.08083 -6.51271 -5.93491 -5.10512 -4.82134 100000 0.999979
b_r[4] -6.09233 0.00281005 0.888616 -7.83205 -7.54945 -6.68611 -6.09106 -5.49716 -4.62802 -4.3491 100000 0.999976
b_r[5] -10.1335 0.00296533 0.93772 -11.9753 -11.6719 -10.7628 -10.1342 -9.49813 -8.59766 -8.30235 100000 0.999964
ab_v[0] -0.607631 0.00258641 0.817893 -2.21712 -1.9559 -1.15031 -0.607124 -0.0654354 0.744995 1.00949 100000 0.999985
ab_v[1] 0.381475 0.00281664 0.8907 -1.37495 -1.08112 -0.212969 0.382057 0.976252 1.84594 2.13356 100000 0.999998
ab_v[2] 0.23221 0.00297549 0.940934 -1.60767 -1.31047 -0.399151 0.229995 0.866498 1.77931 2.08061 100000 0.999969
ab_v[3] -1.05319 0.00268952 0.8505 -2.72909 -2.45038 -1.62572 -1.05216 -0.482988 0.345127 0.618895 100000 0.99998
ab_v[4] 1.81474 0.00281049 0.888756 0.0653521 0.350221 1.22523 1.81277 2.40855 3.27038 3.56064 100000 0.999986
ab_v[5] -0.767603 0.00297959 0.94223 -2.60335 -2.31877 -1.40147 -0.769324 -0.130914 0.778092 1.07769 100000 0.999986
ab_v[6] 0.607631 0.00258641 0.817893 -1.00949 -0.744995 0.0654354 0.607124 1.15031 1.9559 2.21712 100000 0.999985
ab_v[7] -0.381475 0.00281664 0.8907 -2.13356 -1.84594 -0.976252 -0.382057 0.212969 1.08112 1.37495 100000 0.999998
ab_v[8] -0.23221 0.00297549 0.940934 -2.08061 -1.77931 -0.866498 -0.229995 0.399151 1.31047 1.60767 100000 0.999969
ab_v[9] 1.05319 0.00268952 0.8505 -0.618895 -0.345127 0.482988 1.05216 1.62572 2.45038 2.72909 100000 0.99998
ab_v[10] -1.81474 0.00281049 0.888756 -3.56064 -3.27038 -2.40855 -1.81277 -1.22523 -0.350221 -0.0653521 100000 0.999986
ab_v[11] 0.767603 0.00297959 0.94223 -1.07769 -0.778092 0.130914 0.769324 1.40147 2.31877 2.60335 100000 0.999986
prob_b_upper_0[0] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_upper_0[1] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_upper_0[2] 0.62685 0.00152942 0.483644 0 0 0 1 1 1 1 100000 1
prob_b_upper_0[3] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_upper_0[4] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_upper_0[5] 0 0 0 0 0 0 0 0 0 0 100000 nan
eta_a2 0.00168409 1.10406e-05 0.00242859 1.69709e-06 6.70442e-06 0.000168269 0.000753203 0.00220135 0.00650296 0.00855088 48386 1.00007
eta_b2 0.820698 8.31697e-05 0.0263006 0.762582 0.77391 0.804681 0.822937 0.83927 0.859708 0.865594 100000 1.00006
eta_ab2 0.0179692 3.28066e-05 0.00948945 0.00401563 0.00533333 0.0109958 0.0164943 0.0233222 0.0357468 0.0404685 83668 1.00012
eta_t2 0.840352 7.44251e-05 0.0235353 0.788162 0.798497 0.826116 0.842443 0.857014 0.875117 0.880387 100000 1.00003
delta_a 0.0820056 0.000277709 0.0621218 0.0032616 0.00655807 0.032731 0.0691728 0.117879 0.201804 0.231203 50039 1.00011
delta_b 2.28756 0.000631375 0.199658 1.90128 1.96222 2.15243 2.28528 2.42122 2.62019 2.68599 100000 1.00004
delta_ab 0.325882 0.000285382 0.0902456 0.157327 0.181951 0.262588 0.3238 0.386283 0.478362 0.50898 100000 1.00008
prob_b_comparison[0,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[1,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[2,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[3,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[4,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[5,0] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[0,1] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[1,1] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[2,1] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[3,1] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[4,1] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[5,1] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[0,2] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[1,2] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[2,2] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[3,2] 1e-05 1e-05 0.00316228 0 0 0 0 0 0 0 100000 1
prob_b_comparison[4,2] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[5,2] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[0,3] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[1,3] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[2,3] 0.99999 1e-05 0.00316228 1 1 1 1 1 1 1 100000 1
prob_b_comparison[3,3] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[4,3] 0.62121 0.00153398 0.485088 0 0 0 1 1 1 1 100000 0.999983
prob_b_comparison[5,3] 0.00526 0.000240106 0.0723352 0 0 0 0 0 0 0 90760 0.999992
prob_b_comparison[0,4] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[1,4] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[2,4] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[3,4] 0.37879 0.00153398 0.485088 0 0 0 0 1 1 1 100000 0.999983
prob_b_comparison[4,4] 0 0 0 0 0 0 0 0 0 0 100000 nan
prob_b_comparison[5,4] 0.00256 0.000167342 0.0505319 0 0 0 0 0 0 0 91184 1
prob_b_comparison[0,5] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[1,5] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[2,5] 1 0 0 1 1 1 1 1 1 1 100000 nan
prob_b_comparison[3,5] 0.99474 0.000240106 0.0723352 1 1 1 1 1 1 1 90760 0.999992
prob_b_comparison[4,5] 0.99744 0.000167342 0.0505319 1 1 1 1 1 1 1 91184 1
prob_b_comparison[5,5] 0 0 0 0 0 0 0 0 0 0 100000 nan
log_lik[0] -2.73633 0.000935194 0.295734 -3.43849 -3.28893 -2.90328 -2.68784 -2.5182 -2.34766 -2.30554 100000 1.00002
log_lik[1] -2.53453 0.00070078 0.221606 -3.08005 -2.95548 -2.65052 -2.49078 -2.37503 -2.25745 -2.22527 100000 0.999975
log_lik[2] -3.34177 0.00146612 0.463629 -4.39499 -4.1848 -3.6203 -3.28603 -3.00338 -2.6859 -2.6025 100000 0.999994
log_lik[3] -4.21624 0.00207097 0.654899 -5.66336 -5.38642 -4.62001 -4.15503 -3.74499 -3.25356 -3.11148 100000 0.999974
log_lik[4] -2.73634 0.000936502 0.296148 -3.44194 -3.28963 -2.90404 -2.68719 -2.51799 -2.34749 -2.30663 100000 0.999975
log_lik[5] -2.73634 0.000936502 0.296148 -3.44194 -3.28963 -2.90404 -2.68719 -2.51799 -2.34749 -2.30663 100000 0.999975
log_lik[6] -2.33272 0.000426818 0.111143 -2.59392 -2.52945 -2.39027 -2.31979 -2.25852 -2.17872 -2.15482 67807 1.0002
log_lik[7] -2.33272 0.000426818 0.111143 -2.59392 -2.52945 -2.39027 -2.31979 -2.25852 -2.17872 -2.15482 67807 1.0002
log_lik[8] -2.33272 0.000426818 0.111143 -2.59392 -2.52945 -2.39027 -2.31979 -2.25852 -2.17872 -2.15482 67807 1.0002
log_lik[9] -2.53452 0.00069946 0.221189 -3.07629 -2.95565 -2.65031 -2.49165 -2.37471 -2.25725 -2.22585 100000 1.00006
log_lik[10] -2.54659 0.000785481 0.248391 -3.16542 -3.02514 -2.67084 -2.4917 -2.36823 -2.24892 -2.2172 100000 1.00005
log_lik[11] -2.34498 0.000496465 0.125425 -2.64993 -2.57005 -2.40342 -2.3269 -2.26265 -2.18144 -2.15687 63825 1.00009
log_lik[12] -2.34498 0.000496465 0.125425 -2.64993 -2.57005 -2.40342 -2.3269 -2.26265 -2.18144 -2.15687 63825 1.00009
log_lik[13] -2.74833 0.00104876 0.331647 -3.546 -3.37344 -2.93064 -2.68727 -2.50315 -2.32712 -2.28515 100000 1.00003
log_lik[14] -2.41232 0.000613607 0.176642 -2.86322 -2.75229 -2.4868 -2.3749 -2.29388 -2.20328 -2.1765 82872 1.00003
log_lik[15] -3.01794 0.00133705 0.422813 -4.00776 -3.80561 -3.25894 -2.95459 -2.70599 -2.4528 -2.39164 100000 1.00002
log_lik[16] -2.74833 0.00104876 0.331647 -3.546 -3.37344 -2.93064 -2.68727 -2.50315 -2.32712 -2.28515 100000 1.00003
log_lik[17] -2.54692 0.000791085 0.250163 -3.16951 -3.02844 -2.67035 -2.49157 -2.36771 -2.24974 -2.21829 100000 1.00002
log_lik[18] -3.73765 0.00205657 0.650345 -5.21599 -4.92322 -4.1275 -3.66109 -3.2632 -2.81517 -2.69676 100000 0.999991
log_lik[19] -2.5447 0.000828546 0.262009 -3.21331 -3.05791 -2.67094 -2.48239 -2.35768 -2.24115 -2.21039 100000 1.00001
log_lik[20] -2.42805 0.000693641 0.195391 -2.93941 -2.809 -2.50779 -2.38289 -2.29826 -2.2048 -2.17837 79349 1.00005
log_lik[21] -2.3563 0.000560908 0.138402 -2.70641 -2.61165 -2.4159 -2.33316 -2.26661 -2.1839 -2.15899 60884 1.00014
log_lik[22] -4.81377 0.00283885 0.897723 -6.79217 -6.42019 -5.36376 -4.72578 -4.16988 -3.49737 -3.30988 100000 0.999977
log_lik[23] -2.41463 0.000656152 0.184738 -2.90055 -2.776 -2.48901 -2.37315 -2.29235 -2.20083 -2.17451 79269 1.00006
log_lik[24] -2.35183 0.000548798 0.133508 -2.68773 -2.59816 -2.41028 -2.33023 -2.26506 -2.18237 -2.15774 59182 1.00015
log_lik[25] -3.68089 0.00179969 0.569113 -4.95951 -4.71565 -4.02176 -3.61904 -3.26868 -2.86285 -2.7564 100000 0.99997
log_lik[26] -2.39015 0.000531975 0.156079 -2.78711 -2.69127 -2.45852 -2.36012 -2.28619 -2.19753 -2.17168 86081 1.00001
log_lik[27] -2.70464 0.000949753 0.300338 -3.42914 -3.27577 -2.86761 -2.65075 -2.48359 -2.32096 -2.28113 100000 0.999965
log_lik[28] -2.70464 0.000949753 0.300338 -3.42914 -3.27577 -2.86761 -2.65075 -2.48359 -2.32096 -2.28113 100000 0.999965
log_lik[29] -2.3338 0.000448643 0.114007 -2.60405 -2.53915 -2.39156 -2.31964 -2.25791 -2.17838 -2.15378 64574 1.00008
log_lik[30] -2.56839 0.000779171 0.246395 -3.17826 -3.04225 -2.69622 -2.51769 -2.39041 -2.2649 -2.23265 100000 0.999982
log_lik[31] -4.30475 0.00222789 0.70452 -5.86668 -5.56526 -4.73925 -4.23866 -3.80034 -3.26551 -3.12075 100000 0.999971
log_lik[32] -2.51376 0.000705027 0.222949 -3.07066 -2.94619 -2.62376 -2.46574 -2.35573 -2.24338 -2.21269 100000 0.999976
log_lik[33] -3.82294 0.00190039 0.600955 -5.17464 -4.90752 -4.18988 -3.75953 -3.39001 -2.94919 -2.83484 100000 0.999968
log_lik[34] -2.72106 0.00102137 0.322986 -3.50706 -3.33946 -2.89458 -2.66028 -2.48187 -2.31426 -2.27516 100000 0.99998
log_lik[35] -2.34996 0.000499418 0.129235 -2.67219 -2.58849 -2.40913 -2.32994 -2.2649 -2.18353 -2.15836 66963 1.00008
log_lik[36] -2.42536 0.000583436 0.184499 -2.9031 -2.7866 -2.50546 -2.38443 -2.3008 -2.20712 -2.18051 100000 1.00001
log_lik[37] -2.77796 0.00108888 0.344334 -3.61273 -3.42961 -2.96609 -2.71626 -2.52376 -2.33709 -2.29548 100000 0.99997
log_lik[38] -2.56802 0.000821847 0.259891 -3.21964 -3.07226 -2.69821 -2.51079 -2.37974 -2.25574 -2.22447 100000 0.999978
log_lik[39] -2.40097 0.000587076 0.168242 -2.83704 -2.72648 -2.47175 -2.36656 -2.28945 -2.19963 -2.17254 82126 1.00004
log_lik[40] -3.7057 0.00191992 0.607131 -5.07782 -4.81287 -4.07006 -3.63623 -3.26529 -2.83716 -2.72761 100000 0.999971
log_lik[41] -2.56802 0.000821847 0.259891 -3.21964 -3.07226 -2.69821 -2.51079 -2.37974 -2.25574 -2.22447 100000 0.999978
log_lik[42] -2.64915 0.000989698 0.31297 -3.42632 -3.25123 -2.80993 -2.57989 -2.41906 -2.27508 -2.24037 100000 0.999983
log_lik[43] -2.64915 0.000989698 0.31297 -3.42632 -3.25123 -2.80993 -2.57989 -2.41906 -2.27508 -2.24037 100000 0.999983
log_lik[44] -2.48078 0.000720873 0.22796 -3.07134 -2.92787 -2.58129 -2.42439 -2.32413 -2.22039 -2.19265 100000 0.999998
log_lik[45] -2.47998 0.0007191 0.227399 -3.06829 -2.92811 -2.58068 -2.42439 -2.32295 -2.2199 -2.19207 100000 0.999991
log_lik[46] -2.37927 0.000572994 0.1578 -2.78904 -2.68228 -2.44339 -2.34819 -2.27626 -2.19009 -2.16486 75843 1.00002
log_lik[47] -3.1857 0.00158996 0.502788 -4.36319 -4.12359 -3.47549 -3.11292 -2.81382 -2.50661 -2.43336 100000 0.999983
log_lik[48] -2.34584 0.00051513 0.126888 -2.65972 -2.57614 -2.40448 -2.32682 -2.26275 -2.18111 -2.15621 60675 1.00006
log_lik[49] -2.41088 0.000514388 0.162664 -2.82099 -2.72195 -2.48676 -2.37974 -2.29996 -2.20902 -2.18182 100000 0.999999
log_lik[50] -2.76093 0.000962221 0.304281 -3.48081 -3.32983 -2.932 -2.71175 -2.53689 -2.35856 -2.31605 100000 0.999983
log_lik[51] -2.55227 0.000724273 0.229035 -3.11222 -2.98866 -2.67227 -2.5078 -2.38647 -2.26473 -2.23413 100000 0.999982
log_lik[52] -2.32991 0.000421662 0.109997 -2.58737 -2.52469 -2.38765 -2.31737 -2.25642 -2.17761 -2.15289 68051 1.00014
log_lik[53] -2.97519 0.001168 0.369353 -3.83217 -3.65693 -3.18955 -2.92343 -2.70504 -2.47099 -2.41237 100000 1.00001
log_lik[54] -2.32991 0.000421662 0.109997 -2.58737 -2.52469 -2.38765 -2.31737 -2.25642 -2.17761 -2.15289 68051 1.00014
log_lik[55] -2.41088 0.000514388 0.162664 -2.82099 -2.72195 -2.48676 -2.37974 -2.29996 -2.20902 -2.18182 100000 0.999999
log_lik[56] -2.76093 0.000962221 0.304281 -3.48081 -3.32983 -2.932 -2.71175 -2.53689 -2.35856 -2.31605 100000 0.999983
log_lik[57] -2.55227 0.000724273 0.229035 -3.11222 -2.98866 -2.67227 -2.5078 -2.38647 -2.26473 -2.23413 100000 0.999982
log_lik[58] -5.96118 0.00310773 0.982751 -8.08284 -7.69382 -6.58038 -5.88692 -5.26594 -4.47592 -4.2481 100000 1
log_lik[59] -2.52657 0.000757293 0.239477 -3.12965 -2.99042 -2.64358 -2.47232 -2.35648 -2.24149 -2.21098 100000 0.999994
log_lik[60] -2.42542 0.000581973 0.184036 -2.89803 -2.78347 -2.50676 -2.38427 -2.30052 -2.2072 -2.17959 100000 1.00003
log_lik[61] -3.70403 0.00191809 0.606552 -5.07289 -4.8054 -4.06898 -3.63492 -3.26146 -2.84054 -2.73181 100000 1
log_lik[62] -2.34145 0.000491539 0.121126 -2.63625 -2.55886 -2.40019 -2.32462 -2.26172 -2.18046 -2.15654 60724 1.00006
log_lik[63] -2.3498 0.000495617 0.128819 -2.66905 -2.58578 -2.40933 -2.33002 -2.26458 -2.18315 -2.158 67557 1.00007
log_lik[64] -2.34145 0.000491539 0.121126 -2.63625 -2.55886 -2.40019 -2.32462 -2.26172 -2.18046 -2.15654 60724 1.00006
log_lik[65] -3.05587 0.00137053 0.4334 -4.06844 -3.86188 -3.3063 -2.9924 -2.73696 -2.46901 -2.40783 100000 0.999993
log_lik[66] -2.5683 0.000819629 0.259189 -3.21226 -3.06855 -2.70014 -2.51147 -2.38054 -2.25607 -2.22414 100000 1.00001
log_lik[67] -2.49983 0.00075334 0.238227 -3.10952 -2.96599 -2.60912 -2.4412 -2.33334 -2.22697 -2.19767 100000 1
log_lik[68] -2.36999 0.00057377 0.149574 -2.75292 -2.65382 -2.43233 -2.34244 -2.27251 -2.18818 -2.16267 67957 1.00001
log_lik[69] -2.36999 0.00057377 0.149574 -2.75292 -2.65382 -2.43233 -2.34244 -2.27251 -2.18818 -2.16267 67957 1.00001
log_lik[70] -2.36999 0.00057377 0.149574 -2.75292 -2.65382 -2.43233 -2.34244 -2.27251 -2.18818 -2.16267 67957 1.00001
log_lik[71] -2.46125 0.000725389 0.215494 -3.01841 -2.88701 -2.5549 -2.4086 -2.31443 -2.21489 -2.1868 88253 0.999968
log_lik[72] -2.38928 0.000607463 0.165664 -2.81825 -2.70878 -2.45625 -2.35539 -2.28072 -2.19397 -2.16786 74373 1.00004
log_lik[73] -2.38928 0.000607463 0.165664 -2.81825 -2.70878 -2.45625 -2.35539 -2.28072 -2.19397 -2.16786 74373 1.00004
log_lik[74] -3.67851 0.0018051 0.570824 -4.96139 -4.70864 -4.02303 -3.61492 -3.26414 -2.86209 -2.75324 100000 0.999999
log_lik[75] -2.42412 0.000559944 0.17707 -2.87225 -2.76381 -2.50505 -2.38741 -2.30404 -2.20929 -2.18255 100000 1.00006
log_lik[76] -4.85816 0.00258955 0.818888 -6.64133 -6.31093 -5.36921 -4.78662 -4.26925 -3.64582 -3.46943 100000 0.999995
log_lik[77] -2.39 0.000493201 0.155964 -2.78444 -2.6866 -2.4595 -2.36016 -2.28576 -2.1976 -2.17153 100000 1
log_lik[78] -2.9613 0.00121791 0.385138 -3.86024 -3.67563 -3.1825 -2.90355 -2.6775 -2.44311 -2.38634 100000 0.999986
log_lik[79] -2.39 0.000493201 0.155964 -2.78444 -2.6866 -2.4595 -2.36016 -2.28576 -2.1976 -2.17153 100000 1
log_lik[80] -3.41139 0.00160314 0.506957 -4.55559 -4.33322 -3.71439 -3.35026 -3.03988 -2.69708 -2.60639 100000 0.999997
log_lik[81] -3.67851 0.0018051 0.570824 -4.96139 -4.70864 -4.02303 -3.61492 -3.26414 -2.86209 -2.75324 100000 0.999999
log_lik[82] -2.78322 0.00103984 0.328827 -3.56402 -3.39689 -2.96697 -2.72775 -2.54094 -2.35577 -2.31086 100000 1.00001
log_lik[83] -2.84325 0.0011636 0.367961 -3.72012 -3.53811 -3.04991 -2.78014 -2.57055 -2.36481 -2.31905 100000 1
log_lik[84] -2.36323 0.000525796 0.140284 -2.71569 -2.62599 -2.42459 -2.33932 -2.27104 -2.18727 -2.16231 71184 0.999996
log_lik[85] -2.66497 0.000953673 0.301578 -3.40396 -3.24526 -2.82653 -2.60503 -2.44195 -2.29206 -2.25465 100000 1.00001
log_lik[86] -2.37956 0.000542134 0.152967 -2.77117 -2.67185 -2.44493 -2.35106 -2.27869 -2.19225 -2.16674 79613 1
log_lik[87] -2.37956 0.000542134 0.152967 -2.77117 -2.67185 -2.44493 -2.35106 -2.27869 -2.19225 -2.16674 79613 1
log_lik[88] -2.37956 0.000542134 0.152967 -2.77117 -2.67185 -2.44493 -2.35106 -2.27869 -2.19225 -2.16674 79613 1
log_lik[89] -2.61597 0.000890215 0.281511 -3.30934 -3.15762 -2.7638 -2.55755 -2.40974 -2.27316 -2.23913 100000 1
log_lik[90] -2.37956 0.000542134 0.152967 -2.77117 -2.67185 -2.44493 -2.35106 -2.27869 -2.19225 -2.16674 79613 1
log_lik[91] -2.35169 0.000524891 0.13297 -2.6827 -2.595 -2.41106 -2.33035 -2.26492 -2.18286 -2.1578 64175 0.999993
log_lik[92] -3.39297 0.00177845 0.562394 -4.69445 -4.4338 -3.71977 -3.31813 -2.98064 -2.61791 -2.52666 100000 0.999975
log_lik[93] -2.56962 0.000873156 0.276116 -3.26589 -3.1083 -2.7045 -2.50481 -2.37047 -2.24856 -2.21752 100000 0.999975
log_lik[94] -2.73918 0.00110666 0.349957 -3.59245 -3.40636 -2.92461 -2.66867 -2.48011 -2.3068 -2.26847 100000 0.999968
log_lik[95] -2.35169 0.000524891 0.13297 -2.6827 -2.595 -2.41106 -2.33035 -2.26492 -2.18286 -2.1578 64175 0.999993
log_lik[96] -3.73209 0.00205986 0.651384 -5.20444 -4.91743 -4.1244 -3.65239 -3.25571 -2.80926 -2.69637 100000 0.999975
log_lik[97] -2.56962 0.000873156 0.276116 -3.26589 -3.1083 -2.7045 -2.50481 -2.37047 -2.24856 -2.21752 100000 0.999975

水準の差

changeup curve cutball fork slider 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

f:id:ajhjhaf:20170531233333p:plain

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

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

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

4章 対応ある2群の差と相関の推測

内容

  • 対応ある2群のt検定のオルタナとして機能します
  • 対応ある、の意味とは?
    • 同じ観察対象から2回測定しているもの
    • beforeの体重とafterの体重のセット * n個など
  • この解析をするときには実験デザインが大事になります
    • どちらかの群にバイアスがかかることを避ける事が必要です
    • 対応ある2群の実験デザインを行う
      • マッチング : 施策実行前の状態が同じ2つを組にして、ランダムに2群に割り当てる
      • プリテスト・ポストテスト : 施策の前後で同じ対象を観察する
  • 相関関係を表現する要約統計量の例
    • 共分散s
      • 平均偏差の積の平均値
      •  s = \Sigma_{i=1}^n (\frac{v_{1i} * v_{2i}}{n})
      • 正の相関がある時に正、負の相関がある時に負の値になる
      • 相関の強さは表現できない
    • 相関係数 r or  \rho
      • データの正規化を行った後に、積の平均値を計算
      • r = \Sigma_{i=1}^n\frac{z_{1i} * z_{2i}}{n}
      • -1 <= r <= 1
  • 2変量正規分布の導入
    • 共分散: \sigma_{1,2} = \sigma_1 * \sigma_2 * \rhoの関係性を持つ
    • 正規分布に従う2変量が観測される確率をモデル化する際に、あてはめが可能となる理論分布
  • 2群の差異の考察のバリエーション
    • 独立した2群の差の分析
    • 対応ある2群の群間差の分析(Inter)
    • 対応ある2群の個人内差の分析(Intra)
  • 対応ある場合の生成量(2変量に相関が無い場合は、 \rho = 0とすればokです
    • 対応が無い場合に使っていた生成量ももちろん使えます(省略)
    • 対応のある値の差(pair diff)
      •  x_{1i}^\ast - x_{2i}^\astのこと
        • 単に群同士を比較するのではなく、群間である対応関係を利用する
      •  x_{1i}^\ast - x_{2i}^\ast \sim Normal(\mu_1 - \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2})に従う
        • この導出についてはこれを参考に。
    • 対応ある値の差の効果量(pair es)
      •  \delta' = \frac{\mu_1 - \mu_2}{\sigma'}
        •  \sigma' = \sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2}のこと
      • 差の平均値に対するばらつきの平均値を評価できる
    • 対応ある値の差の優越率(pair pod)
      •  \pi_d' = p(x_{1i}^\ast - x_{2i}^\ast) > 0
        • MCMCでは、 F(es_{12}|0, 1)を評価
        • なおここで F(x|0, 1)は平均0、標準偏差1とした時の累積分布関数としています
      • 個人の変動が0より大きい確率を評価できる
    • 対応ある値の差の閾上率(pair pbt)
      •  \pi_c' = p(x_{1i}^\ast - x_{2i})^\ast > c
        • MCMCでは、 F(\frac{\mu_{diff} - c}{\sigma'}|0, 1)を評価
      • 個人の変動がcより大きい確率を評価できる
    • 同順率(pbc)
      •  Con = p((x_{1i} - x_{1j})(x_{2i} - x_{2j}) >0) = 0.5 + \frac{1}{\pi}\sin^{-1}\rho
      • 2つの i, jを選んだ時に、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

f:id:ajhjhaf:20170531011841p:plain

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

f:id:ajhjhaf:20170531011835p:plain

はじめての 統計データ分析 ―ベイズ的〈ポスト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

f:id:ajhjhaf:20170525183254p:plain

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

その4です。今回は第3章の話をしていきます。

3章 独立した2群の差の推測

内容

  • 2群に分けた群比較をする際には、ランダム化による交絡因子の排除が重要!(ランダマイゼーション)
  • 2群の、どんな差を見たいかで仮説は変わる
    • 群1の平均が群2の平均を上回る確率
    • 群1, 2の平均値の差のt年推定、区間推定
    • 群1, 2の平均値の差が基準点cより大きい確率
    • 効果量が基準点cより大きい確率
    • 非重複度、優越率、閾上率…etc
  • 2群でモデルを作る時に、標準偏差は共通させるか?分けるか?
    • 共通させる場合
      • f(\theta) = f(\mu_1, \mu_2, sigma) = f(\mu_1)f(\mu_2)f(\sigma)
      • f(\theta|x) = f(\mu_1, \mu_2, \sigma|x1, x2) ∝ f(x1, x2| u1, u2, \sigma)f(u1, u2, \sigma)
    • 独立に定義する場合
      • f(\theta) = f(\mu_1, \mu_2, \sigma) = f(\mu_1)f(\mu_2)f(\sigma1)f(\sigma2)
      • f(\theta|x) = f(\mu_1, \mu_2, \sigma|x1, x2) ∝ f(x1, x2| \mu_1, \mu_2, \sigma_1, \sigma_2)f(\mu_1, \mu_2, \sigma_1, \sigma_2)
    • かならずしも独立に定義したほうが良いわけではない。モデルの複雑性は増してしまうので。
  • 生成量の例
    • 平均の差
    • 平均の差の効果量(estimated size)
      • es = \frac{\mu_1 - \mu_2} {\sigma}
    • 非重複度(Cohen’s U) : u1が群2で何%点に相当するか
      • U_3 = F(\mu_1| \mu_2, \sigma)
      • 0.5 で2群が完全に重複していることを示す
    • 優越率(pid) : 無作為に選んだ1方の群の測定値が、無作為に選んだもう一方の測定値を上回る確率
      •  pid = p(x_1^\ast - x_2^\ast > 0)
      • 評価をするには、標準正規分布における \frac{es}{\sqrt{2}}の分布関数を調べれば良い
    • 閾上率(pic) : 無作為に選んだ一方の群の測定値と、もういっぽうの測定値の差がcより大きくなる確率
      • pic = p(x^\ast_1 - x^\ast_2 > c)
      • 評価をするには、標準正規分布における \frac{\mu_1 - \mu_2 - c}  {\sqrt{2}\sigma}の分布関数を調べればよい
  • モデルの選択

章末問題

データの取得

とりあえずデータを読み込みます。 今回は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

f:id:ajhjhaf:20170525001759p:plain

f:id:ajhjhaf:20170525010310p:plain

さてこの結果から何が言えるかというと…

  • 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

f:id:ajhjhaf:20170525001803p:plain

問題なく、パラメータ値が収束して計算されました(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検定です。