Easy to type

非正規労働者の職業訓練記録です。ボーナスと福利厚生を勝ち取る夢を持っています。

Stanで生存時間解析(Weibull 回帰)

生存時間解析とは?

生存時間解析は、イベントの時間を解析するための手法です。例えば、

など、様々な応用例が考えられます。生存時間解析はノンパラメトリックな手法で行なうことが多いです。具体的には、

  • カプランマイヤー推定量で生存時間を推定
  • 生存時間の違いをログランク検定で判断
  • 共変量の違いをCoxの比例ハザードモデルやAelanの加法モデルで解析

といった辺りが挙げられるでしょうか。その原因としては、生存時間のデータへ打ち切りという特徴があることが挙げられます。結果として分布にあてはめるだけの解析では、打ち切りを説明することが難しく上手くフィッティング出来ません。並べて、最尤法を使った分布への当てはめでは、値の信頼性が出ません。これに対してノンパラメトリックなカプランマイヤー法を使えば、推定値の範囲(信頼性)も出せますので論文も書きやすいです。

しかし仮に分布を仮定しても、ベイズ的にモデリングすれば、これらの問題を解決することができます。今回はStanを使って生存時間解析を行なうことで、打ち切りを含めてモデリングをしてみました。何分不勉強なもので、ご指導ご鞭撻をいただければ幸いです。

テストデータ

とりあえず解析に使うデータを用意します。今回はpythonlifeline packageに含まれる、rossiのデータを使います。データのソース論文はこちらです。このデータは逮捕された人が再逮捕されるまでの期間(週)と、その他の共変量が記されています。何はともあれ、データを見てみましょう。上から10件だけ表示します。

from lifelines.datasets import load_rossi

df = load_rossi()
df
index week arrest fin age race wexp mar paro prio
0 20 1 0 27 1 0 0 1 3
1 17 1 0 18 1 0 0 1 8
2 25 1 0 19 0 1 0 1 13
3 52 0 1 23 1 1 1 1 1
4 52 0 0 19 0 1 0 1 3
5 52 0 0 24 1 1 0 0 2
6 23 1 0 25 1 1 1 1 0
7 52 0 1 21 1 1 0 1 4
8 52 0 0 22 1 0 0 0 6
9 52 0 0 20 1 1 0 0 0
10 52 0 1 26 1 0 0 1 3

詳細は論文の3.2節に書いてありますが、今回注目するのはweekとarrestです。arrestが1である人は再逮捕されています。0である人は観測期間中では再逮捕はされず、weekの値が52になっています。所謂打ち切りの発生です。

モデル

ヒル本P115 7.8 打ち切り、に従い、現象をモデル化してみます。なお今回のデータはweek単位への丸めも発生しているのですが、本題ではないので飛ばします。

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

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

打ち切りの原理を踏まえて、今回の現象のモデルを次のように考えました。

「再逮捕の区間を表すパラメータと分布に従って、再逮捕までの時間 yが得られる。ただし 52 \lt yである場合は、 y = 52とする 」

分布はワイブル分布を利用します。ワイブル分布は次のような特徴をもち、時間のモデル化へ向いています。

  •  0 \lt t \lt \inftyにおける累積密度関数の和が1になる(イベントの発生確率を無限時間まで考えられる)
  • 現象の傾向を示す形状パラメータ mと、現象の時間スケールを示す尺度パラメータ \etaで表現される
    •  m \lt 1なら、時間とともにハザード関数値(時間 tまでにイベントが発生しなかった場合、その瞬間にイベントが起こる確率)が低くなる
    •  m = 1なら、時間に依らずハザード関数値が一定
    •  m \gt 1なら、時間にとともにハザード関数値が高くなる

ワイブル分布を使って、値に対する尤度を考えます。観測値が52で無い場合は

weibull\_lpdf(y|m, \eta)

として得られます。肝心なのは52の場合です。この時尤度は


\begin{align}
Prob[52 \lt y]   &= \int_{52}^{\infty}weibull(52|m, \eta) \\\
    &= 1 - \int_0^{52} weibull(52|m, \eta) \\\
    &= 1 - weibull\_cdf(52|m, \eta) \\\
    &= weibull\_ccdf(52|m, \eta)\end{align}

となります。

以上を踏まえて、まず共変量に依らず、全ての人が共通の現象に従っているとして、Stanを書いてみました。

data {
    int N ;
    int week[N] ;
    int arrest[N] ;
}

parameters {
    real shape ;
    real scale ;
}

model {
    for(n in 1:N){
        if(arrest[n] == 0){
            target += weibull_lccdf(week[n]| shape, scale) ;
        }else{
            target += weibull_lpdf(week[n]| shape, scale) ;
        }
    }
}

generated quantities {
    real log_lik[N] ;
    for(n in 1:N){
        if(arrest[n] == 0){
            log_lik[n] = weibull_lccdf(week[n]| shape, scale) ;
        }else{
            log_lik[n] = weibull_lpdf(week[n]| shape, scale) ;
        }
    }
}

解析と可視化

データに対して、構築したStanモデルを当てはめてみます。

import numpy as np
import os 
import pystan
import pandas as pd
import pickle
from lifelines.datasets import load_rossi


def stan_compile(stan_file_path, compiled_file_path, recompile=False):
    if os.path.exists(compiled_file_path) is False or recompile is True:
        model = pystan.StanModel(file=stan_file_path)
        with open(compiled_file_path, "wb") as f:
            pickle.dump(model, f)
    else:
        with open(compiled_file_path, "rb") as f:
            model = pickle.load(f)
    return model 

def get_summary_df(fit):
    summary = fit.summary(probs=prob)
    summary_df = pd.DataFrame(summary["summary"],
                              index=summary["summary_rownames"],
                              columns=summary["summary_colnames"])
    return summary_df

def get_waic(fit):
    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))) 
    return waic

# Stanに入れるデータの準備
df = load_rossi()
N = len(df)
week = df["week"]
arrest = df["arrest"]
stan_data = {"N": N, "week": week, "arrest": arrest}
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

# Stanをコンパイル
stan_file = os.path.join("weibull.stan")
stan_file_c = os.path.join("weibull.pkl")
model = stan_compile(stan_file, stan_file_c)

# StanへKick
fit = model.sampling(data=stan_data,
                     iter=3000,
                     chains=5,
                     warmup=1000,
                     seed=1234)

# 結果の集計
summary_df = get_summary_df(fit)
waic = get_waic(fit)
print(waic)
summary_df

output

WAIC: 1397.22843377

parameter mean se_mean sd 2.5% 5% 25% 50% 75% 95% 97.5% n_eff Rhat
shape 1.35177 0.00225134 0.123002 1.12027 1.15716 1.26817 1.347 1.43093 1.56233 1.61164 2985 1.00366
scale 127.559 0.267877 14.5396 103.759 106.833 117.305 125.893 136 153.639 161.079 2946 1.00274
log_lik[0] -5.26822 0.000938137 0.0870498 -5.44863 -5.41606 -5.3248 -5.26577 -5.20816 -5.13211 -5.10419 8610 0.999675
log_lik[1] -5.30895 0.000968553 0.0933687 -5.5044 -5.46746 -5.36939 -5.30494 -5.2457 -5.16148 -5.13464 9293 0.999847
log_lik[2] -5.21905 0.00110178 0.0852228 -5.39228 -5.36399 -5.27516 -5.21646 -5.16062 -5.08346 -5.05936 5983 1.00011

まず集計表です。全パラメータで \hat{R} \lt 1.1が確認できましたので、収束は問題ないようです。続いて、得られたパラメータ値を用いて可視化を行います。

from scipy.stats import weibull_min
import matplotlib.pyplot as plt


# 可視化その1
fig, ax1 = plt.subplots()
x = np.arange(0, 500, 1)
y = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "mean"])
y_min = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "5%"]) 
y_max = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "95%"]) 
# 確率密度関数
ax1.plot(x, y, "-")
ax1.fill_between(x, y_min, y_max, facecolor="pink")
ax1.set_ylabel("probability")
ax1.tick_params("y")
ax1.set_xlabel("Weeks")
# 件数
ax2 = ax1.twinx()
ax2.hist(df["week"])
ax2.set_ylabel("count")
ax2.tick_params("y")
plt.show()

# 可視化その2
fig, ax1 = plt.subplots()
x = np.arange(0, 50, 1)
y = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "mean"])
y_min = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "5%"]) 
y_max = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "95%"]) 
# 生存時間関数
ax1.plot(x, y, "-")
ax1.set_ylabel("Propotion Not Arrested")
ax1.fill_between(x, y_min, y_max, facecolor="pink")
ax1.set_xlabel("Weeks")
ax1.tick_params("y")

f:id:ajhjhaf:20170805173622p:plain

f:id:ajhjhaf:20170805173634p:plain

一枚目のグラフが、確率密度関数とイベント発生のカウントです。MCMCEAPが青線、scale parameterの90%信用区間を用いて描かれた範囲がピンク色になっています。

二枚目のグラフが生存時間関数(survival function)です。見方は1枚目と同様です。元論文のP7には、Coxの比例ハザードモデルで描かれた場合の結果が載っています。今回のWeibull分布を用いて描いた結果と相違ありません。

打ち切り補正をしない場合は?

当然ここで気になってくるのは、補正の効果です。先程の1枚目のグラフは52辺りに最頻値があります。打ち切りをしている効果はどの程度でしょうか?これを検証するために、打ち切りを考慮しないモデルと比較しました。

f:id:ajhjhaf:20170805175603p:plain

打ち切りを考慮しない場合は、めちゃめちゃ過剰適合しているのがわかります。打ち切りのメカニズムをモデルに取り入れることで、52以上でもイベントが生起することを踏まえてモデル化できてるわけですね。ちなみに打ち切りが考慮されない場合、WAICも3411.15613467となり、モデルとしても性能が良くないことが分かります。

比例ハザードモデル

次に共変量の効果を調べます。一番単純な比例ハザードモデルで取り組んでみましょう。比例ハザードモデルは、全時間を通してハザード間数値を一定として、共変量に従ってハザード関数の値が大きくなることを仮定するモデルです。Weibull分布のハザード関数は


\begin{align}
    h(t|m, \eta) &= \frac{m}{\eta}(\frac{x}{\eta})^{m-1} \\\
    &= \frac{m}{\eta^m}x^{m-1} \\\
    &= m \lambda x^{m-1} \\\
\end{align}

と示されます。先行研究[1][2]では、尺度パラメータ \etaへ共変量の効果を加えています。この理由の僕個人による定性的な理解としては、形状パラメータの頑強性に帰するものと考えています。形状パラメータの値はそもそもの現象に依存します。共変量ぐらいでは、中々現象の傾向までは大きく変わりません。

更に、尺度に対する効果も尺度パラメータ \etaへダイレクトに加えるのではなく、 \lambdaへ加えるほうが多いようです。即ち基底状態に対するパラメータを \beta_0、共変量の係数を \beta、共変量のダミー変数をFとして、


    \lambda = exp(\beta_0 + \beta^TF )

ということになります。この式を \etaを中心に変形すれば


    \eta = exp(-\frac{\beta_0 + \beta^TF}{m})

です。Stanで実装するときは、この式で書けばよいかと思います。

今回の解析では、逮捕後の経済的支援(financial aid)の効果を調べてみたいと思います。モデルは次の通りです。

data {
    int N ;
    int week[N] ;
    int arrest[N] ;
    int fin[N] ;
}

parameters {
    real shape ;
    real beta[2] ;
}

model {
    for(n in 1:N){
        if(arrest[n] == 0){
            target += weibull_lccdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ;
        }else{
            target += weibull_lpdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ;
        }
    }
}

generated quantities {
    real log_lik[N] ;
    for(n in 1:N){
        if(arrest[n] == 0){
            log_lik[n] = weibull_lccdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ;
        }else{
            log_lik[n] = weibull_lpdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ;
        }
    }
}

スクリプトは殆ど同様なので飛ばします。結果として次のグラフが得られました。

f:id:ajhjhaf:20170805185919p:plain

f:id:ajhjhaf:20170805185927p:plain

2枚目のグラフが分かりやすいですね。経済的支援を受けた場合は、再犯率が低くなっています。WAICも1395.3474841となって低くなり、予測力が上がっています。更に言うと、経済的支援の効果を表すパラメータ \beta_1に対する、exp(\beta_1)の95%信用区間(片側)は 0.502 \lt exp(\beta_1) \lt 0.948でした。この区間に1(即ち効果が無い)が入っていませんので、経済的支援に効果があったことが期待されます。

まとめ

Python + Stan + Weibull分布を使って、打ち切りを含めたモデル化を行い、共変量の効果の信用区間まで計算しました。…あれ、これはCoxの比例ハザードモデルで十分なのでは(死)。

ベイズ的にモデル化することの利点は自由性です。パッと思いつくものとして、次のような応用が挙げられます。

  • 加法モデルやオッズ比モデルの利用
  • 値の丸めを加味したモデル
  • 時間tに対して、ハザード関数が変化するモデル(比例ハザードモデルでは一定とする)
  • リンク関数の中身へ正規化項を追加して、スパースにする(Lasso + weibull 回帰?)
  • 単純な2群比較のコホートデザインではない場合(利用できる共変量が他にある)

ふぅー、これでやった意味が見いだせました。ぜひ皆様も、Weibull分布に親しんで下さい。

参考

書籍

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

2節「ワイブル分布」がドンピシャでした。

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

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

打ち切りデータのモデリングについて勉強させて貰いました。

Webリソース

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

アヒル本(StanとRでベイズ統計モデリング)のChapter6にPythonで取り組んでいきます。 この章は丁寧に分布を解説していくものなので、内容の復習は飛ばします。おざなりにされそうな章ですが、自分でパラメータをいじって分布からサンプリングしてみると新しい発見もありますので挑戦してみるのをおすすめします。

今回はRどころかStanも出てきません(笑)

練習問題

練習問題(1)

ベルヌーイ分布 or カテゴリカル分布に従う乱数を生成せよ…ということですが、numpyには両方とも存在しません。 二項分布からのサンプリングを一回だけ行うものがベルヌーイ分布であり、多項分布からのサンプリングを一回だけ行うものがカテゴリカル分布であるため、それで代用しました。

カテゴリカル分布についてはRの結果と対応させるために少し手間を加えます。

import numpy as np

random_bernoulli = np.random.binomial(1, 0.5, 100)

random_categorical = np.random.multinomial(1, [0.2, 0.3, 0.5], 100)
np.where(random_categorical == 1)[1]

練習問題(2)

テキストで推奨されている、ベータ分布、ディリクレ分布、ガンマ分布、2変量正規分布、コーシー分布を生成していきます。ついでに可視化します。

ベータ分布

alphaが大きいほど右側に偏り、全体的に値が大きいほど真ん中(0.5辺り)に集合していくことが分かりました。

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

alpha = 0.5
beta = 0.5
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 1 
beta = 0.5
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 1 
beta = 1 
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 5 
beta = 5 
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

g.legend()
plt.show(g)

f:id:ajhjhaf:20170531172702p:plain

ディリクレ分布

トピックモデルで使われる分布です。ベータ分布の多次元拡張なので、alphaに対する挙動はベータ分布と同じです。

今回は3次元で行いました(可視化の限界のためです…)。 可視化の方法は先行事例を参考にしています

import matplotlib.tri as tri
import numpy as np

def plot_points(X, barycentric=True, border=True, **kwargs):
    '''
    Copy from https://gist.github.com/tboggs/8778945
    '''
    corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
    triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
    if barycentric is True:
        X = X.dot(corners)
    plt.plot(X[:, 0], X[:, 1], 'k.', **kwargs)
    plt.axis('equal')
    plt.axis('off')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    if border is True:
        plt.triplot(triangle, linewidth=1)
        
fig = plt.figure()
ax = fig.add_subplot(2,2,1)
alpha=[3,3,3]
r_diri = np.random.dirichlet(alpha, size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,2)
alpha=[0.3,0.3,0.3]
r_diri = np.random.dirichlet([0.3,0.3,0.3], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,3)
alpha=[3,0.3,0.3]
r_diri = np.random.dirichlet([3,0.3,0.3], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,4)
alpha = [1,1,1]
r_diri = np.random.dirichlet([1,1,1], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)
plt.show()

f:id:ajhjhaf:20170531172713p:plain

ガンマ分布

ガンマ分布は離散値を取る分布です。numpyの場合はshapeとscaleの2つのパラメータをとるのですが、scaleは教科書の \betaの逆数になります。 というわけで

になります。ご注意下さい。

shape = 5
scale = 5
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

shape = 5
scale = 20 
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

shape = 40 
scale = 5
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

g.legend()
plt.show()

f:id:ajhjhaf:20170531174245p:plain

2変量正規分布

平均と分散にどう2変量が関わるか示す相関係数を入れ込みます。共分散行列の書き方はもっと綺麗なものがありそうですが。。

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

mu1 = 5
mu2 = 5
sigma1 = 10
sigma2 = 40
rho = 0.6
mu = np.array([mu1, mu2])
cov = np.array([[sigma1**2, sigma1*sigma2*rho], [sigma1*sigma2*rho, sigma2**2]])
r_mnorm = np.random.multivariate_normal(mu, cov, size=1000)

g = sns.jointplot(r_mnorm[:, 0], r_mnorm[:, 1], kind="kde", n_levels=50)
plt.show(g)

f:id:ajhjhaf:20170531172724p:plain

コーシー分布

numpyではパラメータを選べない標準コーシー分布しか使えません。そのためscipyを利用しました。

なおここまでのnumpyで生成した他の分布も、rvs(random variates sampling?)メソッドを使えば同様に得ることが出来ます。 ただnumpyとscipyで仕様が異なり、numpyの正規分布を由来とする乱数はnumpy.random.normalなのに、scipyはscipy.stats.normだったりします(解せぬ)。

コーシー分布は面白い分布ですね。パラメータがあるものの、平均値や分散をそこから推定できない。しかもめちゃめちゃ変なところに外れ値みたいに値が登場する性質を持っています。

import seaborn as sns
import matplotlib.pyploy as plt
from scipy.stats import cauchy

r_cauchy = cauchy.rvs(0, 1, 100)
sns.distplot(r_cauchy)
plt.show()

f:id:ajhjhaf:20170531172819p:plain

練習問題(3)

 x_1 - x_2 \sim Normal(\mu_1 - \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2})のような分布になります。 これは以前に触れた部分ですね。

正規分布の線形結合で表される分布は、記事内のpdfを読めば全部同じように表されることが分かります。

y1 = np.random.normal(50, 20, 2000)
y2 = np.random.normal(20, 15, 2000)
diff = y1 - y2
sns.distplot(diff)
plt.show()

f:id:ajhjhaf:20170531172937p:plain

練習問題(4)

何か自分でテキストに載っていない分布を対象にすること。という問題でした。

今回は負の二項分布を取り扱ってみます。生物学の分野では、遺伝子発現量(RNA-seq or cDNA-seq)の推定モデル(DESeqedgeR)でよく使われている分布です。負の二項分布はポアソン分布とガンマ分布の混合モデルです。ホクソエムさんの分かりやすい解説があるので、そちらをご覧ください。ポアソン分布の \gammaが確率的に変動することで、過分散に対応できるようになっています。

負の二項分布の解説は英語版Wikipediaを見るのが一番いいと思います。以下そこからのざっくりとした抽出です。

「負の」の意味

負の二項分布は定義4パターン存在します。4つはそれぞれ

  • 成功率pの試行をr回失敗(成功)するまでに、成功(失敗)した回数kの分布
  • 成功率pの試行をr回失敗した時の、試行回数nの分布
  • 成功率pの試行をk回成功した時の、失敗回数rの分布
  • 成功率pの試行をk回成功した時の、試行回数nの分布

となります。4つのパターンがありますが、一番良く使われるのは1つ目のものです。

1つ目のパターンを式で表すと、

f:id:ajhjhaf:20170603170330p:plain:w300

となります。このパターンは試行が成功でも失敗でも式に大きな変化がなく、逆で紹介されている場合もあります。その時は逆側の事象が起こる確率が1-pになるように対応が変化するだけです。ところで、成功回数と失敗回数の和が試行回数です。つまり n = k + rとなります。これを使って整理すると

f:id:ajhjhaf:20170603171743p:plain:w300

となりました。ところで普通の二項分布の式は次のようになります。

f:id:ajhjhaf:20170603171236p:plain:w300

つまり負の二項分布だと二項分布と比べてトライアル nの数が-1されていることになります。

平均と分散

負の二項分布の平均と分散は次のようになります。

 \mu = \frac{pr}{1-p}

 var = \mu(1+\frac{\mu}{r})

詳しい導出についてはこちらを参考にするのが適切だと思いました。日本語版Wikipediaの記事は、どこか間違っているような…。

実数への拡張

負の二項分布のパラメータrは正の実数へ拡張することが出来ます。「rが実数の『成功率pの試行をr回失敗(成功)するまでに、成功(失敗)した回数kの分布』ってなんやねん」と僕は思いました。おそらくここで肝心なのは、rを正の実数に拡張しても元々の定義へ矛盾せずに一般化できるということではないかと思います。この時ガンマ関数 \Gamma()を利用して、

 NegativeBinomial(k|n, p) \sim \frac{\Gamma(n)}{k!\Gamma(n-k)}

のように書くことが出来ます。

平均と形状パラメータによる分布の指定

負の二項分布のパラメータを指定する時にn, pではなく、平均パラメータ\muと形状パラメータ \alphaを使って記述されることがあります。形状パラメータとは分散をスッキリさせるためにrの逆数を取って\alpha = \frac{1}{r}としたものです。一応形状パラメータと書きましたが、他にも呼び方が色々あり、

  • shape parameter
  • dispersion parameter
  • clustering coefficient
  • heterogeneity
  • aggregation

のようにとりとめがありません。文字の使い方も適当らしく、によっては \phiが使われているものもありました。

分布の可視化

可視化をする時にbinsを自分で指定しています。何も指定しないと最適なbinsをseabornが決めてくれるのですが、範囲が小数で分割された場合には離散値である負の二項分布の値がダブルカウントされるのか、物凄く変な分布になってしまいます。気をつけて下さい。

import numpy as np
from scipy.stats import nbinom
import seaborn as sns
import matplotlib.pyplot as plt

# r_nb = np.random.negative_binomial(10, 0.3, 10000)
r_nb = nbinom.rvs(10, 0.3, size=10000)

sns.distplot(r_nb, kde=False, bins=np.arange(0,100,5))
plt.show()

f:id:ajhjhaf:20170531193122p:plain


余談ですがedgeRの論文は事実上、RNA-seqの解析に負の二項分布をあてはめだだけの論文です。2ページしかありません。そのわりに1000回を超えれば場外ホームランとなる世界で、2500回も引用されている(Scopus調べ)ので、なんとも夢があります。

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