Rでベイズ統計学の学習記録

【学習動機】

昔、MCMCやらなんやらは学んで、少し実装したことがあったが、

色々忘れていることも多いし、Rで手軽にできるようになっておきたいと思って、

J.アルバート『Rで学ぶベイズ統計学入門』

http://amzn.asia/d/hbKIgqd

をやっていくことにした。

 

【学習記録】

<総括>

LearnBayesパッケージにデータセットが豊富。

パッケージの重要な関数の作り方もあるのでわかりやすい。

練習問題。答えがあればなおよし。

密度や分布、周辺尤度あたりの言葉の使い方が気になった。和訳の問題かも。

やはり、BDA3を読むべき。

 

c1:R入門

LearnBayesパッケージ。

基本的な使い方の紹介。

 

c2:ベイズ的思考への誘い

・割合pのベイズ推定

事前分布の設定の仕方。

  • 離散事前分布の利用

LearnBayesのpdisc関数で計算。p discreteの意味か。

  • 連続事前分布(ベータ事前分布)の利用

共役の事前分布で、簡単に計算可能。

事後密度をカバーする区間でpの値のグリッドを選定する

このグリッド上で尤度と事前密度の積を計算する

それぞれの積を、それらの合計で割って正規化。これによりグリッド上の離散確率分布によって事後密度を近似している

Sample関数を使って、この離散分布からランダムな復元標本を抽出。

こうして得られたシミュレーション標本は事後分布からの標本を近似している。

 

予測

  • 離散事前分布

LearnBayesのpdiscp関数で予測密度計算。

  • ベータ事前分布

LearnBayesのpbetap関数で予測密度計算。

 

任意の事前分布から予測密度を計算できる便利な方法の一つがシミュレーションによるもの。

 

c3、c4 1パラメータモデル、複数パラメータモデル

 

c5:ベイズ計算入門

前2章では事後分布を要約する2タイプの方策。

ひとつは、サンプリング密度が指数分布族など周知の関数形ならば、パラメータを直接シミュレーションできる。もうひとつは周知の関数形でない場合、緻密な点のグリッド上で事後密度の値を計算して、これらの値が集中する離散事後密度によって、連続形の事後密度を近似できる(ブルートフォースと呼ばれる)。

本章では別のアプローチ。一般的なアプローチの一つが、モード周辺で事後分布の挙動に基づく方法で、多変量正規近似となり、最初の近似として役に立つ。続いて、事後分布の要約の計算にシミュレーションを利用する方法。事後分布が標準的な関数形ではない場合、適切に選ばれた提案密度の棄却サンプリング。重点サンプリングとSIR法のアルゴリズム積分計算し、一般の事後分布からシミュレーション。SIR法のアルゴリズムは、事前分布と尤度関数の変更に対して事後分布が敏感か調べられる。

 

事後分布の要約の多くは積分で表現できる。関数h(θ)の事後平均、h(θ)が集合Aに入る事後確率、関心のあるパラメータの周辺密度。これらは積分計算で求める場合、数値的に評価するが求積法は問題が限られる。

 

パラメータを再定式化して、近似しやすい形にする。

近似の一つは、事後モードにもとづく近似。多変量事後分布の要約を行う方法の一つが、そのモード周辺の密度の挙動を調べること

モード周りで多変量正規分布で近似するためにはモードを求める必要がある。

モードを求めるための汎用的な最適化アルゴリズムニュートン法やネルダー・ミード法。

 

棄却サンプリング、重点サンプリング、サンプリング重点リサンプリング。

小澄英男『ベイズ計算統計学

http://amzn.asia/d/6LIzKmB

C.P.ロバート、G.カセーラ『Rによるモンテカルロ法入門』

http://amzn.asia/d/e72cRHI

あたりで勉強しなおしたい。

c6:マルコフ連鎖モンテカルロ

棄却サンプリングは任意の事後分布からシミュレーションする一般的な方法であるが、適切な提案分布を構成する必要があるので、適用が困難。重点サンプリングとSIRもまた汎用的なアルゴリズムであるが、やはり提案分布を必要とし、特に高次元の問題の場合、難しい。そこで、MCMCによって事後分布を要約する方法。

 

c7:階層モデリング

 

心臓移植の死亡率データ。94の病院の心臓移植手術30日以内の死亡数データ。

死亡数y, 曝露数e。曝露1件あたりの死亡率λを推定したい。

すべての病院についての真の死亡率を同時推定したい。三つの方法。

一つ目の選択肢は、個々の死亡率で推定。 y_1/e_1,...,y_{94}/e_{94}

これは曝露数の少ない病院では不正確。ばらつきが大きい。

二つ目の選択肢は、真の死亡率はすべての病院で等しいという仮定で

\Sigma{y_j}/\Sigma{e_j}でプール推定。

三つ目の選択しは、一つ目と二つ目の複合。

(1-\lambda) y_i/e_i + \lambda \Sigma{y_j}/\Sigma{e_j}で、これは個別の推定値をプールされた推定値に縮める縮約推定値。

 

ベイズ分析では、モデルの仮定に変更を加えた場合の推論の感度を評価することが重要。

これにはサンプリング密度と事前密度についての仮定も含まれる。

5章のSIRは、ある事後分布からをシミュレーションした標本を、新しい分布に変換する便利な方法。

 

c8:モデル比較

あるパラメータについて二つの仮説を比較する方法として、ベイズファクター。

二つのベイズモデルの比較を行う。これは事前密度とサンプリング密度の選択に関わる。

この場合、ベイズファクターは二つのモデルの周辺密度の比で表される。

ベイズファクターは仮説の事前オッズに対する事後オッズの比で表される。

事前オッズ比は\pi_0 / \pi_1 = P(\theta_0 \in \Theta_0) / P(\theta_1 \in \Theta_1)

= \int_{\Theta_0}g(\theta)d\theta /  \int_{\Theta_1}g(\theta) d\theta

( g(\theta)は正則事前密度、g(\theta|y) \propto L(\theta)g(\theta)は事後密度 )

事後オッズは

p_0 / p_1 = P(\theta_0 \in \Theta_0 | y) / P(\theta_1 \in \Theta_1 | y)

= \int_{\Theta_0}g(\theta | y)d\theta / \int_{\Theta_1}g(\theta | y) d\theta

ベイズファクターは  BF = 事後オッズ/事前オッズ = (p_0/p_1) / (\pi_0/\pi_1)

統計量BFはデータによる証拠が仮説H0を支持する度合を示す基準。

仮説H0の事後確率はベイズファクターと仮説の事前確率の関数で表される。

 p_0 = \pi_0 BF / (\pi_0 BF + 1 -\pi_0)

 

仮説を比較するベイズの方法は、二つのモデル比較に一般化できる。

仮にyをデータベクトル、θをパラメータとすれば、ベイズモデルはサンプリング密度と事前密度を指定して構成される。

このモデルの周辺尤度は  m(y) = \int f(y|\theta)g(\theta)d\theta

M0とM1二つのベイズモデルを比較したいとする。

パラメータθの定義はモデル間で異なることもありうる。

このときモデルM0を支持するベイズファクターは、二つのモデルそれぞれに対するデータの周辺尤度の比になる。

 BF=m_0(y)/m_1(y)

とがモデルとそれぞれの事前確率を表すとすると、モデルの事前確率は以下で与えられる。

 P(M_0|y) = \pi_0 BF/(\pi_0 BF + \pi_1)

周辺尤度を近似する簡単な方法の一つがラプラス法。

\hat{\theta}が事後モードで、が対数事後密度のヘシアン(行列の2階微分)とする。このとき周辺尤度は以下で近似される。

ここでdはパラメータの数。

  

ベイズファクターの計算。

 

c9:回帰モデル

 

回帰モデル。

Zellnerのg事前分布によるモデル選択。

 

c10:ギブスサンプリング

これは別で学習するのでスルー

 

c11:RでWinBUGSを使うインターフェイス

BUGS言語、好きでないのでスルー

 

【学習予定】

『 StanとRでベイズ統計モデリング

www.amazon.co.jp

『実践ベイズモデリング

www.amazon.co.jp

あたりで、モデリングを習得していきたい。