StanとRでベイズ統計モデリング, ４章の計算をしてみた

概要

StanとRでベイズ統計モデリングの4章のベイズ信頼区間と予測区間をpythonを自分で実装してみました。

github.com

実験したコードはこちらです。

動機

ベイズ統計を実務で使うためにStanとRでベイズ統計モデリングを読んでいます。

StanとRでベイズ統計モデリング (Wonderful R)

この本の4章では，推計統計(頻度論)とベイズ統計の両方のアプローチで信頼区間と予測区間を計算してました。
これらの実装はRで行われていたため，今回はpythonを使ってこれらの計算の描写を行うことを目標にしました。*1
また，統計の定義のややこしい単語を整理しておきます。

推計統計とベイズ統計のアプローチの違い

自分の記憶の整理のために書いているので誤りがあれば指摘してください。
訂正します。

推定統計

図1: 母集団と標本集団

推計統計では，解析対象は「母集団からランダムに抽出された標本」です。(図1)
母集団には全てのデータが存在しており，それらの平均値や分散はただ一つに決定されます。
一方で，標本は，母集団から抽出されるたびに異なるものとなるため，平均や分散は確率的なものになります。

図2:標本集団の推定結果

図2の青色の母集団の正規分布で，黒色が標本集団から抽出した正規分布であるとします。
このように，実際の母集団の結果と，標本集団から推定した結果は微妙に異なることはありえます。
例えば，標本数Nが小さいときの推定した分散は，母集団の分散より小さくなることが計算でも証明されています。

参考 staff.aist.go.jp

N数が十分大きければ，母集団の平均と分散は推定値と完全に一致します。

ベイズ統計

ベイズ統計では，観測されたデータが全てだと視点からスタートします。
推定統計では，標本を抽出するたびに確率分布は異なりますが，ベイズ統計では図1右側の標本集団が全データです。
そのため，推定統計では母集団の平均値や分散はただ一つの値だったのですが，ベイズ統計ではこれを確率的に推定しようという考え方になります。

信頼区間と予測区間

この2つの差を途中まで理解していなかったので記載しておきます。

参考

おっと危ない：信頼区間と予測区間を混同しちゃダメ - Take a Risk：林岳彦の研究メモ

予測区間95%というのは，新しいデータは95%の確率でその区間の中に存在することを指します。*2
一方で，信頼区間は推計統計とベイズ統計で意味合いが異なります。
我々が解釈しやすいのは，ベイズ統計であり推計統計ではないことを覚えておけば

推計統計における信頼区間

推定統計において信頼区間の解釈は，「繰り返し標本を抽出したときに，その抽出した標本の信頼区間の中に母集団の平均値が入っている確率」になります。
例えば，信頼区間95%であるときには，100回標本抽出を繰り返すと95回は推定した信頼区間の中に母平均が含まれます。

下のURLの画像がわかりやすかったです。

参考 bellcurve.jp

　ベイズ統計における信頼区間

これは単純に，母集団の平均値が信頼区間に入っている確率です。
解釈がしやすいですね。

計算と可視化を行ってみた。

推定統計

github.com

コードはこちらになります。

計算はstatsmodelに行わせました。
線形回帰を行うところまでは知っていたのですが，信頼区間と予測区間は，

result = sm.OLS(Y, X).fit()
result95 = result.get_prediction().summary_frame(alpha=0.05)

としてpandasのdataframe書き出すことができます。はじめwls_prediction_stdを使っていたのですが，これは予測区間しか出すことができませんでした。

また，区間を描写するのに，matplotlibのfill_between関数を使いました。大変便利です。

図3: 推定統計による信頼区間と予測区間。区間が細くなっている方が信頼区間

結果として，得られたのが図3になります。これは本のものと対応できていることが確認できました。

ベイズ推定

github.com

stanの部分に関しては，

StanとRでベイズ統計モデリング (Wonderful R)

にかかれています。pystanのAPIに関して言えば，stanファイルに記載する部分をpythonの文字列としてダイレクトに書き込み，

stanmodel='''
data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real<lower=0> sigma;
} 

model {
    for (n in 1:N){
    Y[n] ~ normal(a + b*X[n], sigma);
    }
}
'''

コンパイルすることで，

sm = pystan.StanModel(model_code=stanmodel)

実行できます。今回は初めてということで，generated quantitiesを使わずに，python側で正規分布を生成しています。

結果として，

図4:ベイズ信頼区間と予測区間

図4が得られ，無事に区間推定を行うことができました。

まとめ

pythonを使ってStanとRでベイズ統計モデリングの４章の計算を行ってみました。
信頼区間と予測区間の定義をまとめました。

*1:意外なことに観測範囲内ではあんまり誰もやってませんでした。

*2:いろいろ探ったのですが，推計統計もベイズ統計も同じように思えました。誤っていれば教えてください。