Easy to type

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

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