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