はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その5
その5です。今回は第4章の章末問題に取り組んでいきます。
4章 対応ある2群の差と相関の推測
内容
- 対応ある2群のt検定のオルタナとして機能します
- 対応ある、の意味とは?
- 同じ観察対象から2回測定しているもの
- beforeの体重とafterの体重のセット * n個など
- この解析をするときには実験デザインが大事になります
- どちらかの群にバイアスがかかることを避ける事が必要です
- 対応ある2群の実験デザインを行う
- マッチング : 施策実行前の状態が同じ2つを組にして、ランダムに2群に割り当てる
- プリテスト・ポストテスト : 施策の前後で同じ対象を観察する
- 相関関係を表現する要約統計量の例
- 共分散
- 平均偏差の積の平均値
- 正の相関がある時に正、負の相関がある時に負の値になる
- 相関の強さは表現できない
- 相関係数 or
- データの正規化を行った後に、積の平均値を計算
- -1 <= r <= 1
- 共分散
- 2変量正規分布の導入
- 共分散: の関係性を持つ
- 正規分布に従う2変量が観測される確率をモデル化する際に、あてはめが可能となる理論分布
- 2群の差異の考察のバリエーション
- 独立した2群の差の分析
- 対応ある2群の群間差の分析(Inter)
- 対応ある2群の個人内差の分析(Intra)
- 対応ある場合の生成量(2変量に相関が無い場合は、とすればokです
章末問題
データの取得
取り敢えずデータを準備します。
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 |
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 |
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その4.1
対応のあるt検定を行う前に、3章の内容を等分散モデルで書くとどうなるか、をちゃんと検証します。
番外編 3章の等分散モデルと変分ベイズによる推定
等分散モデルによる差の推定
といっても、Stanのモデルをちょっと変えるだけなのでした。
まずコンパイル
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする # 等分散モデル stan_file = os.path.join("stan", "g2_ind_equ.stan") stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f)
Stanのモデルはこんな感じです。
data { int<lower=0> n_a ; int<lower=0> n_b ; real<lower=0> a[n_a] ; real<lower=0> b[n_b] ; real<lower=0> c_mu_diff ; real<lower=0> c_es ; real<lower=0> c_cohenu ; real<lower=0> c_pod ; real<lower=0> c_pbt ; real<lower=0> cdash_pbt ; } parameters { real<lower=0> mu_a ; real<lower=0> sigma ; real<lower=0> mu_b ; } model { a ~ normal(mu_a, sigma) ; b ~ normal(mu_b, sigma) ; } generated quantities { vector[n_a] log_lik ; real mu_diff ; real es ; real cohenu ; real pod ; real pbt ; int<lower=0, upper=1> prob_mu_diff_upper_0 ; int<lower=0, upper=1> prob_mu_diff_upper_c ; int<lower=0, upper=1> prob_es_upper_c ; int<lower=0, upper=1> prob_cohenu_upper_c ; int<lower=0, upper=1> prob_pod_upper_c ; int<lower=0, upper=1> prob_pbt_upper_cdash ; for(i in 1:n_a){ log_lik[i] = normal_lpdf(a[i] | mu_a, sigma) + normal_lpdf(b[i] | mu_b, sigma) ; } mu_diff = mu_a - mu_b ; es = mu_diff / sigma ; cohenu = normal_cdf(mu_a, mu_b, sigma) ; pod = normal_cdf(es / sqrt(2), 0, 1) ; pbt = normal_cdf((mu_diff - c_pbt) / ( sqrt(2) * sigma ), 0, 1) ; prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ; prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ; prob_es_upper_c = es > c_es ? 1 : 0 ; prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ; prob_pod_upper_c = pod > c_pod ? 1 : 0 ; prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ; }
後はサンプリングを同様に行うだけ。
import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # 等分散モデル # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 pars = ["mu_a", "sigma", "mu_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング fit = model.sampling(data=stan_data, pars=pars, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=pars, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) # 事後分布の可視化 for par in summary_df.index[summary_df["sd"] == 0]: pars.remove(par) fit.traceplot(pars) plt.show() # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
これを実行すると、s.d.についてはEAPが8.63と推定されます。a群とb群のちょうど中間ぐらいですね。 両方の分布に同じ分散を使っているので、推定される値としては納得できるものです。 ただ尤度については716と不等分散モデルより高くなってしまいました。よってデータにはフィッティングしていないと捉えます。
変分ベイズによる推定
これまでは事後分布の推定を、NUTSアルゴリズムによるサンプリングで行ってきました。 ところでStanには、MCMCサンプリング以外にも変分ベイズによる分布の推定も一応実装されています。ついでにこれも試してみましょう。モデルのコンパイル部分までは、上述のものと共通です。
# サンプリングをADVIでやってみる import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display from collections import OrderedDict from pystan.external.pymc import plots %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # 等分散モデル # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 pars = ["mu_a", "sigma", "mu_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind_equ.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでADVIサンプリング fit = model.vb(data=stan_data, pars=pars, iter=100000, seed=1234) # 事後分布の表を取得 # ADVIはAPIがまだ完成されていないので、summaryの表を作る方式が違う # chainを利用したサンプリングもしていないため、Rhatも計算できない vb_sample = pd.read_csv(fit["args"]["sample_file"].decode("utf-8"), comment="#") vb_sample = vb_sample.drop("lp__", 1) summary_df = vb_sample.describe(percentiles=prob).T display(summary_df) # カーネル密度推定出来ないパラメータの削除 for par in summary_df.index[summary_df["std"] == 0]: pars.remove(par) ## 事後分布の可視化 od = OrderedDict() for i, par in enumerate(fit["sampler_param_names"]): par_s = par.split(".") if len(par_s) == 1: od[par] = np.array(fit["sampler_params"][i]) else: par = par_s[0] if par in od.keys(): od[par] = np.vstack([od[par], np.array(fit["sampler_params"][i])]) else: od[par] = np.array(fit["sampler_params"][i]) plots.traceplot(od, pars) plt.show() # WAICの計算 log_lik = od["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
注意すべきところは、PyStanの変分ベイズのAPIが完成されていないためsamplingとやり方がぜんぜん違う点です。 サンプリング結果の集計にはpandasのdescribeメソッドを使っています。
結果は次のとおりです。
count | mean | std | min | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | max | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu_a | 1001 | 56.4035 | 1.14076 | 52.6137 | 54.1746 | 54.6317 | 55.6578 | 56.3594 | 57.2304 | 58.3143 | 58.6232 | 59.8892 |
sigma | 1001 | 7.8148 | 0.593153 | 5.79781 | 6.71698 | 6.90312 | 7.43062 | 7.78376 | 8.18896 | 8.82866 | 9.04493 | 10.1979 |
mu_b | 1001 | 40.3831 | 1.19016 | 36.5435 | 38.1027 | 38.3704 | 39.5665 | 40.4121 | 41.1783 | 42.2879 | 42.7135 | 44.4097 |
log_lik.1 | 1001 | -6.42276 | 0.16823 | -7.06181 | -6.76429 | -6.6899 | -6.53466 | -6.42269 | -6.31037 | -6.15472 | -6.10705 | -5.86712 |
log_lik.2 | 1001 | -6.07869 | 0.150718 | -6.67179 | -6.3649 | -6.31169 | -6.17576 | -6.08422 | -5.98092 | -5.81884 | -5.78335 | -5.58263 |
log_lik.3 | 1001 | -8.86784 | 0.466876 | -10.957 | -9.91909 | -9.69209 | -9.1477 | -8.81625 | -8.52498 | -8.19867 | -8.08506 | -7.54151 |
log_lik.4 | 1001 | -6.23111 | 0.153535 | -6.66425 | -6.54282 | -6.48627 | -6.33461 | -6.2304 | -6.12679 | -5.96821 | -5.93276 | -5.71942 |
log_lik.5 | 1001 | -6.0363 | 0.152214 | -6.58138 | -6.32926 | -6.28372 | -6.13988 | -6.03141 | -5.93468 | -5.77738 | -5.73873 | -5.49402 |
log_lik.6 | 1001 | -9.86775 | 0.610265 | -12.5208 | -11.2156 | -11.0038 | -10.2397 | -9.81467 | -9.41468 | -9.01542 | -8.86263 | -8.08728 |
log_lik.7 | 1001 | -6.96474 | 0.214733 | -7.70444 | -7.41328 | -7.32943 | -7.10373 | -6.95857 | -6.81509 | -6.63545 | -6.57965 | -6.26646 |
log_lik.8 | 1001 | -7.96041 | 0.353478 | -10.3051 | -8.72111 | -8.56551 | -8.16825 | -7.93375 | -7.71524 | -7.44256 | -7.35968 | -7.04291 |
log_lik.9 | 1001 | -6.55416 | 0.179014 | -7.26925 | -6.9163 | -6.838 | -6.66726 | -6.55088 | -6.42996 | -6.27155 | -6.21939 | -6.00683 |
log_lik.10 | 1001 | -6.20992 | 0.155187 | -6.83945 | -6.50715 | -6.45275 | -6.31807 | -6.21339 | -6.10742 | -5.94728 | -5.89824 | -5.662 |
log_lik.11 | 1001 | -7.74777 | 0.313827 | -9.05947 | -8.43614 | -8.28447 | -7.94685 | -7.71843 | -7.52812 | -7.30075 | -7.21452 | -6.86988 |
log_lik.12 | 1001 | -6.23652 | 0.159109 | -6.74122 | -6.5479 | -6.49699 | -6.34863 | -6.22906 | -6.12679 | -5.9821 | -5.94032 | -5.73857 |
log_lik.13 | 1001 | -6.24006 | 0.160162 | -6.76222 | -6.54558 | -6.50324 | -6.35045 | -6.23568 | -6.134 | -5.98247 | -5.93949 | -5.70956 |
log_lik.14 | 1001 | -6.3109 | 0.159837 | -6.75209 | -6.61471 | -6.57248 | -6.42249 | -6.30973 | -6.19668 | -6.06357 | -6.02272 | -5.77449 |
log_lik.15 | 1001 | -6.76248 | 0.194958 | -7.7655 | -7.14874 | -7.08354 | -6.89143 | -6.75921 | -6.62693 | -6.45886 | -6.40091 | -6.21029 |
log_lik.16 | 1001 | -9.50787 | 0.558317 | -12.0117 | -10.7875 | -10.4896 | -9.84995 | -9.46429 | -9.09205 | -8.71395 | -8.58059 | -7.86643 |
log_lik.17 | 1001 | -6.23197 | 0.159713 | -6.77176 | -6.55355 | -6.51252 | -6.33496 | -6.22916 | -6.12582 | -5.97284 | -5.91782 | -5.47086 |
log_lik.18 | 1001 | -10.937 | 0.792785 | -16.4895 | -12.6988 | -12.3467 | -11.3937 | -10.8907 | -10.4062 | -9.72509 | -9.5612 | -8.87464 |
log_lik.19 | 1001 | -6.55719 | 0.178717 | -7.10137 | -6.92829 | -6.86605 | -6.66928 | -6.54786 | -6.43631 | -6.2802 | -6.22669 | -5.77744 |
log_lik.20 | 1001 | -6.92278 | 0.21248 | -7.83707 | -7.3648 | -7.28035 | -7.0519 | -6.91023 | -6.78217 | -6.60749 | -6.55017 | -6.33407 |
log_lik.21 | 1001 | -7.55595 | 0.295133 | -9.58726 | -8.17106 | -8.06674 | -7.72853 | -7.5413 | -7.3452 | -7.10152 | -7.04768 | -6.85252 |
log_lik.22 | 1001 | -7.55767 | 0.281643 | -9.09015 | -8.19547 | -8.05543 | -7.72132 | -7.53802 | -7.36901 | -7.14595 | -7.075 | -6.81468 |
log_lik.23 | 1001 | -7.35302 | 0.270502 | -8.89604 | -7.91456 | -7.82611 | -7.51521 | -7.32474 | -7.17505 | -6.95945 | -6.88067 | -6.60574 |
log_lik.24 | 1001 | -6.94368 | 0.215861 | -8.22401 | -7.36714 | -7.31273 | -7.08547 | -6.93466 | -6.79251 | -6.616 | -6.54842 | -6.35364 |
log_lik.25 | 1001 | -6.71402 | 0.191982 | -7.69412 | -7.09788 | -7.03066 | -6.84066 | -6.71253 | -6.58261 | -6.40948 | -6.34514 | -6.1755 |
log_lik.26 | 1001 | -9.53694 | 0.584227 | -13.4754 | -10.843 | -10.567 | -9.86342 | -9.49629 | -9.13708 | -8.68395 | -8.54657 | -8.1407 |
log_lik.27 | 1001 | -6.78402 | 0.201555 | -7.52225 | -7.20148 | -7.13041 | -6.91297 | -6.78015 | -6.64539 | -6.46899 | -6.42453 | -6.18105 |
log_lik.28 | 1001 | -6.44432 | 0.167778 | -6.96219 | -6.78804 | -6.72768 | -6.55666 | -6.44535 | -6.32604 | -6.18597 | -6.13832 | -5.85676 |
log_lik.29 | 1001 | -7.69239 | 0.315538 | -9.81127 | -8.36852 | -8.23455 | -7.87711 | -7.67732 | -7.4738 | -7.2188 | -7.15845 | -6.87702 |
log_lik.30 | 1001 | -7.7927 | 0.319736 | -9.13938 | -8.50667 | -8.34904 | -7.99599 | -7.76069 | -7.5676 | -7.33473 | -7.24417 | -6.90979 |
log_lik.31 | 1001 | -6.16013 | 0.157118 | -6.71291 | -6.45853 | -6.41783 | -6.26426 | -6.15737 | -6.05708 | -5.90701 | -5.86343 | -5.49712 |
log_lik.32 | 1001 | -6.13017 | 0.151897 | -6.74355 | -6.41578 | -6.37076 | -6.23383 | -6.13207 | -6.0305 | -5.87105 | -5.83521 | -5.59475 |
log_lik.33 | 1001 | -7.00073 | 0.223585 | -7.81207 | -7.48043 | -7.40014 | -7.14758 | -6.99534 | -6.84432 | -6.65199 | -6.60361 | -6.3401 |
log_lik.34 | 1001 | -6.92585 | 0.215031 | -7.65724 | -7.38092 | -7.30914 | -7.06934 | -6.91029 | -6.78167 | -6.59183 | -6.54201 | -6.25093 |
log_lik.35 | 1001 | -6.78402 | 0.201555 | -7.52225 | -7.20148 | -7.13041 | -6.91297 | -6.78015 | -6.64539 | -6.46899 | -6.42453 | -6.18105 |
log_lik.36 | 1001 | -6.45221 | 0.172966 | -7.0252 | -6.79691 | -6.752 | -6.56508 | -6.448 | -6.33613 | -6.17511 | -6.12069 | -5.69448 |
log_lik.37 | 1001 | -6.43874 | 0.169585 | -6.94675 | -6.7849 | -6.73926 | -6.54679 | -6.43412 | -6.32283 | -6.16711 | -6.11841 | -5.65075 |
log_lik.38 | 1001 | -5.98581 | 0.149329 | -6.52304 | -6.27182 | -6.22623 | -6.0843 | -5.98738 | -5.88368 | -5.74347 | -5.69836 | -5.48248 |
log_lik.39 | 1001 | -6.22339 | 0.155309 | -6.81162 | -6.52915 | -6.47324 | -6.33043 | -6.22686 | -6.12094 | -5.9492 | -5.9203 | -5.72975 |
log_lik.40 | 1001 | -6.38489 | 0.163779 | -6.94275 | -6.72587 | -6.64711 | -6.49179 | -6.38321 | -6.27196 | -6.11557 | -6.0736 | -5.90043 |
log_lik.41 | 1001 | -6.66083 | 0.18857 | -7.36905 | -7.064 | -6.97551 | -6.77774 | -6.65256 | -6.5262 | -6.36894 | -6.30808 | -6.14803 |
log_lik.42 | 1001 | -6.03561 | 0.149328 | -6.59567 | -6.32253 | -6.27471 | -6.13574 | -6.03972 | -5.93328 | -5.79028 | -5.74892 | -5.56119 |
log_lik.43 | 1001 | -7.9404 | 0.350003 | -10.3391 | -8.68739 | -8.53891 | -8.1428 | -7.92091 | -7.69262 | -7.41976 | -7.33845 | -7.05618 |
log_lik.44 | 1001 | -11.5909 | 0.862142 | -15.4103 | -13.5076 | -13.1583 | -12.1299 | -11.5137 | -10.9734 | -10.3584 | -10.1675 | -9.0104 |
log_lik.45 | 1001 | -7.77097 | 0.315238 | -9.01645 | -8.47879 | -8.3403 | -7.97236 | -7.74171 | -7.5439 | -7.3051 | -7.22634 | -6.98188 |
log_lik.46 | 1001 | -7.6209 | 0.303068 | -9.7381 | -8.25764 | -8.14321 | -7.80149 | -7.6036 | -7.40669 | -7.16314 | -7.0904 | -6.8879 |
log_lik.47 | 1001 | -6.97629 | 0.217257 | -7.94213 | -7.42008 | -7.33636 | -7.10676 | -6.9649 | -6.82924 | -6.64425 | -6.58657 | -6.38334 |
log_lik.48 | 1001 | -6.66957 | 0.186988 | -7.39538 | -7.06688 | -6.98032 | -6.78594 | -6.65865 | -6.54712 | -6.38723 | -6.3303 | -6.12713 |
log_lik.49 | 1001 | -6.63225 | 0.187684 | -7.31207 | -7.0097 | -6.94733 | -6.75122 | -6.63039 | -6.50275 | -6.33472 | -6.2794 | -6.0357 |
log_lik.50 | 1001 | -6.32385 | 0.159621 | -6.76867 | -6.63689 | -6.58101 | -6.43788 | -6.32294 | -6.2061 | -6.07091 | -6.03085 | -5.80358 |
mu_diff | 1001 | 16.0204 | 1.66137 | 10.9645 | 12.6538 | 13.2746 | 14.9254 | 16.0157 | 17.0902 | 18.7446 | 19.2714 | 20.8654 |
es | 1001 | 2.06108 | 0.259806 | 1.27133 | 1.58195 | 1.65403 | 1.87671 | 2.05535 | 2.23764 | 2.50222 | 2.59257 | 2.90178 |
cohenu | 1001 | 0.977041 | 0.0139309 | 0.898194 | 0.94317 | 0.95094 | 0.969721 | 0.980078 | 0.987378 | 0.993829 | 0.995237 | 0.998145 |
pod | 1001 | 0.924156 | 0.0256932 | 0.815664 | 0.868346 | 0.878915 | 0.907751 | 0.926937 | 0.943204 | 0.961581 | 0.966615 | 0.979909 |
pbt | 1001 | 0.705715 | 0.0529771 | 0.538187 | 0.593225 | 0.61623 | 0.671596 | 0.708586 | 0.741193 | 0.786155 | 0.803224 | 0.849742 |
prob_mu_diff_upper_0 | 1001 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_mu_diff_upper_c | 1001 | 0.879121 | 0.32615 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_es_upper_c | 1001 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
prob_cohenu_upper_c | 1001 | 0.951049 | 0.215874 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
prob_pod_upper_c | 1001 | 0.152847 | 0.36002 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
prob_pbt_upper_cdash | 1001 | 0.97003 | 0.17059 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その4
その4です。今回は第3章の話をしていきます。
3章 独立した2群の差の推測
内容
- 2群に分けた群比較をする際には、ランダム化による交絡因子の排除が重要!(ランダマイゼーション)
- 2群の、どんな差を見たいかで仮説は変わる
- 群1の平均が群2の平均を上回る確率
- 群1, 2の平均値の差のt年推定、区間推定
- 群1, 2の平均値の差が基準点cより大きい確率
- 効果量が基準点cより大きい確率
- 非重複度、優越率、閾上率…etc
- 2群でモデルを作る時に、標準偏差は共通させるか?分けるか?
- 共通させる場合
- 独立に定義する場合
- かならずしも独立に定義したほうが良いわけではない。モデルの複雑性は増してしまうので。
- 共通させる場合
- 生成量の例
- モデルの選択
章末問題
データの取得
とりあえずデータを読み込みます。 今回は2群の対応のない場合の比較です。
import numpy as np from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger # データの準備 a = np.array([56,55,55,62,54,63,47,58,56,56,57,52,53,50,50,57,57,55,60,65,53,43,60,51,52, 60,54,49,56,54,55,57,53,58,54,57,60,57,53,61,60,58,56,52,62,52,66,63,54,50]) b = np.array([33,37,59,41,42,61,46,25,32,35,55,44,45,41,33,61,46,16,48,34,27,37,28,31,32, 20,50,42,26,55,45,36,51,51,50,48,47,39,36,35,32,38,25,66,54,27,35,34,49,39]) logger = get_logger()
各種統計量の計算と可視化
- 各種統計量はpandas.describeを使えば大体出せますね。
- 分布はdistplot以外にも色々ありますが、matplotlibよりseabornの方が簡単に色々でるので使いました。
import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from IPython.core.display import display %matplotlib inline df = pd.DataFrame({"a":a, "b":b}) display(df.describe()) # 分散は出ない sns.distplot(a) plt.show() sns.distplot(b) plt.show()
- 結果
a | b | |
---|---|---|
count | 50 | 50 |
mean | 55.76 | 40.38 |
std | 4.57839 | 11.1427 |
min | 43 | 16 |
25% | 53 | 33 |
50% | 56 | 39 |
75% | 58 | 48 |
max | 66 | 66 |
さてこの結果から何が言えるかというと…
- a群の方が平均値が15ぐらい高くなっている
- 中央値となる50%点もa群の方が高いが、なんといっても最小値が大幅に違う
- この2つの群はともに正規分布っぽい感じにはなっている
- 分散は2群で結構違いそう。
みたいなものがわかります。よくある統計でしたら、t検定で2群の間には有意差があるということを出して終わりですね。 今回は最初に書いた内容を使って、
- 差の大きさがどれぐらい幅があるのか(ちょっとだけの差でしたら、殆ど意味がないですよね)。
- 差は14ぐらいありそうである。それはどれぐらい確かなのか。
- 2つの分布はどれぐらい離れているのか。重複してる部分はどれぐらいあるのか。
- 群としてはaの方が大きくても、この群から適当に1つずつ選んだ時に、どれぐらいの確立でaの方の個体が大きな値になるのか
と言ったことを調べていくわけです。
Stanによる統計モデルの構築
Stanでモデルを作って、2群を比較します。 教科書に従って、とりあえず不偏分散のモデルを作っています。aとbの分散値が先程の要約テーブルで結構違いましたので、等分散よりモデルとしては適当になりそうです(後で検証します)
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "g2_ind.stan") stan_file_c = os.path.join("stan", "g2_ind.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f)
# Stanのモデル data { int<lower=0> n_a ; int<lower=0> n_b ; real<lower=0> a[n_a] ; real<lower=0> b[n_b] ; real<lower=0> c_mu_diff ; real<lower=0> c_es ; real<lower=0> c_cohenu ; real<lower=0> c_pod ; real<lower=0> c_pbt ; real<lower=0> cdash_pbt ; } parameters { real<lower=0> mu_a ; real<lower=0> sigma_a ; real<lower=0> mu_b ; real<lower=0> sigma_b ; } model { a ~ normal(mu_a, sigma_a) ; b ~ normal(mu_b, sigma_b) ; } generated quantities { vector[n_a] log_lik ; real mu_diff ; real es ; real cohenu ; real pod ; real pbt ; int<lower=0, upper=1> prob_mu_diff_upper_0 ; int<lower=0, upper=1> prob_mu_diff_upper_c ; int<lower=0, upper=1> prob_es_upper_c ; int<lower=0, upper=1> prob_cohenu_upper_c ; int<lower=0, upper=1> prob_pod_upper_c ; int<lower=0, upper=1> prob_pbt_upper_cdash ; for(i in 1:n_a){ log_lik[i] = normal_lpdf(a[i] | mu_a, sigma_a) + normal_lpdf(b[i] | mu_b, sigma_b) ; } mu_diff = mu_a - mu_b ; es = mu_diff / sigma_a ; cohenu = normal_cdf(mu_a, mu_b, sigma_b) ; pod = normal_cdf(mu_diff / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; pbt = normal_cdf((mu_diff - c_pbt) / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ; prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ; prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ; prob_es_upper_c = es > c_es ? 1 : 0 ; prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ; prob_pod_upper_c = pod > c_pod ? 1 : 0 ; prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ; }
統計モデルによる事後分布のサンプリング
注意点は2つあります。
まず1つは可視化をする部分です。
PyStanのモジュールをそのまま使ったサンプリング値の可視化では、サンプリング値が全て同値である場合(=分散が0の時)はtraceplotでエラーが発生してしまいます。
これは可視化するために行っている、カーネル密度推定の行列計算で失敗するためです。
今回、僕は自分でPyStanの一部を改造したものを利用しました。
が、これは試してみたかっただけなので、こんな面倒なことをせず、summary_df
のsd columnを見て、0だったらtrace plotを書かなければ良いと思います。
もう1つはWAICの計算の部分です。 計算をする際に、データ数分だけfor文を回して対数尤度を計算しています この実装をしているものが検索してみると見つかるし、豊田本上級でもこの紹介がされています(定義的にも正しい?)。 ただ今回の場合は、モデルが単純で50個のデータが共通の母数を持つ正規分布から生成されていので、対数尤度は1次元配列で計算し、waicの計算でsumをする部分が無くても大丈夫です。 本の付録のモデルでは、こちらの簡潔版で計算されています。
import pandas as pd import pickle import pystan import matplotlib import os import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 50) # Stanで使うデータの用意 stan_data = {"n_a": a.size, "n_b": b.size, "a": a, "b": b, "c_mu_diff": 14, #標本平均の差を使ってみる "c_es": 3.0, # 標本効果量で、a群からみたもの "c_cohenu": 0.95, # a群から見た非重複度 "c_pod": 0.95, # 優越率 "c_pbt": 10, # 閾上率の基準値 "cdash_pbt": 0.60} #閾上率 # 興味のあるパラメータの設定 par = ["mu_a", "sigma_a", "mu_b", "sigma_b", "log_lik", "mu_diff", "es", "cohenu", "pod", "pbt", "prob_mu_diff_upper_0", "prob_mu_diff_upper_c", "prob_es_upper_c", "prob_cohenu_upper_c", "prob_pod_upper_c", "prob_pbt_upper_cdash"] prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # モデルの読み込み stan_file_c = os.path.join("stan", "g2_ind.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="NUTS") # 事後分布の表を取得 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) # 事後分布の可視化 fit.traceplot(par, {"prob_mu_diff_upper_0":np.int}) plt.show() # WAICの計算 log_lik = fit.extract("log_lik")["log_lik"] waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0)) logger.info("WAICの値は{0}です".format(waic))
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu_a | 55.7604 | 0.00210572 | 0.665888 | 54.4483 | 54.6698 | 55.3161 | 55.7601 | 56.2025 | 56.8535 | 57.0771 | 100000 | 0.99999 |
sigma_a | 4.69936 | 0.00163133 | 0.492013 | 3.85845 | 3.97235 | 4.35244 | 4.65696 | 5.00083 | 5.56823 | 5.78077 | 90964 | 0.999993 |
mu_b | 40.3768 | 0.00513303 | 1.62321 | 37.1566 | 37.7085 | 39.2942 | 40.3763 | 41.4601 | 43.0385 | 43.5743 | 100000 | 0.999979 |
sigma_b | 11.4351 | 0.00377584 | 1.19403 | 9.39226 | 9.66662 | 10.5891 | 11.3316 | 12.1638 | 13.5548 | 14.0673 | 100000 | 1.00003 |
log_lik[0] | -6.04735 | 0.000490513 | 0.152051 | -6.35677 | -6.30426 | -6.14733 | -6.04291 | -5.94225 | -5.80452 | -5.76237 | 96090 | 0.999988 |
log_lik[1] | -5.88977 | 0.000482951 | 0.146824 | -6.18808 | -6.13823 | -5.98677 | -5.88593 | -5.78838 | -5.65612 | -5.61159 | 92424 | 1.00001 |
log_lik[2] | -7.21332 | 0.000987215 | 0.312185 | -7.8894 | -7.76195 | -7.40848 | -7.19029 | -6.99224 | -6.74368 | -6.67145 | 100000 | 0.999982 |
log_lik[3] | -6.74241 | 0.000737068 | 0.233082 | -7.24164 | -7.15023 | -6.89021 | -6.72715 | -6.57867 | -6.38654 | -6.32838 | 100000 | 0.999992 |
log_lik[4] | -5.91407 | 0.000486231 | 0.147301 | -6.21372 | -6.16309 | -6.01166 | -5.90981 | -5.8126 | -5.67868 | -5.63577 | 91775 | 1.00002 |
log_lik[5] | -8.73421 | 0.00139258 | 0.440372 | -9.66647 | -9.49826 | -9.01462 | -8.7067 | -8.42342 | -8.06105 | -7.95471 | 100000 | 0.999977 |
log_lik[6] | -7.74988 | 0.00123701 | 0.391176 | -8.59949 | -8.44025 | -7.99428 | -7.7188 | -7.47338 | -7.16121 | -7.07202 | 100000 | 1.00001 |
log_lik[7] | -6.88161 | 0.000752406 | 0.237932 | -7.39441 | -7.29783 | -7.03121 | -6.86651 | -6.71561 | -6.51828 | -6.45829 | 100000 | 0.999984 |
log_lik[8] | -6.10953 | 0.000492635 | 0.155785 | -6.42726 | -6.37336 | -6.21194 | -6.10435 | -6.00181 | -5.86187 | -5.81865 | 100000 | 0.999984 |
log_lik[9] | -5.94667 | 0.00048357 | 0.148166 | -6.24786 | -6.19765 | -6.04461 | -5.9428 | -5.8443 | -5.70916 | -5.66716 | 93881 | 0.999998 |
log_lik[10] | -6.71085 | 0.000707805 | 0.223828 | -7.19011 | -7.1025 | -6.85252 | -6.6971 | -6.55379 | -6.36867 | -6.30991 | 100000 | 0.999989 |
log_lik[11] | -6.21359 | 0.000504989 | 0.159691 | -6.54217 | -6.48589 | -6.31792 | -6.20807 | -6.10355 | -5.96017 | -5.91534 | 100000 | 1.00002 |
log_lik[12] | -6.09369 | 0.00049108 | 0.150984 | -6.40111 | -6.34936 | -6.19326 | -6.08984 | -5.98977 | -5.85288 | -5.80928 | 94527 | 1.00002 |
log_lik[13] | -6.60844 | 0.000670976 | 0.212181 | -7.05939 | -6.97869 | -6.7431 | -6.59501 | -6.461 | -6.28327 | -6.22806 | 100000 | 1.00001 |
log_lik[14] | -6.82174 | 0.000683128 | 0.216024 | -7.28002 | -7.19688 | -6.95918 | -6.80882 | -6.67223 | -6.48941 | -6.43331 | 100000 | 1.00001 |
log_lik[15] | -7.5454 | 0.0011654 | 0.368531 | -8.34713 | -8.19487 | -7.77589 | -7.51674 | -7.28485 | -6.99185 | -6.90637 | 100000 | 0.999982 |
log_lik[16] | -5.99181 | 0.000501445 | 0.14957 | -6.29642 | -6.24462 | -6.09066 | -5.9877 | -5.8892 | -5.7523 | -5.71028 | 88970 | 1.00003 |
log_lik[17] | -8.19007 | 0.00156933 | 0.496266 | -9.27191 | -9.06673 | -8.49926 | -8.15016 | -7.83863 | -7.44239 | -7.33214 | 100000 | 1 |
log_lik[18] | -6.48041 | 0.000548936 | 0.173589 | -6.83737 | -6.77563 | -6.59425 | -6.47343 | -6.35979 | -6.20713 | -6.15985 | 100000 | 1.00001 |
log_lik[19] | -7.98662 | 0.00135638 | 0.428926 | -8.91464 | -8.74592 | -8.25546 | -7.95257 | -7.68182 | -7.34138 | -7.24086 | 100000 | 0.999979 |
log_lik[20] | -6.71568 | 0.000647957 | 0.204902 | -7.14689 | -7.07082 | -6.84698 | -6.70495 | -6.57247 | -6.39928 | -6.34525 | 100000 | 0.999987 |
log_lik[21] | -9.68219 | 0.00249132 | 0.787825 | -11.3858 | -11.0712 | -10.1776 | -9.6205 | -9.12183 | -8.49679 | -8.32186 | 100000 | 1.00001 |
log_lik[22] | -6.85578 | 0.000652113 | 0.206216 | -7.2915 | -7.21196 | -6.98743 | -6.84574 | -6.71229 | -6.53602 | -6.47974 | 100000 | 0.999978 |
log_lik[23] | -6.70807 | 0.000604833 | 0.191265 | -7.10608 | -7.03765 | -6.83147 | -6.69943 | -6.57574 | -6.40914 | -6.35915 | 100000 | 1 |
log_lik[24] | -6.43882 | 0.000529426 | 0.167419 | -6.78252 | -6.724 | -6.54813 | -6.43286 | -6.32297 | -6.17362 | -6.12905 | 100000 | 0.999994 |
log_lik[25] | -7.88994 | 0.0011714 | 0.37043 | -8.69871 | -8.54522 | -8.12159 | -7.86367 | -7.6282 | -7.33072 | -7.24402 | 100000 | 0.999994 |
log_lik[26] | -6.26907 | 0.000533125 | 0.164099 | -6.60602 | -6.54939 | -6.37576 | -6.26344 | -6.15559 | -6.00966 | -5.96361 | 94744 | 1.00001 |
log_lik[27] | -6.90996 | 0.000818656 | 0.258882 | -7.4684 | -7.36431 | -7.0731 | -6.8915 | -6.72761 | -6.51943 | -6.45387 | 100000 | 1.00001 |
log_lik[28] | -6.64837 | 0.000691255 | 0.218594 | -7.11348 | -7.02718 | -6.7864 | -6.63571 | -6.49597 | -6.31147 | -6.25681 | 100000 | 0.999982 |
log_lik[29] | -6.74747 | 0.000709018 | 0.224211 | -7.22664 | -7.13774 | -6.88949 | -6.73451 | -6.59022 | -6.40498 | -6.34879 | 100000 | 0.999986 |
log_lik[30] | -5.92905 | 0.000495427 | 0.148163 | -6.23007 | -6.17892 | -6.02743 | -5.92463 | -5.82661 | -5.69268 | -5.64909 | 89438 | 1.00002 |
log_lik[31] | -5.9427 | 0.000481788 | 0.147672 | -6.24341 | -6.19239 | -6.04019 | -5.93912 | -5.84076 | -5.70632 | -5.66325 | 93948 | 1 |
log_lik[32] | -6.45466 | 0.000551094 | 0.174271 | -6.81692 | -6.7533 | -6.56731 | -6.4479 | -6.33399 | -6.1815 | -6.13349 | 100000 | 1 |
log_lik[33] | -6.39362 | 0.00054583 | 0.172607 | -6.75128 | -6.68942 | -6.50549 | -6.38655 | -6.27466 | -6.12209 | -6.07451 | 100000 | 1.00001 |
log_lik[34] | -6.26907 | 0.000533125 | 0.164099 | -6.60602 | -6.54939 | -6.37576 | -6.26344 | -6.15559 | -6.00966 | -5.96361 | 94744 | 1.00001 |
log_lik[35] | -6.09634 | 0.000515248 | 0.154159 | -6.40996 | -6.35791 | -6.19751 | -6.09192 | -5.99059 | -5.85023 | -5.80673 | 89517 | 1.00002 |
log_lik[36] | -6.42419 | 0.000541343 | 0.171188 | -6.77728 | -6.71547 | -6.53626 | -6.41756 | -6.30555 | -6.15428 | -6.1082 | 100000 | 1.00001 |
log_lik[37] | -5.87456 | 0.000486567 | 0.147509 | -6.1755 | -6.12359 | -5.97152 | -5.87056 | -5.77306 | -5.63864 | -5.59665 | 91908 | 1.00002 |
log_lik[38] | -6.08502 | 0.000474243 | 0.149969 | -6.38969 | -6.33788 | -6.18427 | -6.08081 | -5.98184 | -5.84684 | -5.80406 | 100000 | 1 |
log_lik[39] | -6.58677 | 0.000616932 | 0.195091 | -6.9986 | -6.92179 | -6.71234 | -6.57743 | -6.45114 | -6.2843 | -6.23202 | 100000 | 0.999988 |
log_lik[40] | -6.52813 | 0.000557474 | 0.176289 | -6.89347 | -6.82835 | -6.64266 | -6.52131 | -6.40592 | -6.25027 | -6.20306 | 100000 | 0.999983 |
log_lik[41] | -5.97066 | 0.00048391 | 0.148823 | -6.27342 | -6.2218 | -6.06889 | -5.9667 | -5.86849 | -5.73262 | -5.69075 | 94583 | 1.00001 |
log_lik[42] | -6.7658 | 0.000749145 | 0.2369 | -7.27327 | -7.17892 | -6.91408 | -6.75068 | -6.60047 | -6.40354 | -6.34414 | 100000 | 0.999985 |
log_lik[43] | -8.75266 | 0.00173329 | 0.548115 | -9.93939 | -9.71752 | -9.09476 | -8.70939 | -8.36307 | -7.93373 | -7.80486 | 100000 | 0.999986 |
log_lik[44] | -7.4732 | 0.000864603 | 0.273411 | -8.05226 | -7.94825 | -7.64826 | -7.45718 | -7.28092 | -7.05404 | -6.98571 | 100000 | 0.99998 |
log_lik[45] | -6.86811 | 0.000670866 | 0.212146 | -7.31424 | -7.2334 | -7.00434 | -6.85678 | -6.71946 | -6.53968 | -6.48454 | 100000 | 0.999993 |
log_lik[46] | -8.39545 | 0.00163259 | 0.516269 | -9.51436 | -9.31025 | -8.71891 | -8.35481 | -8.02785 | -7.61944 | -7.50049 | 100000 | 0.999979 |
log_lik[47] | -7.21641 | 0.00090867 | 0.287347 | -7.83597 | -7.72231 | -7.398 | -7.19558 | -7.01405 | -6.78238 | -6.71172 | 100000 | 0.99998 |
log_lik[48] | -6.19707 | 0.000517114 | 0.15826 | -6.52008 | -6.46738 | -6.30066 | -6.19209 | -6.08782 | -5.94499 | -5.90124 | 93663 | 1.00001 |
log_lik[49] | -6.61441 | 0.000670758 | 0.212112 | -7.0674 | -6.98435 | -6.74872 | -6.60168 | -6.46637 | -6.28944 | -6.23353 | 100000 | 1.00001 |
mu_diff | 15.3836 | 0.00555407 | 1.75635 | 11.9336 | 12.5127 | 14.207 | 15.3843 | 16.5573 | 18.2745 | 18.853 | 100000 | 0.999983 |
es | 3.30844 | 0.00160306 | 0.506931 | 2.36452 | 2.50846 | 2.95786 | 3.28843 | 3.63776 | 4.17218 | 4.35433 | 100000 | 0.99998 |
cohenu | 0.90849 | 0.000107709 | 0.0335828 | 0.831913 | 0.847353 | 0.888497 | 0.912552 | 0.932937 | 0.955767 | 0.961685 | 97214 | 0.999988 |
pod | 0.891195 | 0.00010542 | 0.0333368 | 0.816538 | 0.831072 | 0.87091 | 0.89468 | 0.915269 | 0.939324 | 0.945675 | 100000 | 0.999983 |
pbt | 0.667754 | 0.000166733 | 0.0527255 | 0.560294 | 0.578451 | 0.632602 | 0.669254 | 0.704582 | 0.751648 | 0.766491 | 100000 | 0.999974 |
prob_mu_diff_upper_0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 100000 | nan |
prob_mu_diff_upper_c | 0.78579 | 0.00138494 | 0.410275 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 87758 | 0.999979 |
prob_es_upper_c | 0.72119 | 0.00151356 | 0.448416 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 87773 | 0.999991 |
prob_cohenu_upper_c | 0.0856 | 0.000957681 | 0.279774 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 85344 | 0.999981 |
prob_pod_upper_c | 0.01434 | 0.00040641 | 0.118889 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 85576 | 0.999976 |
prob_pbt_upper_cdash | 0.89648 | 0.00108438 | 0.304638 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 78924 | 1.00001 |
問題なく、パラメータ値が収束して計算されました(Rhat<1.1で判定しています)
WAIC = 683.147として計算されます。 なお等分散モデルだとWAIC=716.444なので、実際に不偏分散モデルのほうが感覚としても、情報量としてもデータとモデルがフィットしていることがわかります。
仮説と検証
- A群の平均値がB群より高い確率は?
- prob_mu_diff_upper_0のサンプリング結果より、100%の確率で、A群の方がB群より高いことがわかります。
- AぐんとB群の平均値の差の点推定、95%信頼区間推定
- 15.384(1.756)[11.934, 18.853]となりました。
- 平均値の差の片側95%信頼区間推定の上限、下限はどのあたりか
- mu_diff > 12.513、あるいはmu_diff < 18.275の区間に95%の確率で入る
- 平均値の差が14より大きい確率は?
- prob_mu_diff_upper_cより、78.6%の確率で、14より大きいことがわかります。
- A群から見たB群の効果量の点推定、95%信頼区間推定 両側/片側
- A群から見たB群の効果量が3より大きい確率は?
- prob_es_upper_cより、72.1%の確率で、3より大きいことがわかります。
- A群の分布とB群の分布の重複する割合(非重複度)の点推定、95%信頼区間推定 両側/片側
- 非重複度が0.95より大きい確率
- prob_cohenu_upper_cより、8.6% の確率で、0.95より
- 無作為に選んだ際に、A群の値がB群より大きくなる確率(優越率)はどのぐらいか
- 0.891(0.033)[0.817, 0.946]の確率で高くなることが推定されます。
- 優越率が0.95より大きい確率
- prob_pod_upper_cより 1.4%の確率で0.95より高くなることが推定されます。
- 無作為に選んだそれぞれの測定値の差が10より大きい確率(閾上率)
- 0.668(0.053)[0.560, 0.766]の確率で、大きいことが推定されます。
- 閾上率が0.60(60%)より大きい確率
- prob_pbt_upper_cdashより89.6%の確率で60%より高くなることが推定されます。
以上のように、MCMCによるサンプリングの試行結果を使うことで、様々な情報を抜き出すことが出来ました。 次回は第4章の、対応のあるt検定です。
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その3
その3です。今回は第2章の話をしていきます。
2章 MCMCと正規分布の推測
内容
- 解析的に事後分布を推定するのは難しい!
事後分布に従う乱数を発生させる、という発想転換を行う(マルコフ連鎖モンテカルロ法; MCMC)
- 乱数群の結果を比較して、解析する
- HM法とGS法が具体的なアルゴリズムとしては代表的な存在
- GS法は特定の事前分布であることが必要(やや恣意的)
- HM法は自由な事前分布を使うことが出来る
- 事前分布の設定 : 常識的な知識から決める
- 乱数の発生 : 初期乱数は使用しない(burn in、あるいはwarm upと呼ばれる)
- 値が登ったり下りたり(ドリフトした)形状が観測される場合には、好ましくない乱数生成となる
- 複数chainを発生させて、全部で統一的な傾向が見られてない場合、好ましくない乱数生成が行われていると考えられる
- 収束判定指標(Rcap) : 複数chainの分散, 1.1以下が理想
- 有効標本数(neff) : 乱数の内、いくつが無関係な乱数として扱うことが出来るかを示す指標
- (以上2つの指標の導出が気になった場合、基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門に結構丁寧に書いてあるので、そちらをお薦めします。)
事後分布の解釈をどう行うか?
- 点推定量
- 確信区間(信用区間)
- 予測分布
- 新たなデータが得られうる分布のこと
- 事後予測分布
- 推定された母数の事後分布を用いて、未知データを予測する統計モデルの平均を計算
- 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
- 条件付き予測分布
- 事後分布の特定の点推定量を利用する
- 簡単だが、せっかく推定したモデルを有効活用出来ない
- (これについては、その1で紹介した「はじめての統計データ分析」 豊田秀樹のメモ - StatModeling Memorandumさんの方で問題点が指摘されています)
- (が、豊田先生もこのような指摘を結構受けたのか、本に付録でついてくる追加Q&Aに解説が追記されています(後述))
- ベイズ的推測
- 何らかの仮説を定めて、データから推定すること
- 統計的仮説検定では帰無仮説に基づいたが、そんなことはしなくてよいのがベイズ統計
- 生成量
- 推定した母数を利用して、各種統計量の分布も計算することが出来る
- 母数より生成される値なので、生成量という
- 推定した母数を利用して、各種統計量の分布も計算することが出来る
- 研究仮説が正しい確率
- 真偽を判定する二値変数を利用して、研究仮説が正しい確率を評価することが出来る
- 0 or 1しか使わないので、分散や信頼区間は使えない
条件付き予測分布の取扱について(2017-05-28追記)
ここまでで書いた概要部分でも説明しましたが、議論のされている部分なので補足説明です。 このテキストではデータが得られた状態で、次に得られるであろう未知データについての予測方法を2つ紹介しています。
1つは事後予測分布です。事後予測分布の数式はと記述されます。これを踏まえたやり方です。 数式だけ見ると混乱しがちなのですが、MCMCで計算する時にはそんなに難しくありません。 MCMCの各サンプリングにおいて得られるがを示していますし、このを構築したモデルに代入して得られる値が分布からサンプリングされたうちの1つということになります。 MCMCが積分の代わりをしている、という認識ですね。
もう1つは条件付き予測分布です。これはMCMCの結果として得られるのEAPやMAPといった点推定値をモデルに代入して、予測分布を得るというものです。 これは理解もし易いし、簡単です。分布の中で最もそれっぽいと理解出来た値をモデルに代入して予測をすることができます。 著者の豊田先生としては、
- 理解が簡単
- 速度が速い(生成量や予測量として計算しなくて済む)
- オンライン性(事後予測分布は過去データが必要だが、条件付き予測分布は点推定値だけでよい)
という利点を挙げています。
しかしながらリンク先にもあるように、条件付き予測分布には色々問題点もありますし、何より個人的には「もったいない」という点が目立つと考えています。 せっかくMCMCまでして分布を計算したんですから、それを有効に活用しましょう。 それに「仮説が正しい確率」を計算するには、どっちにしろ事後予想分布が必要になります。
以上の理由で、以後の説明では条件付き予測分布は使用せず、事後予想分布で推定を行っています。
章末問題
データの準備
- 1章のデータをそのまま利用します
import numpy as np from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger # データの準備 x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33]) logger = get_logger()
Stanによる統計モデルの構築
- MCMCをStanで行います
- Stan言語で書かれたファイルを読み込んで、C++にコンパイルすることで、高速にサンプリングを行うことが可能となります。
- 標本データや問題文を見た時に、値の範囲について情報が全く得られなかったので、パラメータ分布には0以上という制限しか利用していません。
import os import pystan import pickle # Stanのモデルを読み込んでコンパイルする stan_file = os.path.join("stan", "g1_mean.stan") stan_file_c = os.path.join("stan", "g1_mean.pkl") model = pystan.StanModel(file=stan_file) with open(stan_file_c, "wb") as f: pickle.dump(model, f)
# Stanのモデル data { int<lower=0> n ; real<lower=0> x[n] ; int<lower=0> c ; real<lower=0, upper=1> percent ; real ces ; } parameters { real<lower=0> mu ; real<lower=0> sigma ; } model { x ~ normal(mu,sigma); } generated quantities{ real xaste ; real sigma2 ; real cv ; real es; real percent_point ; real<lower=0, upper=1> prob_under_c_cdf ; real rate_c ; real<lower=0, upper=1> prob_mu_under_c ; real<lower=0, upper=1> prob_under_c_bool ; real<lower=0, upper=1> prob_es_under_ces ; real<lower=0, upper=1> prob_prob_under_c_bool_over_percent ; xaste = normal_rng(mu,sigma) ; sigma2 = sigma * sigma ; cv = sigma / mu ; es = (mu - c) / sigma ; percent_point = mu + normal_cdf(percent, 0, 1) * sigma ; rate_c = xaste / c ; prob_under_c_cdf = normal_cdf(c, mu, sigma) ; prob_under_c_bool = xaste < c ? 1 : 0 ; prob_mu_under_c = mu < c ? 1 : 0 ; prob_es_under_ces = es < ces ? 1 : 0 ; prob_prob_under_c_bool_over_percent = prob_under_c_bool > percent ? 1 : 0 ; }
統計モデルによる事後分布のサンプリング
- コンパイルしたモデルを使って、正規分布からデータが得られた際に、正規分布の母数がどのように分布するか調べます
- 母数の分布から、各種効果量を推定しています。
- 仮説が得られる確率については、各サンプリング値において判定を行い、その結果のbool値の頻度で確認します。
- 付録でついてきたスクリプトにはここらへんの記述がありませんので、何か問題が有った場合には指摘頂けると助かります。
import pandas as pd import pickle import pystan import matplotlib import matplotlib.pyplot as plt from IPython.core.display import display %matplotlib inline matplotlib.rcParams['figure.figsize'] = (10, 40) # Stanで使うデータの用意 stan_data = {"n": x.size, "x": x, "c": 41, "percent": 0.70, "ces": -0.05} # 興味のあるパラメータの設定 par = ["mu", "sigma", "xaste", "sigma2", "cv", "es", "percent_point", "rate_c", "prob_under_c_cdf", "prob_under_c_bool", "prob_mu_under_c", "prob_es_under_ces", "prob_prob_under_c_bool_over_percent"] # 興味のある母数の%点を書く prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # コンパイルしたモデルの読み込み stan_file_c = os.path.join("stan", "g1_mean.pkl") with open(stan_file_c, "rb") as f: model = pickle.load(f) # MCMCでサンプリング fit = model.sampling(data=stan_data, pars=par, iter=21000, chains=5, warmup=1000, seed=1234, algorithm="HMC") # 事後分布の表を取得、可視化 summary = fit.summary(pars=par, probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) display(summary_df) # 事後分布の可視化 fit.plot(pars=par) plt.show()
- 結果(表と分布)
mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
mu | 40.660 | 0.010 | 0.649 | 39.383 | 39.603 | 40.224 | 40.659 | 41.091 | 41.735 | 41.947 | 4120.000 | 1.004 |
sigma | 6.523 | 0.014 | 0.468 | 5.660 | 5.797 | 6.200 | 6.499 | 6.828 | 7.321 | 7.488 | 1154.000 | 1.008 |
xaste | 40.693 | 0.021 | 6.560 | 27.784 | 29.927 | 36.314 | 40.679 | 45.108 | 51.424 | 53.495 | 96434.000 | 1.000 |
sigma2 | 42.772 | 0.182 | 6.178 | 32.033 | 33.604 | 38.445 | 42.237 | 46.623 | 53.602 | 56.063 | 1153.000 | 1.008 |
cv | 0.160 | 0.000 | 0.012 | 0.139 | 0.142 | 0.152 | 0.160 | 0.168 | 0.181 | 0.185 | 1372.000 | 1.006 |
es | -0.052 | 0.002 | 0.099 | -0.248 | -0.215 | -0.120 | -0.053 | 0.014 | 0.112 | 0.144 | 4065.000 | 1.004 |
percent_point | 45.605 | 0.019 | 0.744 | 44.204 | 44.418 | 45.099 | 45.583 | 46.092 | 46.862 | 47.132 | 1509.000 | 1.009 |
rate_c | 0.993 | 0.001 | 0.160 | 0.678 | 0.730 | 0.886 | 0.992 | 1.100 | 1.254 | 1.305 | 96434.000 | 1.000 |
prob_under_c_cdf | 0.521 | 0.001 | 0.039 | 0.443 | 0.455 | 0.494 | 0.521 | 0.548 | 0.585 | 0.598 | 4055.000 | 1.004 |
prob_under_c_bool | 0.520 | 0.002 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 97815.000 | 1.000 |
prob_mu_under_c | 0.703 | 0.007 | 0.457 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 4861.000 | 1.001 |
prob_es_under_ces | 0.510 | 0.007 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 5481.000 | 1.002 |
prob_prob_under_c_bool_over_percent | 0.520 | 0.002 | 0.500 | 0.000 | 0.000 | 0.000 | 1.000 | 1.000 | 1.000 | 1.000 | 97815.000 | 1.000 |
Rhat値が全パラメータで1.1以下になっていることや、Traceplotの様子からMCMCが収束していることが確認できると思います。 というわけで推定値は一応問題がありません。実際に推定値を用いて推論を行うことが出来ます。
仮説と検証
平均値
- このデータにおける平均値はいくつか
- 事後分布のEAPが40.660(0.650)であるため、40.660と推定される
- このデータにおいて、平均値が95%の確率で分布する区間はどこか
- このデータにおいて平均値が片側で95%の信頼区間を持つ区間はどこか
- 事後分布の平均値の95%信頼区間より、41.735以下になる確率が95%である
- もしくは39.603以上になる確率が95%である
標準偏差
予測分布
- 新しい人からデータをとったときに、95%の確率で何点を取るでしょうか
- 27.769から53.520点の間のいずれかの点数をとる
生成量
- 分散の点推定と95%信頼区間推定
- 42.772(6.178)[32.033, 56.063]
- 基準値を41とした際の、変動係数の点推定と95%信頼区間推定
- 0.160(0.012)[0.139, 0.185]
- 基準値を41とした際の、効果量の点推定
- -0.524(0.099)
- 基準値を41とした際の、効果量の95%区間推定(両側/片側)
- 両側[-0.248, 0.144]、片側[-0.215, ] or [, 0.112]
- 70%点の点推定と区間推定
- 45.605(0.744)[44.204, 47.132]
- 41以下の値が観測される確率と95%信頼区間
- 分布関数を利用 : 0.521(0.039)[0.443, 0.598]
- 確率的命題を利用 : 0520
- 41と測定値の比の、点推定値と95%信頼区間推定値
- 0.993
研究仮説が正しい確率
- 平均値が41より小さい確率
- 0.703
- 基準値を41とした際の効果量が、-0.05より小さい確率
- 0.510
- 41未満の測定値が観測される確率が、70%より大きい確率
- 0.512
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その2
前置き
環境構築
まずPythonを実行する環境を作ります。 一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。 以下のものはOSXを想定しています。
自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。 導入は各リポジトリのページに書いてあるとおりです。
pythonは3系を使います。今回利用するライブラリはそんなに変なものではないので、2系である必要はありません。 解析はスクリプトでやるよりも、jupyterを使う方が僕は好きです。Stan言語をいじれない問題点はありますが、パラメータを変えるのが楽です。
pyenv install 3.6.1 pyenv virtualenv 3.6.1 3.6.1_postp mkdir -p ~/Desktop/postp # 作業ディレクトリは適当に変えてください。 cd ~/Desktop postp pyenv local 3.6.1_postp pip install numpy scipy pandas matplotlib cython seaborn jupyter pip install pystan
1章 データの整理とベイズの定理
内容
一応箇条書きで書いていますが、メモです。 ちゃんと勉強したいときは、本文を参考にしてくださいね。
- 分布の紹介
- 正規分布
- 一様分布
- データの解釈
章末問題
- pythonにて解析をしていきます。
データの準備
僕の好みですが、printではなくloggingモジュールを利用して途中経過を書いていきます。計算時間なども把握できるのが便利です。
import numpy as np from logging import getLogger, Formatter, StreamHandler, DEBUG # printではなくloggerを使う def get_logger(): logger = getLogger(__name__) logger.setLevel(DEBUG) log_fmt = '%(asctime)s : %(name)s : %(levelname)s : %(message)s' formatter = Formatter(log_fmt) stream_handler = StreamHandler() stream_handler.setLevel(DEBUG) stream_handler.setFormatter(formatter) logger.addHandler(stream_handler) return logger # データの準備 x = np.array([36,38,51,40,41,52,43,31,35,37,49,43,43,41,36,53,43,26,45,37,33,38,33,35,36,28,46,41,32,49,43,38,46,46,46,45,44,40,38,37,35,39,31,55,48,32,37,37,45,39,42,40,40,50,38,51,29,44,41,42,43,36,38,33,32,42,43,40,46,54,37,24,47,35,35,47,38,31,41,39,40,43,37,45,38,42,48,43,38,48,47,44,42,36,50,36,55,51,38,33]) logger = get_logger()
度数分布表
いきなりですが、pandasでは階級幅を決めるメソッドが無いので超めんどくさい… 手動でデータを入れるbinsを作成した後、pd.cutを使って仕分けます。
import pandas as pd df = pd.DataFrame(x, columns=["weight"]) width = 5 minim_range = int((df.min() / width).apply(np.floor) * width) maxim_range = int((df.max() / width).apply(np.ceil) * width) bins = range(minim_range, maxim_range + 1, width) c = pd.cut(df["weight"], bins) d_df = df.groupby(c).count() d_df.columns = ["count"] d_df.columns.name = "range" d_df.index.name = None d_df
ヒストグラムの作成
seabornを使って可視化します。matplotlibではないのは、やはり好みです pandas.plotと違ってseaborn / matplotlibのhistメソッドはbinsに配列を受け取ってくれます。
import seaborn as sns import pandas as pd %matplotlib inline df = pd.DataFrame(x, columns=["weight"]) width = 5 minim_range = int((df.min() / width).apply(np.floor) * width) maxim_range = int((df.max() / width).apply(np.ceil) * width) bins = range(minim_range, maxim_range + 1, width) sns.distplot(df["weight"], kde=False, bins=bins)
統計量
最頻値だけはメソッドが無いので、countをとって最大の値を持つものを返しています。 pandasのvar, stdメソッドはデフォルトパラメータだと母分散、母標準偏差を返す事に注意してください。 numpyは標本分散、標本標準偏差を返す。謎の挙動。
import pandas as pd df = pd.DataFrame(x, columns=["weight"]) # 標本平均 df_mean = df["weight"].mean() logger.info("平均値は{0}".format(df_mean)) # 標本分散 df_var = df["weight"].var(ddof=0) logger.info("標本分散は{0}".format(df_var)) # 標本標準偏差 df_std = df["weight"].std(ddof=1) logger.info("標本標準偏差は{0}".format(df_std)) # データのソート df["weight"].sort_values() # 最大値 df_max = df["weight"].max() logger.info("最大値は{0}".format(df_max)) # 最小値 df_min = df["weight"].min() logger.info("最小値は{0}".format(df_min)) # 中央値 df_medium = df["weight"].median() logger.info("中央値は{0}".format(df_min)) # 最頻値 df_most_freq = df["weight"].value_counts().iloc[0] logger.info("最頻値は{0}".format(df_most_freq))
正規分部の扱い
scipy.statsのnormオブジェクトを利用します。 片側だけから特定の確率になる値を求めるメソッドが無いですが、intervalメソッドで両側の場合の値を求めて、2で割って考えれば良いです。 これは正規分布の対称性を利用したものですね。
import pandas as pd from scipy.stats import norm df = pd.DataFrame(x, columns=["weight"]) # パラメータを設定 df_mean = df["weight"].mean() df_std = df["weight"].std(ddof=1) logger.info("正規分布のパラメータ: 平均は{0}、標準偏差は{1}".format(df_mean, df_std)) # 分布の平均が40.64なので、40付近の値が観測されやすい # 45以上の値が観測される確率 # 1から45までの値が観測される確率を引く(1-cdf) # scipyでは、sfを使えば良い from scipy.stats import norm prob_upper_45 = norm.sf([45], loc=df_mean, scale=df_std)[0] logger.info("45以上の値が観測される確率は{0}".format(prob_upper_45)) # 35以上45未満の値が観測される確率 # 45以下が観測される確率から、35以下が観測される確率を引く prob_upper_35 = norm.cdf([35], loc=df_mean, scale=df_std)[0] prob_upper_45 = norm.cdf([45], loc=df_mean, scale=df_std)[0] prob_between_35_45 = prob_upper_45 - prob_upper_35 logger.info("35以上45未満の値が観測される確率は{0}".format(prob_between_35_45)) # 95%信頼区間 # 正規分布において、平均 ± 1.96 * 標準偏差の範囲が該当することを利用する #min_edge = df_mean - df_std * 1.96 #max_edge = df_mean + df_std * 1.96 #prob_min = norm.cdf([min_edge], loc=df_mean, scale=df_std)[0] #prob_max = norm.cdf([max_edge], loc=df_mean, scale=df_std)[0] # prob_sig_range = prob_max - prob_min # が、scipyだとintervalを利用すれば良い min_edge, max_edge = norm.interval(0.95, loc=df_mean, scale=df_std) logger.info("95%信頼区間は{0}から{1}".format(min_edge, max_edge)) # p(x>a) = 0.05であるような点a # F(平均 + 1.64 * 標準偏差)が大体0.95になる #prob_less_095 = norm.sf([df_mean + df_std * 1.64], loc=df_mean, scale=df_std)[0] # 何も情報が無いとしたら、正規分布性を元に、90%信頼区間を出して、上端を利用すれば良い _, less_095 = norm.interval(0.9, loc=df_mean, scale=df_std) prob_less_095 = norm.cdf([less_095], loc=df_mean, scale=df_std)[0] logger.info("a = {0} であれば、p(x>a)となる確率は{1}である".format(less_095, prob_less_095)) # 3つの4分位点 # 50%信頼区間と中央値を出せば良い quad_first, quad_third = norm.interval(0.5, loc=df_mean, scale=df_std) quad_second = df_medium logger.info("第1四分位点は{0}, 第2四分位点は{1}, 第3四分位点は{2}".format(quad_first, quad_second, quad_third))
Rで出来ることがPythonだとめんどくさかったりしますね。 特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります!
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学― その1
このブログの開設目的の、
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―
- 作者: 豊田秀樹
- 出版社/メーカー: 朝倉書店
- 発売日: 2016/06/02
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (11件) を見る
についての記事インデックスです。
目次
- その1 ←このページ
- その2 環境構築と第1章 データ整理とベイズ概要
- その3 第2章 : MCMCによる正規分布のパラメータ推定
- その4 第3章 : 独立した2群の差の推定
- その4.1 第3章おまけ : 変分ベイズによる推定
- その5 第4章 : 対応ある2群の差の推定
- その6 第5章 : 要因毎の効果推定
- その7 第6章 : カウントデータと比率データの推定
本の紹介
この本は、昨今はびこる統計データ分析についてベイジアンの知見から解説した本です。 他の方の感想では、 statmodeling.hatenablog.com
のようなものがあります。個人的に一冊まるまるよんだ感想としては、次のような印象を持ちました。
良い点
MCMCについての細かい知識は要らない
ベイズの定理やMCMCを利用した多くの解析の本では、数式がガッツリ出てきて頭を抱える方も多いと思います。 しかし、この本ではStan言語を利用して解析を行っていきます。そのため細かい数式の変形や、推論の部分についてはStanがなんとかしてくれますので、本にも書いてありませんし、やる必要もありません。
Research Questionベースの解析
実際にデータ分析をする際には、様々なデータ提供者からの要望に答える形で解析をすることが多いと思います。 この本では、そのような問いをResearch Questionと呼び、各章で該当するResearch Questionを最初に挙げて、解決する形で進めていきます。 結果として、どのような形で今やっている解析が活きるものなのか、具体的にイメージしながら進める事ができます。
ベイズベースでのデータ解析への「入門」
豊田先生が後書きにも記していますが、普通データ解析や統計学を勉強するときには、頻度主義に基づく統計的仮説検定を勉強します(僕もそうしました)。 一方でベイズの定理を基づくベイジアン統計学は、最近計算機処理が高速化したことやデータ量が多くなっていることが原因で、頻繁に使われています。 しかしベイジアン統計学は、頻度主義とはぜんぜん違うため、入門が難しいと僕も思います。その中で本書は、統計的仮説検定の手法でやるような内容をベイジアン統計学で行った場合にどうなるかを解説しています。 その点で、最初に取り組む際に何が利点となるのか分かりやすいと思います。
答えが用意されている
最近は多いですが、解析に使っているスクリプトがR + Stanで用意されています。
朝倉書店|はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―
更に、各章末問題は巻末に答えがかなり丁寧に書いてあります。 大学の教科書もそうですが、答えがない場合自分が間違っていても気が付かないことがあります。答えがあるのは非常に助かりますね。
悪い点
入門書ではない
幾つもの指摘がありましたが、入門書ではありません。 たとえば、
といった点が目立ちました。1つ目、2つ目については久保先生の通称緑本、ベイズの定理について入門レベルから丁寧に勉強したい場合には超入門が良かったと思います。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (29件) を見る
図解・ベイズ統計「超」入門 あいまいなデータから未来を予測する技術 (サイエンス・アイ新書)
- 作者: 涌井貞美
- 出版社/メーカー: SBクリエイティブ
- 発売日: 2013/12/17
- メディア: 新書
- この商品を含むブログ (15件) を見る
なにをするのか
以降からは、各章の簡単な解説をした上で、Pythonで章末問題を解いていきたいと思います。 Stan言語については、PythonのインターフェイスであるPyStanを利用します。 といっても、Stan言語で書かれている部分についてはRと共通して使うことが可能なので、PyStanを使うためにどのような小手先のテクニックがあるかを書いていくつもりです。
ブログ始めました。
自己紹介
ブログ始めました。
普段は次世代シーケンサー(NGS)を用いた遺伝統計学や、ゲノム配列から得られる知見の解析、それらのツール・DB開発などを行っています。
ブログの目的
自身で勉強したことを内部のWikiに書いていたのですが、もっと外部に公開してマサカリを投げてもらったり、何かしらの役に立つことを期待してブログを開設します。 基本的には、
Javascriptを用いたデータ可視化(D3.js)、アプリケーション開発
次世代シーケンサーから得られるデータの解析について
を主に書いていきたいと思います。特にPythonについては日本語知識があまり転がっていないので、Rや海外サイトから輸入することが多いかもしれません。 3番目のNGS解析については、需要がなさそうなのであまり書かない気がします笑
それでは宜しくお願いします。Amazon.co.jpアソシエイトです。