前回
開発環境
ぱっと思いつく選択肢としては以下のようなものがありました。そもそも以下のような理由からDockerを使用するのが良いかなと思っていたこともあり、はじめに(3)の方法で進めました。しかし、途中でエラーにはまってしまい、(2)の方法に変えました。
Dockerである理由
- 既存のローカルの環境（Windowsのバージョンやエディション、pythonの環境など）との齟齬による、「環境構築できない」、「特定の関数の実行できない」等のエラーやコンフリクトを防げる。
- 手元のデータを使いやすい。
- ローカルの環境をきれいに保てる。
選択肢
-
PyMC公式のインストール方法（ローカル環境にAnacondaをインストールし、Anaconda経由でPyMCをインストール）に従う。
-
condaコマンドが使えるDockerコンテナ（Anaconda入りイメージ）を探し、コンテナ内においてPyMC公式のインストール方法に従う。
-
PyMC公式のContributingの為のリポジトリをクローンして、その中に用意されているPyMC公式のDockerコンテナを立ち上げる
-
Google Colabを利用する
例題の推論結果について
事前予測分布のチェック
以下は、例題（ケーススタディ）と同じコードです。
数理モデルの設計として、事前分布の設定が妥当そうかどうかを可視化してチェックします。あまりにお門違いなデータを出力したりしていないかどうかを見極めます。
with test_score_model:
prior_samples = pm.sample_prior_predictive(100)
az.plot_dist(
test_scores["score"].values,
kind="hist",
color="C1",
hist_kwargs={"alpha": 0.6},
label="observed",)
az.plot_dist(
prior_samples.prior_predictive["scores"],
kind="hist",
hist_kwargs={"alpha": 0.6,},
label="simulated",)
plt.xticks(rotation=45);
<出力結果>
縦軸はデフォルトではの確率密度です。縦軸をサンプル数にするには以下のように引数を変更します。
az.plot_dist(
test_scores["score"].values,
kind="hist",
color="C1",
- hist_kwargs={"alpha"},
+ hist_kwargs={"alpha": 0.6,"density": False},
label="observed",
)
まず、pm.sample_prior_predictive(100)により、事前分布からサンプリングを実行します。事前分布からサンプリングした後に、サンプリング値と入力変数をモデルに代入して順方向に計算しています。
一般的な理論式と実装上の式を調べてみると、以下のような式を多く見かけるかと思います。しかし本実装は、モンテカルロ法的にパラメータに対し多くのサンプルセットから平均を取り、予測精度を上げた予測値を出してプロットしているわけではなく、サンプリングしたパラメータを全セット分プロットしている実装になっています。
・事前予測分布（理論式）
\begin{aligned}
&p(y) = \int p(y \mid \theta) \pi(\theta) \, d\theta \\
\end{aligned}
・事前予測値（実装可能な式）
\begin{aligned}
&p(y) \approx \frac{1}{N} \sum_{i=1}^N p(y \mid \theta^{(i)}), \quad \theta^{(i)} \sim \pi(\theta) \\
\end{aligned}
テストデータは101人分ありますので、パラメータを100セットサンプリングした今回の場合、101×100の出力（score）が生成されます。また、試しにpm.sample_prior_predictive(100)の引数を1にしてプロットすると、100データ予測されていたので、この理解であっていると思います。
<入力データ（説明変数）>
ベイズ推定を理論から学び、理論上の演習問題などを解いた後に、実際の数式とコードと推論結果を追いながら理論を理解しようとするとサンプリング後に何が求まっているのか混乱しそうになります。しかし、当然のことですが、MCMCによるサンプリング実行後に何か値が求まったわけではないのがベイズ推論です。 MCMCによりサンプル列が取得できただけ であり、そのサンプル列を描画しているに過ぎません。そのサンプル列から、所望の分布（予測分布など）や値（期待値など）などを計算して何かしらのインサイトを得るのがベイズ推定です。
実際のライブラリの実装
もっと詳しく知りたい方は以下から追えます。
コード内において、変数やattributeがグローバルに使われているため解読が難しいです。オブジェクトが色々な値や構造を内部的に持っているようです。結局、予測値の生成箇所に関しては現在のところ追いきれませんでした。また時間と熱量がある時に追いたいと思います。
・compile_forward_sampling_functionがsampler_fnを作成している。
↓
・sampler_fnがサンプリング実行し、予測値も生成している。
↓
・compile_pymcの実装の中身を知る必要がありそう。
ライブラリの実装の参照先
sample_prior_predictive：事前予測分布の実装。
compile_forward_sampling_function：サンプリング関数の実体。
sample_posterior_predictive：事後予測分布の実装。sampler_fnの使われ方の参考に出来ます。事前予測分布のメソッドでは、sampler_fnの引数がありませんが、事後予測分布の実装には事後分布のサンプル列が必要なため引数があります。
参考：事前予測分布チェックのチュートリアル
事後分布のプロット
こちらは、ベイズを扱う方に対しては釈迦に説法かもしれませんね。
あえて説明するとすれば、左の図は、右の図を90度右回転させて、横軸上にプロットしたものを確率密度関数として表したもの（実際にはカーネル密度推定）になります。MCMCサンプリングの軌跡です。上手く説明できているかわかりませんが…
各チェーンも似たような結果を示せているので、問題なさそうです。
エネルギー診断プロット
az.plot_energy(idata);
<出力>
チュートリアルには、かなりざっくりとした説明が少し書いてあるのみでよくわかりません。
エネルギー診断プロットについて調べましたが、こちらに関してはかなりしっかりと学ぶ必要がありそうです。私も現在のところ数式レベルでの理解は出来ていません。
標準的なベイズ統計の教科書では、MCMCのアルゴリズムの理解として主にギブスサンプリングとメトロポリスヘイスティング法が採用され解説されていると思います。しかし、このエネルギーの概念は、 HMC（ハミルトニアンモンテカルロ法） とその改良版である NUTS (No-U-Turn Sampler） などのハミルトニアン力学を利用したサンプリング手法にしか登場しない概念です。そもそもランダムウォーク型のパラメータ提案手法には「ハミルトニアン」や「エネルギー保存」という概念がないのです。そのため、標準的な教材で勉強しただけではこちらのエネルギーという概念を理解することは難しいかと思います。
詳細な理解には及びませんが、調べた限り分かったことと要点だと思った箇所を下記に記します。イメージのレベルでもわかっていた方が良いと思います。
-
所望のサンプリングが得やすい計算を出来るようにするため、系（パラメータ空間）に対して 運動量$\mathbf{p}$ という便宜的な変数を導入し、擬似力学系を構成する。物理の力学になぞらえて効率よくサンプリングできる枠組みを作ります。
→ 便宜的な変数というと、「ラプラス変換の虚数」や「モーメント母関数の数列」が思い浮かびます。似たような発想で次元（空間）を恣意的に拡張する感じだと思います。 -
パラメータ$\boldsymbol{\theta}$ の勾配を利用する。従来のMCMCでは乱数発生によるランダムウォーク探索であった。HMCでは、ハミルトニアンが出来るだけ保存されるようハミルトニアンの運動方程式を解いたものを利用し、提案を受理すると、無駄な棄却を減らすことができて効率的になります。 ランダムウォークを可能な限り回避しよう とする試みです。
→ 勾配法的に勾配を利用しようというのは自然な発想だと思います。 -
ハミルトニアンの等エネルギー面上を運動させるようにサンプルを生成する。更新する次の値がどこまで飛ぶのかは運動量に基づいて決定される。運動量は位置エネルギーとトレードオフの関係にあるので、位置エネルギーによって決定される。従って、高次元空間上の超平面の高さが低い場所では、運動量を大きく持っているはずなので、遠くへ飛ぶことが出来るため効率は良くなる。勾配に関しては、勾配が方向を決めるし、また勾配が急なところでは加速する。
しかし、 最適化と違い、運動量のランダム生成と受理率の確率的な振る舞いによって、分布全体の探索を行える ようにはなっている。 -
受理率$α$はハミルトニアンの値ではなく、ハミルトニアンの差分（変化量）によって、計算される。
解析力学や統計力学から着想を得た考え方ですので、物理系出身者（主に機電系）だと理解はしやすいのかなと思います。余談ですが、Wikipediaに2次元確率分布のHMCシミュレーションのgifが載っているので、一度見てみるのも面白いです。
（参考）
ハミルトニアン
\begin{aligned}
&H(\boldsymbol{\theta}, \mathbf{p}) = U(\boldsymbol{\theta}) + K(\mathbf{p}),\\\\
&U(\boldsymbol{\theta}) = -\log P(\boldsymbol{\theta}),\\
&K(\mathbf{p}) = \frac{1}{2} \mathbf{p}^\top M^{-1} \mathbf{p}
\end{aligned}
ハミルトンの運動方程式
\begin{aligned}
&\frac{\mathrm{d}\boldsymbol{\theta}}{\mathrm{d}t} = \frac{\partial H}{\partial \mathbf{p}}\\
&\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t} = -\frac{\partial H}{\partial \boldsymbol{\theta}}
\end{aligned}
受理率
\begin{aligned}
&\alpha = \min\left(1, \exp(-\Delta H)\right)\\
&\Delta H = H(\boldsymbol{\theta}^*, \mathbf{p}^*) - H(\boldsymbol{\theta}^{(t)}, \mathbf{p}^{(t)})
\end{aligned}
雑談
確率統計において一般に言えることかもしれませんが、全く違う関数なのに関数を表す文字が同じであるのが、ややこしく混乱しがちになります。
入力変数による関数の区別、条件付きかどうかによる区別。これらは全く違う関数になりますが、同じ文字を使って別の関数を表して理論を説明しているのを見かけます。尤度関数、事前分布、事後分布、事前予測分布、事後予測分布を $p()$ で統一された日には堪ったものではありません(笑)。確率を表す $p()$ は便利な表現ですが、ある種全て確率でもあるのでややこしいです。
出来れば、以下のように統一するのが良いと密かに思っています。（結局は慣れの問題でしょうか...）
- 関数が違う時点で別の文字を使う
- 事前分布と事後分布は$π()$を使う
- モデル（尤度関数）は全て条件付きの表現をする（予測分布を$p(y)$で表すものとかを見かける）
どれがどれだかわからなくなる、様々な表現
\begin{aligned}
&p(x), \quad p(x \mid \boldsymbol{\theta}), \quad p(\boldsymbol{\theta} \mid x), \quad p(y), \quad p(y \mid D), \quad p(D \mid \boldsymbol{\theta})
\end{aligned}
提案したい表現（事前、尤度、事前予測、事後、事後予測）
\begin{aligned}
\pi(\boldsymbol{\theta}),\quad
f(x \mid \boldsymbol{\theta}),\quad
f(x),\quad
\pi(\boldsymbol{\theta} \mid x),\quad
f(x_{\text{new}} \mid x_{\text{old}})
\end{aligned}
エネルギー診断プロットについては、調べるのにここまで時間をかけることになるとは思いませんでした（良い意味です）。
調べていると思わず気になって深入りしてしまいましたが、まだ腑に落ちるほどの理解には及んでいません。そのうちまたやりましょうかね。
最後に、間違いなどお気づきになりましたら、ご連絡いただけますと幸いです。