その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

f:id:ajhjhaf:20170531011835p:plain

2017-05-27

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

ベイジアンモデリング

対応のあるt検定を行う前に、3章の内容を等分散モデルで書くとどうなるか、をちゃんと検証します。

番外編 3章の等分散モデルと変分ベイズによる推定

等分散モデルによる差の推定

といっても、Stanのモデルをちょっと変えるだけなのでした。

まずコンパイル

import os
import pystan
import pickle

# Stanのモデルを読み込んでコンパイルする
# 等分散モデル
stan_file = os.path.join("stan", "g2_ind_equ.stan")
stan_file_c = os.path.join("stan", "g2_ind_equ.pkl")
model = pystan.StanModel(file=stan_file)
with open(stan_file_c, "wb") as f:
    pickle.dump(model, f)

Stanのモデルはこんな感じです。

data {
    int<lower=0> n_a ;
    int<lower=0> n_b ;
    real<lower=0> a[n_a] ;
    real<lower=0> b[n_b] ;
    real<lower=0> c_mu_diff ;
    real<lower=0> c_es ;
    real<lower=0> c_cohenu ;
    real<lower=0> c_pod ;
    real<lower=0> c_pbt ;
    real<lower=0> cdash_pbt ;
}

parameters {
    real<lower=0> mu_a ;
    real<lower=0> sigma ;
    real<lower=0> mu_b ;
}

model {
    a ~ normal(mu_a, sigma) ;
    b ~ normal(mu_b, sigma) ;
}

generated quantities {
    vector[n_a] log_lik ;
    real mu_diff ;
    real es ;
    real cohenu ;
    real pod ;
    real pbt ;
    int<lower=0, upper=1> prob_mu_diff_upper_0 ;
    int<lower=0, upper=1> prob_mu_diff_upper_c ;
    int<lower=0, upper=1> prob_es_upper_c ;
    int<lower=0, upper=1> prob_cohenu_upper_c ;
    int<lower=0, upper=1> prob_pod_upper_c ;
    int<lower=0, upper=1> prob_pbt_upper_cdash ;

    for(i in 1:n_a){
        log_lik[i] = normal_lpdf(a[i] | mu_a, sigma) + normal_lpdf(b[i] | mu_b, sigma) ;
    }
    mu_diff = mu_a - mu_b ;
    es = mu_diff / sigma ;
    cohenu = normal_cdf(mu_a, mu_b, sigma) ;
    pod = normal_cdf(es / sqrt(2), 0, 1) ;
    pbt = normal_cdf((mu_diff - c_pbt) / ( sqrt(2) * sigma ), 0, 1) ;
    prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ;
    prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ;
    prob_es_upper_c = es > c_es ? 1 : 0 ;
    prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ;
    prob_pod_upper_c = pod > c_pod ? 1 : 0 ;
    prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ;
}

後はサンプリングを同様に行うだけ。

import pandas as pd
import pickle
import pystan
import matplotlib
import os
import matplotlib.pyplot as plt
from IPython.core.display import display
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (10, 50)
# 等分散モデル

# Stanで使うデータの用意
stan_data = {"n_a": a.size,
             "n_b": b.size,
             "a": a,
             "b": b,
             "c_mu_diff": 14, #標本平均の差を使ってみる
             "c_es": 3.0, # 標本効果量で、a群からみたもの
             "c_cohenu": 0.95, # a群から見た非重複度
             "c_pod": 0.95, # 優越率
             "c_pbt": 10, # 閾上率の基準値
             "cdash_pbt": 0.60} #閾上率
# 興味のあるパラメータの設定
pars = ["mu_a",
       "sigma",
       "mu_b",
       "log_lik",
       "mu_diff",
       "es",
       "cohenu",
       "pod",
       "pbt",
       "prob_mu_diff_upper_0",
       "prob_mu_diff_upper_c",
       "prob_es_upper_c",
       "prob_cohenu_upper_c",
       "prob_pod_upper_c",
       "prob_pbt_upper_cdash"]
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

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

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


# 事後分布の表を取得
summary = fit.summary(pars=pars, probs=prob)
summary_df = pd.DataFrame(summary["summary"],
                          index=summary["summary_rownames"],
                          columns=summary["summary_colnames"])
display(summary_df)

# 事後分布の可視化
for par in summary_df.index[summary_df["sd"] == 0]:
    pars.remove(par)
fit.traceplot(pars)
plt.show()

# WAICの計算
log_lik = fit.extract("log_lik")["log_lik"]
waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))
logger.info("WAICの値は{0}です".format(waic))

これを実行すると、s.d.についてはEAPが8.63と推定されます。a群とb群のちょうど中間ぐらいですね。両方の分布に同じ分散を使っているので、推定される値としては納得できるものです。ただ尤度については716と不等分散モデルより高くなってしまいました。よってデータにはフィッティングしていないと捉えます。

変分ベイズによる推定

これまでは事後分布の推定を、NUTSアルゴリズムによるサンプリングで行ってきました。ところでStanには、MCMCサンプリング以外にも変分ベイズによる分布の推定も一応実装されています。ついでにこれも試してみましょう。モデルのコンパイル部分までは、上述のものと共通です。

# サンプリングをADVIでやってみる

import pandas as pd
import pickle
import pystan
import matplotlib
import os
import matplotlib.pyplot as plt
from IPython.core.display import display
from collections import OrderedDict
from pystan.external.pymc import plots

%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (10, 50)
# 等分散モデル

# Stanで使うデータの用意
stan_data = {"n_a": a.size,
             "n_b": b.size,
             "a": a,
             "b": b,
             "c_mu_diff": 14, #標本平均の差を使ってみる
             "c_es": 3.0, # 標本効果量で、a群からみたもの
             "c_cohenu": 0.95, # a群から見た非重複度
             "c_pod": 0.95, # 優越率
             "c_pbt": 10, # 閾上率の基準値
             "cdash_pbt": 0.60} #閾上率
# 興味のあるパラメータの設定
pars = ["mu_a",
       "sigma",
       "mu_b",
       "log_lik",
       "mu_diff",
       "es",
       "cohenu",
       "pod",
       "pbt",
       "prob_mu_diff_upper_0",
       "prob_mu_diff_upper_c",
       "prob_es_upper_c",
       "prob_cohenu_upper_c",
       "prob_pod_upper_c",
       "prob_pbt_upper_cdash"]
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

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

# MCMCでADVIサンプリング
fit = model.vb(data=stan_data,
               pars=pars,
               iter=100000,
               seed=1234)


# 事後分布の表を取得
# ADVIはAPIがまだ完成されていないので、summaryの表を作る方式が違う
# chainを利用したサンプリングもしていないため、Rhatも計算できない
vb_sample = pd.read_csv(fit["args"]["sample_file"].decode("utf-8"), comment="#")
vb_sample = vb_sample.drop("lp__", 1)
summary_df = vb_sample.describe(percentiles=prob).T
display(summary_df)

# カーネル密度推定出来ないパラメータの削除
for par in summary_df.index[summary_df["std"] == 0]:
    pars.remove(par)
## 事後分布の可視化
od = OrderedDict()
for i, par in enumerate(fit["sampler_param_names"]):
    par_s = par.split(".")
    if len(par_s) == 1:
        od[par] = np.array(fit["sampler_params"][i])
    else:
        par = par_s[0]
        if par in od.keys():
            od[par] = np.vstack([od[par], np.array(fit["sampler_params"][i])])
        else:
            od[par] = np.array(fit["sampler_params"][i])
plots.traceplot(od, pars)
plt.show()

# WAICの計算
log_lik = od["log_lik"] 
waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))
logger.info("WAICの値は{0}です".format(waic))

注意すべきところは、PyStanの変分ベイズのAPIが完成されていないためsamplingとやり方がぜんぜん違う点です。サンプリング結果の集計にはpandasのdescribeメソッドを使っています。

結果は次のとおりです。

	count	mean	std	min	2.5%	5%	25%	50%	75%	95%	97.5%	max
mu_a	1001	56.4035	1.14076	52.6137	54.1746	54.6317	55.6578	56.3594	57.2304	58.3143	58.6232	59.8892
sigma	1001	7.8148	0.593153	5.79781	6.71698	6.90312	7.43062	7.78376	8.18896	8.82866	9.04493	10.1979
mu_b	1001	40.3831	1.19016	36.5435	38.1027	38.3704	39.5665	40.4121	41.1783	42.2879	42.7135	44.4097
log_lik.1	1001	-6.42276	0.16823	-7.06181	-6.76429	-6.6899	-6.53466	-6.42269	-6.31037	-6.15472	-6.10705	-5.86712
log_lik.2	1001	-6.07869	0.150718	-6.67179	-6.3649	-6.31169	-6.17576	-6.08422	-5.98092	-5.81884	-5.78335	-5.58263
log_lik.3	1001	-8.86784	0.466876	-10.957	-9.91909	-9.69209	-9.1477	-8.81625	-8.52498	-8.19867	-8.08506	-7.54151
log_lik.4	1001	-6.23111	0.153535	-6.66425	-6.54282	-6.48627	-6.33461	-6.2304	-6.12679	-5.96821	-5.93276	-5.71942
log_lik.5	1001	-6.0363	0.152214	-6.58138	-6.32926	-6.28372	-6.13988	-6.03141	-5.93468	-5.77738	-5.73873	-5.49402
log_lik.6	1001	-9.86775	0.610265	-12.5208	-11.2156	-11.0038	-10.2397	-9.81467	-9.41468	-9.01542	-8.86263	-8.08728
log_lik.7	1001	-6.96474	0.214733	-7.70444	-7.41328	-7.32943	-7.10373	-6.95857	-6.81509	-6.63545	-6.57965	-6.26646
log_lik.8	1001	-7.96041	0.353478	-10.3051	-8.72111	-8.56551	-8.16825	-7.93375	-7.71524	-7.44256	-7.35968	-7.04291
log_lik.9	1001	-6.55416	0.179014	-7.26925	-6.9163	-6.838	-6.66726	-6.55088	-6.42996	-6.27155	-6.21939	-6.00683
log_lik.10	1001	-6.20992	0.155187	-6.83945	-6.50715	-6.45275	-6.31807	-6.21339	-6.10742	-5.94728	-5.89824	-5.662
log_lik.11	1001	-7.74777	0.313827	-9.05947	-8.43614	-8.28447	-7.94685	-7.71843	-7.52812	-7.30075	-7.21452	-6.86988
log_lik.12	1001	-6.23652	0.159109	-6.74122	-6.5479	-6.49699	-6.34863	-6.22906	-6.12679	-5.9821	-5.94032	-5.73857
log_lik.13	1001	-6.24006	0.160162	-6.76222	-6.54558	-6.50324	-6.35045	-6.23568	-6.134	-5.98247	-5.93949	-5.70956
log_lik.14	1001	-6.3109	0.159837	-6.75209	-6.61471	-6.57248	-6.42249	-6.30973	-6.19668	-6.06357	-6.02272	-5.77449
log_lik.15	1001	-6.76248	0.194958	-7.7655	-7.14874	-7.08354	-6.89143	-6.75921	-6.62693	-6.45886	-6.40091	-6.21029
log_lik.16	1001	-9.50787	0.558317	-12.0117	-10.7875	-10.4896	-9.84995	-9.46429	-9.09205	-8.71395	-8.58059	-7.86643
log_lik.17	1001	-6.23197	0.159713	-6.77176	-6.55355	-6.51252	-6.33496	-6.22916	-6.12582	-5.97284	-5.91782	-5.47086
log_lik.18	1001	-10.937	0.792785	-16.4895	-12.6988	-12.3467	-11.3937	-10.8907	-10.4062	-9.72509	-9.5612	-8.87464
log_lik.19	1001	-6.55719	0.178717	-7.10137	-6.92829	-6.86605	-6.66928	-6.54786	-6.43631	-6.2802	-6.22669	-5.77744
log_lik.20	1001	-6.92278	0.21248	-7.83707	-7.3648	-7.28035	-7.0519	-6.91023	-6.78217	-6.60749	-6.55017	-6.33407
log_lik.21	1001	-7.55595	0.295133	-9.58726	-8.17106	-8.06674	-7.72853	-7.5413	-7.3452	-7.10152	-7.04768	-6.85252
log_lik.22	1001	-7.55767	0.281643	-9.09015	-8.19547	-8.05543	-7.72132	-7.53802	-7.36901	-7.14595	-7.075	-6.81468
log_lik.23	1001	-7.35302	0.270502	-8.89604	-7.91456	-7.82611	-7.51521	-7.32474	-7.17505	-6.95945	-6.88067	-6.60574
log_lik.24	1001	-6.94368	0.215861	-8.22401	-7.36714	-7.31273	-7.08547	-6.93466	-6.79251	-6.616	-6.54842	-6.35364
log_lik.25	1001	-6.71402	0.191982	-7.69412	-7.09788	-7.03066	-6.84066	-6.71253	-6.58261	-6.40948	-6.34514	-6.1755
log_lik.26	1001	-9.53694	0.584227	-13.4754	-10.843	-10.567	-9.86342	-9.49629	-9.13708	-8.68395	-8.54657	-8.1407
log_lik.27	1001	-6.78402	0.201555	-7.52225	-7.20148	-7.13041	-6.91297	-6.78015	-6.64539	-6.46899	-6.42453	-6.18105
log_lik.28	1001	-6.44432	0.167778	-6.96219	-6.78804	-6.72768	-6.55666	-6.44535	-6.32604	-6.18597	-6.13832	-5.85676
log_lik.29	1001	-7.69239	0.315538	-9.81127	-8.36852	-8.23455	-7.87711	-7.67732	-7.4738	-7.2188	-7.15845	-6.87702
log_lik.30	1001	-7.7927	0.319736	-9.13938	-8.50667	-8.34904	-7.99599	-7.76069	-7.5676	-7.33473	-7.24417	-6.90979
log_lik.31	1001	-6.16013	0.157118	-6.71291	-6.45853	-6.41783	-6.26426	-6.15737	-6.05708	-5.90701	-5.86343	-5.49712
log_lik.32	1001	-6.13017	0.151897	-6.74355	-6.41578	-6.37076	-6.23383	-6.13207	-6.0305	-5.87105	-5.83521	-5.59475
log_lik.33	1001	-7.00073	0.223585	-7.81207	-7.48043	-7.40014	-7.14758	-6.99534	-6.84432	-6.65199	-6.60361	-6.3401
log_lik.34	1001	-6.92585	0.215031	-7.65724	-7.38092	-7.30914	-7.06934	-6.91029	-6.78167	-6.59183	-6.54201	-6.25093
log_lik.35	1001	-6.78402	0.201555	-7.52225	-7.20148	-7.13041	-6.91297	-6.78015	-6.64539	-6.46899	-6.42453	-6.18105
log_lik.36	1001	-6.45221	0.172966	-7.0252	-6.79691	-6.752	-6.56508	-6.448	-6.33613	-6.17511	-6.12069	-5.69448
log_lik.37	1001	-6.43874	0.169585	-6.94675	-6.7849	-6.73926	-6.54679	-6.43412	-6.32283	-6.16711	-6.11841	-5.65075
log_lik.38	1001	-5.98581	0.149329	-6.52304	-6.27182	-6.22623	-6.0843	-5.98738	-5.88368	-5.74347	-5.69836	-5.48248
log_lik.39	1001	-6.22339	0.155309	-6.81162	-6.52915	-6.47324	-6.33043	-6.22686	-6.12094	-5.9492	-5.9203	-5.72975
log_lik.40	1001	-6.38489	0.163779	-6.94275	-6.72587	-6.64711	-6.49179	-6.38321	-6.27196	-6.11557	-6.0736	-5.90043
log_lik.41	1001	-6.66083	0.18857	-7.36905	-7.064	-6.97551	-6.77774	-6.65256	-6.5262	-6.36894	-6.30808	-6.14803
log_lik.42	1001	-6.03561	0.149328	-6.59567	-6.32253	-6.27471	-6.13574	-6.03972	-5.93328	-5.79028	-5.74892	-5.56119
log_lik.43	1001	-7.9404	0.350003	-10.3391	-8.68739	-8.53891	-8.1428	-7.92091	-7.69262	-7.41976	-7.33845	-7.05618
log_lik.44	1001	-11.5909	0.862142	-15.4103	-13.5076	-13.1583	-12.1299	-11.5137	-10.9734	-10.3584	-10.1675	-9.0104
log_lik.45	1001	-7.77097	0.315238	-9.01645	-8.47879	-8.3403	-7.97236	-7.74171	-7.5439	-7.3051	-7.22634	-6.98188
log_lik.46	1001	-7.6209	0.303068	-9.7381	-8.25764	-8.14321	-7.80149	-7.6036	-7.40669	-7.16314	-7.0904	-6.8879
log_lik.47	1001	-6.97629	0.217257	-7.94213	-7.42008	-7.33636	-7.10676	-6.9649	-6.82924	-6.64425	-6.58657	-6.38334
log_lik.48	1001	-6.66957	0.186988	-7.39538	-7.06688	-6.98032	-6.78594	-6.65865	-6.54712	-6.38723	-6.3303	-6.12713
log_lik.49	1001	-6.63225	0.187684	-7.31207	-7.0097	-6.94733	-6.75122	-6.63039	-6.50275	-6.33472	-6.2794	-6.0357
log_lik.50	1001	-6.32385	0.159621	-6.76867	-6.63689	-6.58101	-6.43788	-6.32294	-6.2061	-6.07091	-6.03085	-5.80358
mu_diff	1001	16.0204	1.66137	10.9645	12.6538	13.2746	14.9254	16.0157	17.0902	18.7446	19.2714	20.8654
es	1001	2.06108	0.259806	1.27133	1.58195	1.65403	1.87671	2.05535	2.23764	2.50222	2.59257	2.90178
cohenu	1001	0.977041	0.0139309	0.898194	0.94317	0.95094	0.969721	0.980078	0.987378	0.993829	0.995237	0.998145
pod	1001	0.924156	0.0256932	0.815664	0.868346	0.878915	0.907751	0.926937	0.943204	0.961581	0.966615	0.979909
pbt	1001	0.705715	0.0529771	0.538187	0.593225	0.61623	0.671596	0.708586	0.741193	0.786155	0.803224	0.849742
prob_mu_diff_upper_0	1001	1	0	1	1	1	1	1	1	1	1	1
prob_mu_diff_upper_c	1001	0.879121	0.32615	0	0	0	1	1	1	1	1	1
prob_es_upper_c	1001	0	0	0	0	0	0	0	0	0	0	0
prob_cohenu_upper_c	1001	0.951049	0.215874	0	0	1	1	1	1	1	1	1
prob_pod_upper_c	1001	0.152847	0.36002	0	0	0	0	0	0	1	1	1
prob_pbt_upper_cdash	1001	0.97003	0.17059	0	0	1	1	1	1	1	1	1

f:id:ajhjhaf:20170525183254p:plain

2017-05-25

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

ベイジアンモデリング

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

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

内容

2群に分けた群比較をする際には、ランダム化による交絡因子の排除が重要！(ランダマイゼーション)
2群の、どんな差を見たいかで仮説は変わる
- 群1の平均が群2の平均を上回る確率
- 群1, 2の平均値の差のt年推定、区間推定
- 群1, 2の平均値の差が基準点cより大きい確率
- 効果量が基準点cより大きい確率
- 非重複度、優越率、閾上率…etc
2群でモデルを作る時に、標準偏差は共通させるか？分けるか？
- 共通させる場合
  - $f(\theta) = f(\mu_1, \mu_2, sigma) = f(\mu_1)f(\mu_2)f(\sigma)$
  - $f(\theta|x) = f(\mu_1, \mu_2, \sigma|x1, x2) ∝ f(x1, x2| u1, u2, \sigma)f(u1, u2, \sigma)$
- 独立に定義する場合
  - $f(\theta) = f(\mu_1, \mu_2, \sigma) = f(\mu_1)f(\mu_2)f(\sigma1)f(\sigma2)$
  - $f(\theta|x) = f(\mu_1, \mu_2, \sigma|x1, x2) ∝ f(x1, x2| \mu_1, \mu_2, \sigma_1, \sigma_2)f(\mu_1, \mu_2, \sigma_1, \sigma_2)$
- かならずしも独立に定義したほうが良いわけではない。モデルの複雑性は増してしまうので。
生成量の例
- 平均の差
- 平均の差の効果量(estimated size)
  - $es = \frac{\mu_1 - \mu_2} {\sigma}$
- 非重複度(Cohen’s U) : u1が群2で何%点に相当するか
  - $U_3 = F(\mu_1| \mu_2, \sigma)$
  - 0.5 で2群が完全に重複していることを示す
- 優越率(pid) : 無作為に選んだ1方の群の測定値が、無作為に選んだもう一方の測定値を上回る確率
  - - (zは標準正規分布に従う変数)
  - 評価をするには、標準正規分布における $\frac{es}{\sqrt{2}}$ の分布関数を調べれば良い
- 閾上率(pic) : 無作為に選んだ一方の群の測定値と、もういっぽうの測定値の差がcより大きくなる確率
  - $pic = p(x^\ast_1 - x^\ast_2 > c)$
  - 評価をするには、標準正規分布における $\frac{\mu_1 - \mu_2 - c} {\sqrt{2}\sigma}$ の分布関数を調べればよい
モデルの選択
- 標準偏差は共通させたほうがよい？独立のほうがよい？
  - 基本的に独立の方が制約がなく、子細に利用可能
  - しかし、制約が正しい場合には、安定しやすい
- WAICの導入
  - 真の分布が得られたときの、未知データに対する平均対数尤度の推定値(豊田実践ベイズ本より)
    - 未知データに対する予測力
    - 小さい方が良い
  - - 詳細はhttp://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html
    - 実践ベイズモデリング -解析技法と認知モデル-も詳しかったです。
  - WAICはMCMCのChain中にある全サンプル値を使って未知データの予測を行い、過学習している分を補正する感じ？

章末問題

データの取得

とりあえずデータを読み込みます。今回は2群の対応のない場合の比較です。

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

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

# データの準備
a = np.array([56,55,55,62,54,63,47,58,56,56,57,52,53,50,50,57,57,55,60,65,53,43,60,51,52,
60,54,49,56,54,55,57,53,58,54,57,60,57,53,61,60,58,56,52,62,52,66,63,54,50])
b = np.array([33,37,59,41,42,61,46,25,32,35,55,44,45,41,33,61,46,16,48,34,27,37,28,31,32,
20,50,42,26,55,45,36,51,51,50,48,47,39,36,35,32,38,25,66,54,27,35,34,49,39])
logger = get_logger()

各種統計量の計算と可視化

各種統計量はpandas.describeを使えば大体出せますね。
分布はdistplot以外にも色々ありますが、matplotlibよりseabornの方が簡単に色々でるので使いました。

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.core.display import display
%matplotlib inline

df = pd.DataFrame({"a":a, "b":b})
display(df.describe()) # 分散は出ない

sns.distplot(a)
plt.show()
sns.distplot(b)
plt.show()

結果

	a	b
count	50	50
mean	55.76	40.38
std	4.57839	11.1427
min	43	16
25%	53	33
50%	56	39
75%	58	48
max	66	66

f:id:ajhjhaf:20170525001759p:plain

f:id:ajhjhaf:20170525010310p:plain

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

a群の方が平均値が15ぐらい高くなっている
中央値となる50%点もa群の方が高いが、なんといっても最小値が大幅に違う
この2つの群はともに正規分布っぽい感じにはなっている
分散は2群で結構違いそう。

みたいなものがわかります。よくある統計でしたら、t検定で2群の間には有意差があるということを出して終わりですね。今回は最初に書いた内容を使って、

差の大きさがどれぐらい幅があるのか(ちょっとだけの差でしたら、殆ど意味がないですよね)。
差は14ぐらいありそうである。それはどれぐらい確かなのか。
2つの分布はどれぐらい離れているのか。重複してる部分はどれぐらいあるのか。
群としてはaの方が大きくても、この群から適当に1つずつ選んだ時に、どれぐらいの確立でaの方の個体が大きな値になるのか

と言ったことを調べていくわけです。

Stanによる統計モデルの構築

Stanでモデルを作って、2群を比較します。教科書に従って、とりあえず不偏分散のモデルを作っています。aとbの分散値が先程の要約テーブルで結構違いましたので、等分散よりモデルとしては適当になりそうです(後で検証します)

import os
import pystan
import pickle

# Stanのモデルを読み込んでコンパイルする
stan_file = os.path.join("stan", "g2_ind.stan")
stan_file_c = os.path.join("stan", "g2_ind.pkl")
model = pystan.StanModel(file=stan_file)
with open(stan_file_c, "wb") as f:
    pickle.dump(model, f)

# Stanのモデル
data {
    int<lower=0> n_a ;
    int<lower=0> n_b ;
    real<lower=0> a[n_a] ;
    real<lower=0> b[n_b] ;
    real<lower=0> c_mu_diff ;
    real<lower=0> c_es ;
    real<lower=0> c_cohenu ;
    real<lower=0> c_pod ;
    real<lower=0> c_pbt ;
    real<lower=0> cdash_pbt ;
}

parameters {
    real<lower=0> mu_a ;
    real<lower=0> sigma_a ;
    real<lower=0> mu_b ;
    real<lower=0> sigma_b ;
}

model {
    a ~ normal(mu_a, sigma_a) ;
    b ~ normal(mu_b, sigma_b) ;
}

generated quantities {
    vector[n_a] log_lik ;
    real mu_diff ;
    real es ;
    real cohenu ;
    real pod ;
    real pbt ;
    int<lower=0, upper=1> prob_mu_diff_upper_0 ;
    int<lower=0, upper=1> prob_mu_diff_upper_c ;
    int<lower=0, upper=1> prob_es_upper_c ;
    int<lower=0, upper=1> prob_cohenu_upper_c ;
    int<lower=0, upper=1> prob_pod_upper_c ;
    int<lower=0, upper=1> prob_pbt_upper_cdash ;

    for(i in 1:n_a){
        log_lik[i] = normal_lpdf(a[i] | mu_a, sigma_a) + normal_lpdf(b[i] | mu_b, sigma_b) ;
    }
    mu_diff = mu_a - mu_b ;
    es = mu_diff / sigma_a ;
    cohenu = normal_cdf(mu_a, mu_b, sigma_b) ;
    pod = normal_cdf(mu_diff / sqrt(pow(sigma_a, 2) + pow(sigma_b, 2)), 0, 1) ;
    pbt = normal_cdf((mu_diff - c_pbt) / sqrt(pow(sigma_a, 2)  + pow(sigma_b, 2)), 0, 1) ;
    prob_mu_diff_upper_0 = mu_diff > 0 ? 1 : 0 ;
    prob_mu_diff_upper_c = mu_diff > c_mu_diff ? 1 : 0 ;
    prob_es_upper_c = es > c_es ? 1 : 0 ;
    prob_cohenu_upper_c = cohenu > c_cohenu ? 1 : 0 ;
    prob_pod_upper_c = pod > c_pod ? 1 : 0 ;
    prob_pbt_upper_cdash = pbt > cdash_pbt ? 1 : 0 ;
}

統計モデルによる事後分布のサンプリング

注意点は2つあります。まず1つは可視化をする部分です。 PyStanのモジュールをそのまま使ったサンプリング値の可視化では、サンプリング値が全て同値である場合(=分散が0の時)はtraceplotでエラーが発生してしまいます。これは可視化するために行っている、カーネル密度推定の行列計算で失敗するためです。今回、僕は自分でPyStanの一部を改造したものを利用しました。が、これは試してみたかっただけなので、こんな面倒なことをせず、summary_dfのsd columnを見て、0だったらtrace plotを書かなければ良いと思います。

もう1つはWAICの計算の部分です。計算をする際に、データ数分だけfor文を回して対数尤度を計算していますこの実装をしているものが検索してみると見つかるし、豊田本上級でもこの紹介がされています(定義的にも正しい？)。ただ今回の場合は、モデルが単純で50個のデータが共通の母数を持つ正規分布から生成されていので、対数尤度は1次元配列で計算し、waicの計算でsumをする部分が無くても大丈夫です。本の付録のモデルでは、こちらの簡潔版で計算されています。

import pandas as pd
import pickle
import pystan
import matplotlib
import os
import matplotlib.pyplot as plt
from IPython.core.display import display
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (10, 50)

# Stanで使うデータの用意
stan_data = {"n_a": a.size,
             "n_b": b.size,
             "a": a,
             "b": b,
             "c_mu_diff": 14, #標本平均の差を使ってみる
             "c_es": 3.0, # 標本効果量で、a群からみたもの
             "c_cohenu": 0.95, # a群から見た非重複度
             "c_pod": 0.95, # 優越率
             "c_pbt": 10, # 閾上率の基準値
             "cdash_pbt": 0.60} #閾上率
# 興味のあるパラメータの設定
par = ["mu_a",
       "sigma_a",
       "mu_b",
       "sigma_b",
       "log_lik",
       "mu_diff",
       "es",
       "cohenu",
       "pod",
       "pbt",
       "prob_mu_diff_upper_0",
       "prob_mu_diff_upper_c",
       "prob_es_upper_c",
       "prob_cohenu_upper_c",
       "prob_pod_upper_c",
       "prob_pbt_upper_cdash"]
prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975]

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

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


# 事後分布の表を取得
summary = fit.summary(pars=par, probs=prob)
summary_df = pd.DataFrame(summary["summary"],
                          index=summary["summary_rownames"],
                          columns=summary["summary_colnames"])
display(summary_df)

# 事後分布の可視化
fit.traceplot(par, {"prob_mu_diff_upper_0":np.int})
plt.show()

# WAICの計算
log_lik = fit.extract("log_lik")["log_lik"]
waic = -2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))
logger.info("WAICの値は{0}です".format(waic))

結果(表と分布)

	mean	se_mean	sd	2.5%	5%	25%	50%	75%	95%	97.5%	n_eff	Rhat
mu_a	55.7604	0.00210572	0.665888	54.4483	54.6698	55.3161	55.7601	56.2025	56.8535	57.0771	100000	0.99999
sigma_a	4.69936	0.00163133	0.492013	3.85845	3.97235	4.35244	4.65696	5.00083	5.56823	5.78077	90964	0.999993
mu_b	40.3768	0.00513303	1.62321	37.1566	37.7085	39.2942	40.3763	41.4601	43.0385	43.5743	100000	0.999979
sigma_b	11.4351	0.00377584	1.19403	9.39226	9.66662	10.5891	11.3316	12.1638	13.5548	14.0673	100000	1.00003
log_lik[0]	-6.04735	0.000490513	0.152051	-6.35677	-6.30426	-6.14733	-6.04291	-5.94225	-5.80452	-5.76237	96090	0.999988
log_lik[1]	-5.88977	0.000482951	0.146824	-6.18808	-6.13823	-5.98677	-5.88593	-5.78838	-5.65612	-5.61159	92424	1.00001
log_lik[2]	-7.21332	0.000987215	0.312185	-7.8894	-7.76195	-7.40848	-7.19029	-6.99224	-6.74368	-6.67145	100000	0.999982
log_lik[3]	-6.74241	0.000737068	0.233082	-7.24164	-7.15023	-6.89021	-6.72715	-6.57867	-6.38654	-6.32838	100000	0.999992
log_lik[4]	-5.91407	0.000486231	0.147301	-6.21372	-6.16309	-6.01166	-5.90981	-5.8126	-5.67868	-5.63577	91775	1.00002
log_lik[5]	-8.73421	0.00139258	0.440372	-9.66647	-9.49826	-9.01462	-8.7067	-8.42342	-8.06105	-7.95471	100000	0.999977
log_lik[6]	-7.74988	0.00123701	0.391176	-8.59949	-8.44025	-7.99428	-7.7188	-7.47338	-7.16121	-7.07202	100000	1.00001
log_lik[7]	-6.88161	0.000752406	0.237932	-7.39441	-7.29783	-7.03121	-6.86651	-6.71561	-6.51828	-6.45829	100000	0.999984
log_lik[8]	-6.10953	0.000492635	0.155785	-6.42726	-6.37336	-6.21194	-6.10435	-6.00181	-5.86187	-5.81865	100000	0.999984
log_lik[9]	-5.94667	0.00048357	0.148166	-6.24786	-6.19765	-6.04461	-5.9428	-5.8443	-5.70916	-5.66716	93881	0.999998
log_lik[10]	-6.71085	0.000707805	0.223828	-7.19011	-7.1025	-6.85252	-6.6971	-6.55379	-6.36867	-6.30991	100000	0.999989
log_lik[11]	-6.21359	0.000504989	0.159691	-6.54217	-6.48589	-6.31792	-6.20807	-6.10355	-5.96017	-5.91534	100000	1.00002
log_lik[12]	-6.09369	0.00049108	0.150984	-6.40111	-6.34936	-6.19326	-6.08984	-5.98977	-5.85288	-5.80928	94527	1.00002
log_lik[13]	-6.60844	0.000670976	0.212181	-7.05939	-6.97869	-6.7431	-6.59501	-6.461	-6.28327	-6.22806	100000	1.00001
log_lik[14]	-6.82174	0.000683128	0.216024	-7.28002	-7.19688	-6.95918	-6.80882	-6.67223	-6.48941	-6.43331	100000	1.00001
log_lik[15]	-7.5454	0.0011654	0.368531	-8.34713	-8.19487	-7.77589	-7.51674	-7.28485	-6.99185	-6.90637	100000	0.999982
log_lik[16]	-5.99181	0.000501445	0.14957	-6.29642	-6.24462	-6.09066	-5.9877	-5.8892	-5.7523	-5.71028	88970	1.00003
log_lik[17]	-8.19007	0.00156933	0.496266	-9.27191	-9.06673	-8.49926	-8.15016	-7.83863	-7.44239	-7.33214	100000	1
log_lik[18]	-6.48041	0.000548936	0.173589	-6.83737	-6.77563	-6.59425	-6.47343	-6.35979	-6.20713	-6.15985	100000	1.00001
log_lik[19]	-7.98662	0.00135638	0.428926	-8.91464	-8.74592	-8.25546	-7.95257	-7.68182	-7.34138	-7.24086	100000	0.999979
log_lik[20]	-6.71568	0.000647957	0.204902	-7.14689	-7.07082	-6.84698	-6.70495	-6.57247	-6.39928	-6.34525	100000	0.999987
log_lik[21]	-9.68219	0.00249132	0.787825	-11.3858	-11.0712	-10.1776	-9.6205	-9.12183	-8.49679	-8.32186	100000	1.00001
log_lik[22]	-6.85578	0.000652113	0.206216	-7.2915	-7.21196	-6.98743	-6.84574	-6.71229	-6.53602	-6.47974	100000	0.999978
log_lik[23]	-6.70807	0.000604833	0.191265	-7.10608	-7.03765	-6.83147	-6.69943	-6.57574	-6.40914	-6.35915	100000	1
log_lik[24]	-6.43882	0.000529426	0.167419	-6.78252	-6.724	-6.54813	-6.43286	-6.32297	-6.17362	-6.12905	100000	0.999994
log_lik[25]	-7.88994	0.0011714	0.37043	-8.69871	-8.54522	-8.12159	-7.86367	-7.6282	-7.33072	-7.24402	100000	0.999994
log_lik[26]	-6.26907	0.000533125	0.164099	-6.60602	-6.54939	-6.37576	-6.26344	-6.15559	-6.00966	-5.96361	94744	1.00001
log_lik[27]	-6.90996	0.000818656	0.258882	-7.4684	-7.36431	-7.0731	-6.8915	-6.72761	-6.51943	-6.45387	100000	1.00001
log_lik[28]	-6.64837	0.000691255	0.218594	-7.11348	-7.02718	-6.7864	-6.63571	-6.49597	-6.31147	-6.25681	100000	0.999982
log_lik[29]	-6.74747	0.000709018	0.224211	-7.22664	-7.13774	-6.88949	-6.73451	-6.59022	-6.40498	-6.34879	100000	0.999986
log_lik[30]	-5.92905	0.000495427	0.148163	-6.23007	-6.17892	-6.02743	-5.92463	-5.82661	-5.69268	-5.64909	89438	1.00002
log_lik[31]	-5.9427	0.000481788	0.147672	-6.24341	-6.19239	-6.04019	-5.93912	-5.84076	-5.70632	-5.66325	93948	1
log_lik[32]	-6.45466	0.000551094	0.174271	-6.81692	-6.7533	-6.56731	-6.4479	-6.33399	-6.1815	-6.13349	100000	1
log_lik[33]	-6.39362	0.00054583	0.172607	-6.75128	-6.68942	-6.50549	-6.38655	-6.27466	-6.12209	-6.07451	100000	1.00001
log_lik[34]	-6.26907	0.000533125	0.164099	-6.60602	-6.54939	-6.37576	-6.26344	-6.15559	-6.00966	-5.96361	94744	1.00001
log_lik[35]	-6.09634	0.000515248	0.154159	-6.40996	-6.35791	-6.19751	-6.09192	-5.99059	-5.85023	-5.80673	89517	1.00002
log_lik[36]	-6.42419	0.000541343	0.171188	-6.77728	-6.71547	-6.53626	-6.41756	-6.30555	-6.15428	-6.1082	100000	1.00001
log_lik[37]	-5.87456	0.000486567	0.147509	-6.1755	-6.12359	-5.97152	-5.87056	-5.77306	-5.63864	-5.59665	91908	1.00002
log_lik[38]	-6.08502	0.000474243	0.149969	-6.38969	-6.33788	-6.18427	-6.08081	-5.98184	-5.84684	-5.80406	100000	1
log_lik[39]	-6.58677	0.000616932	0.195091	-6.9986	-6.92179	-6.71234	-6.57743	-6.45114	-6.2843	-6.23202	100000	0.999988
log_lik[40]	-6.52813	0.000557474	0.176289	-6.89347	-6.82835	-6.64266	-6.52131	-6.40592	-6.25027	-6.20306	100000	0.999983
log_lik[41]	-5.97066	0.00048391	0.148823	-6.27342	-6.2218	-6.06889	-5.9667	-5.86849	-5.73262	-5.69075	94583	1.00001
log_lik[42]	-6.7658	0.000749145	0.2369	-7.27327	-7.17892	-6.91408	-6.75068	-6.60047	-6.40354	-6.34414	100000	0.999985
log_lik[43]	-8.75266	0.00173329	0.548115	-9.93939	-9.71752	-9.09476	-8.70939	-8.36307	-7.93373	-7.80486	100000	0.999986
log_lik[44]	-7.4732	0.000864603	0.273411	-8.05226	-7.94825	-7.64826	-7.45718	-7.28092	-7.05404	-6.98571	100000	0.99998
log_lik[45]	-6.86811	0.000670866	0.212146	-7.31424	-7.2334	-7.00434	-6.85678	-6.71946	-6.53968	-6.48454	100000	0.999993
log_lik[46]	-8.39545	0.00163259	0.516269	-9.51436	-9.31025	-8.71891	-8.35481	-8.02785	-7.61944	-7.50049	100000	0.999979
log_lik[47]	-7.21641	0.00090867	0.287347	-7.83597	-7.72231	-7.398	-7.19558	-7.01405	-6.78238	-6.71172	100000	0.99998
log_lik[48]	-6.19707	0.000517114	0.15826	-6.52008	-6.46738	-6.30066	-6.19209	-6.08782	-5.94499	-5.90124	93663	1.00001
log_lik[49]	-6.61441	0.000670758	0.212112	-7.0674	-6.98435	-6.74872	-6.60168	-6.46637	-6.28944	-6.23353	100000	1.00001
mu_diff	15.3836	0.00555407	1.75635	11.9336	12.5127	14.207	15.3843	16.5573	18.2745	18.853	100000	0.999983
es	3.30844	0.00160306	0.506931	2.36452	2.50846	2.95786	3.28843	3.63776	4.17218	4.35433	100000	0.99998
cohenu	0.90849	0.000107709	0.0335828	0.831913	0.847353	0.888497	0.912552	0.932937	0.955767	0.961685	97214	0.999988
pod	0.891195	0.00010542	0.0333368	0.816538	0.831072	0.87091	0.89468	0.915269	0.939324	0.945675	100000	0.999983
pbt	0.667754	0.000166733	0.0527255	0.560294	0.578451	0.632602	0.669254	0.704582	0.751648	0.766491	100000	0.999974
prob_mu_diff_upper_0	1	0	0	1	1	1	1	1	1	1	100000	nan
prob_mu_diff_upper_c	0.78579	0.00138494	0.410275	0	0	1	1	1	1	1	87758	0.999979
prob_es_upper_c	0.72119	0.00151356	0.448416	0	0	0	1	1	1	1	87773	0.999991
prob_cohenu_upper_c	0.0856	0.000957681	0.279774	0	0	0	0	0	1	1	85344	0.999981
prob_pod_upper_c	0.01434	0.00040641	0.118889	0	0	0	0	0	0	0	85576	0.999976
prob_pbt_upper_cdash	0.89648	0.00108438	0.304638	0	0	1	1	1	1	1	78924	1.00001

f:id:ajhjhaf:20170525001803p:plain

問題なく、パラメータ値が収束して計算されました(Rhat<1.1で判定しています)

WAIC = 683.147として計算されます。なお等分散モデルだとWAIC=716.444なので、実際に不偏分散モデルのほうが感覚としても、情報量としてもデータとモデルがフィットしていることがわかります。

仮説と検証

A群の平均値がB群より高い確率は？
- prob_mu_diff_upper_0のサンプリング結果より、100%の確率で、A群の方がB群より高いことがわかります。
AぐんとB群の平均値の差の点推定、95%信頼区間推定
- 15.384(1.756)[11.934, 18.853]となりました。
平均値の差の片側95%信頼区間推定の上限、下限はどのあたりか
- mu_diff > 12.513、あるいはmu_diff < 18.275の区間に95%の確率で入る
平均値の差が14より大きい確率は？
- prob_mu_diff_upper_cより、78.6%の確率で、14より大きいことがわかります。
A群から見たB群の効果量の点推定、95%信頼区間推定両側/片側
- 3.308(0.507)[2.365, 4.354]となりました。
- 片側区間なら、es > 2.508、あるいは es < 4.172の区間に95%の確率で入りました。
A群から見たB群の効果量が3より大きい確率は？
- prob_es_upper_cより、72.1%の確率で、3より大きいことがわかります。
A群の分布とB群の分布の重複する割合(非重複度)の点推定、95%信頼区間推定両側/片側
- 0.908(0.034)[0.832, 0.962]となりました。
- 片側区間なら、cohenu > 2.508、あるいはcohenu < 0.956の区間に95%の確率で入ります。
- これはB群から見た時に、典型的なA群の測定値は40.8%情報に存在することを示す
非重複度が0.95より大きい確率
- prob_cohenu_upper_cより、8.6% の確率で、0.95より
無作為に選んだ際に、A群の値がB群より大きくなる確率(優越率)はどのぐらいか
- 0.891(0.033)[0.817, 0.946]の確率で高くなることが推定されます。
優越率が0.95より大きい確率
- prob_pod_upper_cより 1.4%の確率で0.95より高くなることが推定されます。
無作為に選んだそれぞれの測定値の差が10より大きい確率(閾上率)
- 0.668(0.053)[0.560, 0.766]の確率で、大きいことが推定されます。
閾上率が0.60(60%)より大きい確率
- prob_pbt_upper_cdashより89.6%の確率で60%より高くなることが推定されます。

以上のように、MCMCによるサンプリングの試行結果を使うことで、様々な情報を抜き出すことが出来ました。次回は第4章の、対応のあるt検定です。

2017-05-23

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

ベイジアンモデリング

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

2章 MCMCと正規分布の推測

内容

解析的に事後分布を推定するのは難しい!
- ベイズの定理を利用しても、分母に位置する正則化定数の積分が複雑で解析的に解けない
事後分布に従う乱数を発生させる、という発想転換を行う(マルコフ連鎖モンテカルロ法; MCMC)
- 乱数群の結果を比較して、解析する
- HM法とGS法が具体的なアルゴリズムとしては代表的な存在
  - GS法は特定の事前分布であることが必要(やや恣意的)
  - HM法は自由な事前分布を使うことが出来る
- 事前分布の設定 : 常識的な知識から決める
- 乱数の発生 : 初期乱数は使用しない(burn in、あるいはwarm upと呼ばれる)
  - 値が登ったり下りたり(ドリフトした)形状が観測される場合には、好ましくない乱数生成となる
  - 複数chainを発生させて、全部で統一的な傾向が見られてない場合、好ましくない乱数生成が行われていると考えられる
    - 収束判定指標(Rcap) : 複数chainの分散, 1.1以下が理想
    - 有効標本数(neff) : 乱数の内、いくつが無関係な乱数として扱うことが出来るかを示す指標
    - (以上２つの指標の導出が気になった場合、基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門に結構丁寧に書いてあるので、そちらをお薦めします。)
事後分布の解釈をどう行うか？
- 点推定量
  - 事後期待値(EAP) = 平均値として扱うことが出来る
  - 事後中央値(MED) = 中央値として(ry
  - 事後確率最大値(MAP) = 最頻値(ry
  - 事後分散
  - 事後標準偏差
    - 与えられたデータを元に推定した分布自体の広がりがどれぐらいあるかを示したもの
  - 事後期待値の標準偏差(標準誤差; se)
    - EAPがどれぐらいの広がりを持つか示したもの
- 確信区間(信用区間)
  - 母数が存在する確率がx % である区間を確信区間とする
- 予測分布
  - 新たなデータが得られうる分布のこと
  - 事後予測分布
    - $f(x^*|x)$
    - 推定された母数の事後分布を用いて、未知データ $x^{*}$ を予測する統計モデルの平均を計算
    - $x^*(t) \sim f(θ(t))$
    - 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
  - 条件付き予測分布
    - 事後分布の特定の点推定量を利用する
    - 簡単だが、せっかく推定したモデルを有効活用出来ない
    - (これについては、その1で紹介した「はじめての統計データ分析」豊田秀樹のメモ - StatModeling Memorandumさんの方で問題点が指摘されています)
    - (が、豊田先生もこのような指摘を結構受けたのか、本に付録でついてくる追加Q&Aに解説が追記されています(後述))
- ベイズ的推測
  - 何らかの仮説を定めて、データから推定すること
  - 統計的仮説検定では帰無仮説に基づいたが、そんなことはしなくてよいのがベイズ統計
- 生成量
  - 推定した母数を利用して、各種統計量の分布も計算することが出来る
    - 母数より生成される値なので、生成量という
- 研究仮説が正しい確率
  - 真偽を判定する二値変数を利用して、研究仮説が正しい確率を評価することが出来る
  - 0 or 1しか使わないので、分散や信頼区間は使えない

条件付き予測分布の取扱について(2017-05-28追記)

ここまでで書いた概要部分でも説明しましたが、議論のされている部分なので補足説明です。このテキストではデータ $x$ が得られた状態で、次に得られるであろう未知データ $x^{\ast}$ についての予測方法を2つ紹介しています。

1つは事後予測分布です。事後予測分布の数式は $p_{pred}(x^\ast|x) = \int p(x^\ast|\theta)p(\theta|x)d\theta$ と記述されます。これを踏まえたやり方です。数式だけ見ると混乱しがちなのですが、MCMCで計算する時にはそんなに難しくありません。 MCMCの各サンプリングにおいて得られる $\theta$ が $p(\theta|x^\ast)$ を示していますし、この $\theta$ を構築したモデルに代入して得られる値が分布 $p_{pred}(x^\ast|x)$ からサンプリングされたうちの1つということになります。 MCMCが積分の代わりをしている、という認識ですね。

もう1つは条件付き予測分布です。これはMCMCの結果として得られる $\theta$ のEAPやMAPといった点推定値をモデルに代入して、予測分布を得るというものです。これは理解もし易いし、簡単です。分布の中で最もそれっぽいと理解出来た値をモデルに代入して予測をすることができます。著者の豊田先生としては、

理解が簡単
速度が速い(生成量や予測量として計算しなくて済む)
オンライン性(事後予測分布は過去データ $x$ が必要だが、条件付き予測分布は点推定値だけでよい)

という利点を挙げています。

しかしながらリンク先にもあるように、条件付き予測分布には色々問題点もありますし、何より個人的には「もったいない」という点が目立つと考えています。せっかく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

f:id:ajhjhaf:20170523020636p:plain

Rhat値が全パラメータで1.1以下になっていることや、Traceplotの様子からMCMCが収束していることが確認できると思います。というわけで推定値は一応問題がありません。実際に推定値を用いて推論を行うことが出来ます。

仮説と検証

平均値

このデータにおける平均値はいくつか
- 事後分布のEAPが40.660(0.650)であるため、40.660と推定される
このデータにおいて、平均値が95%の確率で分布する区間はどこか
- 事後分布の平均値の2.5%区間と97.5%信頼区間は39.383から41.947である
このデータにおいて平均値が片側で95%の信頼区間を持つ区間はどこか
- 事後分布の平均値の95%信頼区間より、41.735以下になる確率が95%である
- もしくは39.603以上になる確率が95%である

標準偏差

標準偏差の平均値、両側/片側95%信頼区間はどこか
- 6.523(0.468)、5.660から7.488、5.797以上あるいは7.321以下である

予測分布

新しい人からデータをとったときに、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

2017-05-22

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

ベイジアンモデリング

前置き

環境構築

まずPythonを実行する環境を作ります。一応メモとして書きますが、他にもPythonの導入の記事はたくさんありますので、そちらも参考ください。以下のものはOSXを想定しています。

自分は基本的にはpyenv + pyenv-virtualenvを利用しています。pyenvで切り替えつつ、学習する本の内容に合わせてpyenv-virtualenvで環境を作っています。導入は各リポジトリのページに書いてあるとおりです。

github.com github.com

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章データの整理とベイズの定理

内容

一応箇条書きで書いていますが、メモです。ちゃんと勉強したいときは、本文を参考にしてくださいね。

分布の紹介
- 正規分布
- 一様分布
データの解釈
- なぜ正規分布にあてはめる?
  - 観測データの実際の分布は、より複雑な謎の分布に従うと予測される
  - 正規分布はパラメータが少なく、扱いが容易
```
 * データが分布にあてはめられることで、分布の数式を元にしたマニピュレーションが可能
```
- どのようにあてはめる?
  - 頻度主義 : 標本平均、標本偏差を推定して、あてはめる
  - ベイズ統計 : データからパラメータの分布を推定し、その曖昧さを踏まえて分布を推定
    - ベイズの定理より、データが得られた時のパラメータ分布が得られる
```
      * ベイズの定理のために必要なもの
```
      - 尤度
      - 事前分布
        
        無情報事前分布の利用
      - 正則化定数(扱いが難しい!)
      - そこでMCMCですよ。

章末問題

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だとめんどくさかったりしますね。特定の幅長での度数表や、最頻値について簡単に求めることが出来る方法を知っていらっしゃる方が居たら、是非コメントを頂けると助かります！

2017-05-22

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

ベイジアンモデリング

このブログの開設目的の、

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

作者: 豊田秀樹
出版社/メーカー: 朝倉書店
発売日: 2016/06/02
メディア: 単行本（ソフトカバー）
この商品を含むブログ (11件) を見る

についての記事インデックスです。

本の紹介

この本は、昨今はびこる統計データ分析についてベイジアンの知見から解説した本です。他の方の感想では、 statmodeling.hatenablog.com

d.hatena.ne.jp

のようなものがあります。個人的に一冊まるまるよんだ感想としては、次のような印象を持ちました。

良い点

MCMCについての細かい知識は要らない

ベイズの定理やMCMCを利用した多くの解析の本では、数式がガッツリ出てきて頭を抱える方も多いと思います。しかし、この本ではStan言語を利用して解析を行っていきます。そのため細かい数式の変形や、推論の部分についてはStanがなんとかしてくれますので、本にも書いてありませんし、やる必要もありません。

Research Questionベースの解析

実際にデータ分析をする際には、様々なデータ提供者からの要望に答える形で解析をすることが多いと思います。この本では、そのような問いをResearch Questionと呼び、各章で該当するResearch Questionを最初に挙げて、解決する形で進めていきます。結果として、どのような形で今やっている解析が活きるものなのか、具体的にイメージしながら進める事ができます。

ベイズベースでのデータ解析への「入門」

豊田先生が後書きにも記していますが、普通データ解析や統計学を勉強するときには、頻度主義に基づく統計的仮説検定を勉強します(僕もそうしました)。一方でベイズの定理を基づくベイジアン統計学は、最近計算機処理が高速化したことやデータ量が多くなっていることが原因で、頻繁に使われています。しかしベイジアン統計学は、頻度主義とはぜんぜん違うため、入門が難しいと僕も思います。その中で本書は、統計的仮説検定の手法でやるような内容をベイジアン統計学で行った場合にどうなるかを解説しています。その点で、最初に取り組む際に何が利点となるのか分かりやすいと思います。

答えが用意されている

最近は多いですが、解析に使っているスクリプトがR + Stanで用意されています。

朝倉書店｜はじめての統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

更に、各章末問題は巻末に答えがかなり丁寧に書いてあります。大学の教科書もそうですが、答えがない場合自分が間違っていても気が付かないことがあります。答えがあるのは非常に助かりますね。

悪い点

入門書ではない

幾つもの指摘がありましたが、入門書ではありません。たとえば、

Stanの文法や使い方について丸々省略されている
統計モデリングの必要性が良くわからない
ベイズの定理の導出なども駆け足気味

といった点が目立ちました。1つ目、2つ目については久保先生の通称緑本、ベイズの定理について入門レベルから丁寧に勉強したい場合には超入門が良かったと思います。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

図解・ベイズ統計「超」入門あいまいなデータから未来を予測する技術 (サイエンス・アイ新書)

作者: 涌井貞美
出版社/メーカー: SBクリエイティブ
発売日: 2013/12/17
メディア: 新書
この商品を含むブログ (15件) を見る

なにをするのか

以降からは、各章の簡単な解説をした上で、Pythonで章末問題を解いていきたいと思います。 Stan言語については、PythonのインターフェイスであるPyStanを利用します。といっても、Stan言語で書かれている部分についてはRと共通して使うことが可能なので、PyStanを使うためにどのような小手先のテクニックがあるかを書いていくつもりです。

2017-05-22

ブログ始めました。

自己紹介

ブログ始めました。

普段は次世代シーケンサー(NGS)を用いた遺伝統計学や、ゲノム配列から得られる知見の解析、それらのツール・DB開発などを行っています。

ブログの目的

自身で勉強したことを内部のWikiに書いていたのですが、もっと外部に公開してマサカリを投げてもらったり、何かしらの役に立つことを期待してブログを開設します。基本的には、

Pythonによるモデリング、機械学習
Javascriptを用いたデータ可視化(D3.js)、アプリケーション開発
次世代シーケンサーから得られるデータの解析について

を主に書いていきたいと思います。特にPythonについては日本語知識があまり転がっていないので、Rや海外サイトから輸入することが多いかもしれません。 3番目のNGS解析については、需要がなさそうなのであまり書かない気がします笑

それでは宜しくお願いします。Amazon.co.jpアソシエイトです。

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

内容

章末問題

データの取得

各種統計量の計算と可視化

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

番外編 3章の等分散モデルと変分ベイズによる推定

等分散モデルによる差の推定

変分ベイズによる推定

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

内容

章末問題

データの取得

各種統計量の計算と可視化

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

仮説と検証

2章 MCMCと正規分布の推測

内容

条件付き予測分布の取扱について(2017-05-28追記)

章末問題

データの準備

Stanによる統計モデルの構築

統計モデルによる事後分布のサンプリング

仮説と検証

平均値

予測分布

生成量

研究仮説が正しい確率

前置き

環境構築

1章 データの整理とベイズの定理

内容

章末問題

データの準備

度数分布表

ヒストグラムの作成

統計量

正規分部の扱い

目次

本の紹介

良い点

MCMCについての細かい知識は要らない

Research Questionベースの解析

ベイズベースでのデータ解析への「入門」

答えが用意されている

悪い点

入門書ではない

なにをするのか

自己紹介

ブログの目的

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

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

1章データの整理とベイズの定理