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 メモ
データ駆動科学
- データを出発点と考える=因果律を遡る
- 従来:変形係数を確定(モデルを確定)→応力、ひずみがノイズを含み確率的分布(→安全率を含めた設計)
- データ駆動:応力、ひずみが確定→変形係数が確率的分布(逆解析的発想)
- データx,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法)
- モンテカルロ法:乱数を使ったシミュレーション手法
- マルコフ連鎖:現在の状態が、前時刻の状態のみに依存するモデル→メトロポリス法などの利用
- θ0を仮定
- 乱数εを加え、θ1を作成。
- 確率密度の比rを計算r=p(D|θ1)p(θ1)/ p(D|θ0)p(θ0)
ノイズを正規分布と仮定した場合、p(D|θ) ∝ exp(-NE(θ)/σ2)
あるθ0、θ1における exp(-NE(θ)/σ2)p(θ) を求める
←obs=D を利用 ※点での推定 - 確率密度の比rを評価
r>1でθ1を採用
r<=1かつ乱数U(1,0)<rでθ1を採用、U>=rでx0を採用 - 1~4の繰り返し←exp(-NE(θ)/σ2)p(θ)の分布形状を推定
他手法との違い
- 最小二乗法:残差の平方和 E(θ) を最小化
- 最尤推定:a,b,σについて導関数=0の位置(パラメータの値)を求める(global minimum とは限らない)
- ベイズ推定:パラメータを確率的に扱う→分布(形)を求める
- 作成した予測分布モデル+観測データ(obs)に対し、MAP(maximum a posterior :最大事後確率)推定でパラメータを同定
- pyro.condition で同定パラメータをモデルに反映
- svi = SVI(model, guide, optimizer, loss=Trace_ELBO())でモデル構造を作成
- step method でデータDをモデルに渡して調整←obsの役目
- predictive で予測分布を作成(モンテカルロ近似)
数式は理解できていません。残念ながら後回しです。