Easy to type

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

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

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

2章 MCMC正規分布の推測

内容

  • 解析的に事後分布を推定するのは難しい!
    • ベイズの定理を利用しても、分母に位置する正則化定数の積分が複雑で解析的に解けない
  • 事後分布に従う乱数を発生させる、という発想転換を行う(マルコフ連鎖モンテカルロ法; MCMC)

    • 乱数群の結果を比較して、解析する
    • HM法とGS法が具体的なアルゴリズムとしては代表的な存在
      • GS法は特定の事前分布であることが必要(やや恣意的)
      • HM法は自由な事前分布を使うことが出来る
    • 事前分布の設定 : 常識的な知識から決める
    • 乱数の発生 : 初期乱数は使用しない(burn in、あるいはwarm upと呼ばれる)
      • 値が登ったり下りたり(ドリフトした)形状が観測される場合には、好ましくない乱数生成となる
      • 複数chainを発生させて、全部で統一的な傾向が見られてない場合、好ましくない乱数生成が行われていると考えられる
  • 事後分布の解釈をどう行うか?

    • 点推定量
      • 事後期待値(EAP) = 平均値として扱うことが出来る
      • 事後中央値(MED) = 中央値として(ry
      • 事後確率最大値(MAP) = 最頻値(ry
      • 事後分散
      • 事後標準偏差
        • 与えられたデータを元に推定した分布自体の広がりがどれぐらいあるかを示したもの
      • 事後期待値の標準偏差(標準誤差; se)
        • EAPがどれぐらいの広がりを持つか示したもの
    • 確信区間(信用区間)
      • 母数が存在する確率がx % である区間を確信区間とする
    • 予測分布
      • 新たなデータが得られうる分布のこと
      • 事後予測分布
        • f(x^*|x)
        • 推定された母数の事後分布を用いて、未知データx^{*}を予測する統計モデルの平均を計算
        • x^*(t) \sim f(θ(t))
        • 事後分布から母数をサンプリングし、それを統計モデルにあてはめて予測をする
      • 条件付き予測分布
    • ベイズ的推測
      • 何らかの仮説を定めて、データから推定すること
      • 統計的仮説検定では帰無仮説に基づいたが、そんなことはしなくてよいのがベイズ統計
    • 生成量
      • 推定した母数を利用して、各種統計量の分布も計算することが出来る
        • 母数より生成される値なので、生成量という
    • 研究仮説が正しい確率
      • 真偽を判定する二値変数を利用して、研究仮説が正しい確率を評価することが出来る
      • 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の各サンプリングにおいて得られる\thetap(\theta|x^\ast)を示していますし、この\thetaを構築したモデルに代入して得られる値が分布p_{pred}(x^\ast|x)からサンプリングされたうちの1つということになります。 MCMC積分の代わりをしている、という認識ですね。

もう1つは条件付き予測分布です。これはMCMCの結果として得られる \thetaEAPや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