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

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

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

内容

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