MCMC

MCMCを実装するライブラリはいくつかあるようです。

PyMC3:Theano
PyMC4:TensorFlow
Pyro:PyTorch
NumPyro:JAX

手元の環境で動いたのはPyroのみ。環境構築に関してはシビアなようです。

web上のサンプルを見てみますと、書き方はほぼPyMC3と同じでした。動かすのは簡単（traceplot がないのは残念。seaborn で書きました）。

が、写経し終わると理解できていない箇所が露わに。

 y_model = pyro.sample('y_model', dist.Normal(a*x + b, sd), obs=y)

ここの obs の計算内容がわかりません。ここがキモなのはわかっているのですが、マニュアルに詳細が書かれていませんし、ソースも追えません。尤度×事前確率にどのようにかかわっているのでしょう？
おそらく、理解している方々が必要として使われているので、学び方が逆なのでしょうね。

ま、他の部分は理解できていることが分かりました。メモしておきましょう。

 **************************************
20201123 MCMC メモ
 
データ駆動科学

• データを出発点と考える=因果律を遡る
• 従来：変形係数を確定(モデルを確定)→応力、ひずみがノイズを含み確率的分布（→安全率を含めた設計）
• データ駆動：応力、ひずみが確定→変形係数が確率的分布（逆解析的発想）
• データｘ,y（ひずみ、応力）が得られたとき、パラメータa（変形係数）が取り得る確率分布（事後確率）を求めたい。
• 必要なのは、①パラメータaを与えたときにyが取り得る確率分布(尤度)、②aの確立分布（事前確立）

取得データD[x,y]

• 確定x：ひずみx
• 確定y：応力y

予測分布モデル　（例）物理モデル+ノイズ
　物理モデル y=ax+b

• 未確定a:平均μa,標準偏差σaとして設定(事前確率P(a)=N(μa,σa2)を仮定)
• 未確定b:平均μb,標準偏差σbとして設定(事前確率P(b)=N(μb,σb2)を仮定)
• 未確定σn:一様分布(U([n])を仮定)
  ←いずれもデータDを見て取りうる範囲を仮定
  

　ノイズ

• 重畳ノイズを正規分布として設定（尤度P(y_model)=N(ax+b,σn2)を仮定
  

予測分布モデルから θ[a:変形係数、b：切片、sd:標準偏差]の事後確率 p(θ|D)を求める。

• p(θ|D) = p(D|θ)p(θ)/p(D) ∝ p(D|θ)p(θ)
  
• 事後確率 ∝ 尤度×事前確率
• 事前確率、尤度に正規分布を仮定した場合、事後分布を解析的に求めることが可能（共役事前分布）。
• しかし、一般的には計算が困難。→MCMC法などで帰納的に事後確率（最大）を求める。

マルコフ連鎖モンテカルロ法（MCMC法）

• モンテカルロ法：乱数を使ったシミュレーション手法
• マルコフ連鎖：現在の状態が、前時刻の状態のみに依存するモデル→メトロポリス法などの利用

1. θ0を仮定
2. 乱数εを加え、θ1を作成。
3. 確率密度の比ｒを計算r=p(D|θ1)p(θ1)/ p(D|θ0)p(θ0)
   ノイズを正規分布と仮定した場合、p(D|θ) ∝ exp(-NE(θ)/σ2)
   あるθ0、θ1における exp(-NE(θ)/σ2)p(θ) を求める
   ←obs=D を利用 ※点での推定
   
4. 確率密度の比ｒを評価
   r>1でθ1を採用
   r<=1かつ乱数U（1,0）＜ｒでθ1を採用、U>=ｒでx0を採用
5. １～４の繰り返し←exp(-NE(θ)/σ2)p(θ)の分布形状を推定

 

他手法との違い

• 最小二乗法:残差の平方和 E(θ) を最小化
  
• 最尤推定:a,b,σについて導関数=0の位置（パラメータの値）を求める(global minimum とは限らない)
  
• ベイズ推定:パラメータを確率的に扱う→分布(形)を求める 

******************************************

20201129

MAP推定：パラメータ値の同定のみの場合（パラメータの分布を求めない）

• 作成した予測分布モデル＋観測データ（obs）に対し、MAP（maximum a posterior :最大事後確率）推定でパラメータを同定
• pyro.condition で同定パラメータをモデルに反映

 

20201202

変分推論：分布も推定

• svi = SVI(model, guide, optimizer, loss=Trace_ELBO())でモデル構造を作成
  
• step method でデータDをモデルに渡して調整←obsの役目
• predictive で予測分布を作成（モンテカルロ近似）
  

概ね、obs の役割について理解しました。
数式は理解できていません。残念ながら後回しです。

 

20201204

数式もOK（おそらく）。

これで MCMC（メトロポリス法利用）は理解できたと思います。