雑記ブログ【HELLO REALWORLD】始めました。

変分ベイズ法の心２

はじめに

下記記事の続きで、お気持ちは理解している前提で進みます。

www.hellocybernetics.tech

変分ベイズ法の戦略

基本の復習

データ $D = {x_1, \cdots, x_N}$ が手元にあるときに確率モデル $p(x|\theta)$ と事前分布 $p(\theta)$ を設計するのがモデリングの第一歩でした。するとベイズの定理（あるいは乗法定理）から、下記の事後分布を獲得することがベイズ推論の目標になります。

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

さて、この分布を推論したときに最も嬉しい結果は $\theta$ を表す真の分布 $\hat p(\theta) $ にピッタリ一致することでしょう。しかし、私達は真の分布を知らないので、どのようにこれに近づけて行くかというのが問題になります。

一つの解決法として、真の分布に収束するような確率遷移核を上手に設計し、その確率遷移核を持つマルコフ連鎖でサンプリングを繰り返し、ヒストグラムを描かせることで分布を得るMCMCがあります。

もう1つが今回紹介する変分ベイズ法です。変分ベイズ法では推論したい分布をパラメトリックな関数として仮定し、適切な評価指標のもと仮定した近似分布の出来栄えを評価しながら、パラメータを少しずつ変えていく手法です（推論が最適化問題に帰着されている）。

すなわち、事後分布が

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

と表されており、こいつの分布の形状を表す関数を $q(\eta)$ と書き表すことにします。$\eta$ を上手く調整して事後分布を近似してやりたいという話です。実際にはこいつは $\theta$ を発生させる分布であるので、 $q(\theta ;\eta)$ と書いたり、 $q(\theta | \eta)$ と書いたり $q(\theta)_{\eta}$ と書いたりする場合もあります。この $q$ は兎にも角にも $\theta$ を発生させる分布であり、私達がこれから調整するのは $\eta$ であるということを明確に意識しておきましょう。

分布の評価指標

近似分布の出来栄えを評価するのに使われるのが、KLダイバージェンスです。天下り的ですが、KLダイバージェンスは確率分布 $q$ と 確率分布 $p$ の離れ具合を評価することが出来ます。そこで下記を最小化することで、 $q(\theta; \eta)$ の $\eta$ をいろいろ調整してみて $p(\theta | D)$ に最も近い分布を見つけてやろうと試みます。

$$ KL [q(\theta; \eta) : p(\theta | D)] = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta $$

今、どういう状況になったのかを振り返ります。確率分布 $p(\theta |D )$ を求めたいというのがベイズ推論でした。そこで、こいつを一先ずでたらめな $q (\theta)$ で表して少しずつ分布の形を変えながら良さげなものを探そうと考えます。ところがこれでは抽象的すぎるし、分布の変え方をどうするのかもよくわからんので $q(\theta ; \eta)$ とパラメトリックな関数にして、 $\eta$ の調整に問題をすり替えました。そして、どんな分布の形が良いのか（どんな $\eta$ が良いのかを）をKLダイバージェンスで図ろうと決めたのです。

KLダイバージェンスは見た目が複雑に見えますが、もはや $\eta$ の関数に過ぎません。

ただし、ここで注意が必要です。私達は今 $q(\theta ; \eta)$ と事後分布の形状を制限してしまいました。パラメトリックに表された関数で表現できる形状以外のものは $\eta$ をどれだけ調整したところで獲得できないことになります。しかし、それでも、表しうる形状の中ではKLダイバージェンスの意味で最も近い分布を得ることが出来ますので一先ず良いとしましょう。

もう1つ注意すべきは、パラメトリックに関数を置いたことで最適化の進み方も決めてしまっている点です。何度も申してきたとおり、分布の形状を知りたいというのが問題であってその形状を調整するために $\eta$ というパラメータを導入したのでした。言わばこれは分布のイジりやすさを向上するために導入したものであって、それ以外の意味は基本的にはありません（仮定した分布によほどの意図がない限り、こいつはただの調整役に過ぎない）。

もしもパラメータの置き方をもっと工夫してみたり、あるいはパラメータを置かずとも近似分布 $q$ を調整する手段があるならばそれでも構わないのです。よほど良い $q$ の選び方をすれば、上手に調整する方法が何かしら見つかるかもしれませんね。

ELBO

引き続いて、KLダイバージェンスを最小化することをもう少し詳しく見ます。

$$ \begin{align} KL [q(\theta; \eta) : p(\theta | D)] & = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta \\ & = \int q(\theta; \eta) \log \frac{q(\theta; \eta) p(D)}{p(\theta, D)} d\theta \\ & = \int q(\theta; \eta) \log p(D)d\theta + \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta, D)} d\theta \\ & = \int q(\theta; \eta) \log p(D)d\theta - \int q(\theta; \eta) \log \frac{p(\theta, D)}{q(\theta; \eta)} d\theta \\ \end{align} $$

ひたすら式変形を実施しました。ここで第一項に関しては $q(\theta; \eta)$ は確率分布であるので $\log p(D)$ と表すことができます。これは対数周辺尤度、あるいはエビデンスと呼ばれる値です。引数がすでに手元にある確率変数の実現値 $D$ であるので、$p(D)$ は具体的に値を計算できる定数にすぎません（ただし $p$ の形を私達は知らないので今は保留…）。

第二項の方がEvidence Lower Bound、略してELBOと呼ばれるものになります。ELBOは名前の通り、エビデンス（対数周辺尤度）の下界を示しており必ず、

$$ \log p(D) \geq \log \frac{p(\theta, D)}{q(\theta; \eta)} d\theta $$

となることが知られています（話は逆で、そうなることがわかったのでELBOと呼ばれている）。このエビデンスとELBOが答式で結ばれる条件が、

$$ KL [q(\theta; \eta) : p(\theta | D)] = 0 $$

となること、すなわちベイズ推論の結果、$q$ が 事後分布 $p$ をピタリと当てられた時になります。なので変分推論を行う際には「ELBOを最大化する」という表現が使われることも多々あります。エビデンスの方は直接は計算できない保留中の定数でしたので、 ${\rm KL} = {\rm evidence} - {\rm ELBO}$ となっており、KLを最小化することと ELBOを最大化することは同じことをしていることになります。

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

$p(D)$ はデータと $p$ が決まれば自動的に定数となります。当然データは手元にあるものを使うのでこの時点で確率変数の確定値となっており変数として扱う必要はありません。言わば単なる正規化定数（ちゃんと確率の和が1になるように割り算するための定数）のように振る舞います。事後分布の形状を決めているのは分子の方です。

分子の形状を上手く近似してやることができれば求めたい事後分布が求まりそうです。その分子の形状に相当するものを $q(\eta)$ という関数で近似することにします。

$$ q(\eta) = p(\theta | D) = \frac{p(D|\theta) $$

変分ベイズ法の具体的手段

関数 $q$ をどのように置くのか

さて、上記までの話で兎にも角にも

$$ KL [q(\theta; \eta) : p(\theta | D)] = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta $$

を最小化しましょうという話になりました。 さて $q(\theta ; \eta)$ はどのように置くと良いのでしょうか。一つは正規分布のような分かりやすい分布を置いて、 $\eta = (\mu, \sigma)$ みたいにしてしまうことです。そうすると、実際の $\theta$ の事後分布が正規分布でなかったとしても、正規分布の中で最も実際のサンプルに近いものを探すことになります。

もう１つは

$$ \q(\theta, \eta) = q(\theta_1, \eta_1)\cdots q(\theta_K, \eta_K) $$

と置いてみることです（ここで求めたい $\theta$ が$K$次元ベクトル、あるいは $K$個のパラメータを持っていたとしています）。これは各々のパラメータが（データの条件付き）独立であることを仮定していることになります。独立性を仮定しているので、各パラメータ間の相関はつかめなくなることに注意してください。仮に $\theta$ が行列ならば、 各成分それぞれを独立とみなしても良いですし、それがやりすぎだと思ったら、各列が独立とみなして各列が多次元正規分布であると仮定するような方法もあるでしょう。

このように複数のパラメータを有する分布を、個々のパラメータの単純な分布の積に分解してしまう方法を平均場近似と言います。

平均場近似は、事後分布を推論する際に「うまい分布の変え方」を上手くひねり出すときに使えます。例えば、$\eta$ を変分パラメータとして分布の調整役を担ってもらう方法を集中的に取りざたしてきましたが、単純な分布の集まりに変えてしまえば、うまい分布の調整の仕方をパラメータに頼らず考案できる可能性もあるのです（分布をイジる主たる手法が $\eta$ の導入であるというだけで、実は必ずしも使わなくていい…）。

個々らへんをゴリゴリに手計算しながら学びたい場合は

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

をおすすめします。

変分ベイズ法の心

ベイズ推論の基本

ベイズモデリングの概要については下記の記事を参考にしてください。

www.hellocybernetics.tech

概要をさらっとなぞると、ベイズ推論の基本的な話としては、観測データ $x$ の真の確率分布 $\hat p(x)$ を知る由もないので、確率モデル $p(x | \theta)$ でモデル化し、更にパラメータ $\theta$ にも事前分布 $p(\theta)$ を仮定します。

$$ p(x, \theta) = p(x | \theta)p(\theta) = p(\theta | x) p(x) $$

という確率分布に対していつでも成り立っている乗法定理から、

$$ p(\theta | x) = \frac{p(x|\theta)p(\theta)}{p(x)} $$

とできます。そこで 各 $i$ に対して $x_i$ が互いに独立でかつ同じ確率分布から生起しているとすれば、$D = \{x_1, \cdots, x_N\}$ というデータが手元にあるとして、事後分布

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

と書くことが出来ます。 この事後分布を使って、下記のベイズ予測分布を計算し

$$ p(x_{new} | D) = \int_{\theta} p(x_{new} | \theta) p(\theta | D)d\theta $$

真の分布を近似してみたり、あるいは予測に使ってみたりします。これがベイズの基本的な流れになるわけです。

ベイズ予測分布を計算するためには、当然事後分布 $p(\theta | D)$ を予め求めておく必要があり、これを求める（推論する）ことを学習とひとまず呼ぶことにします（実は事後分布を周辺化して消してしまうことを考えれば、予測分布を直接求めるようなアルゴリズムを考えても良いわけで、ベイズ推論は実はもっとフレキシブルです）。

変分ベイズ学習

変分法の心

一旦、ベイズ推論における学習のことに集中しましょう。今達成したいことは、$p(\theta | D)$ という関数を求めるということです。この関数は事後分布と呼ばれてます。この関数を決めることによって$\theta$ の実現値の発生源を決めようとしているわけです（$\theta$を求めようとしているのではなく、 $\theta$ がどのようにバラつくのかを知ろうとしている）。$\theta$ の発生源として使える関数をどのように決定すると良いのでしょうか。

関数を決定する上で大事だと思われるのは、この関数がどんなふうになれば「良い（あるいは悪い）」と言えるかの指標です。そのような指標はすでに与えられているとしましょう。指標を $M(p)$ と書くことにします。仮に良さを表す$M(p)$がどういうものかわかっているなら、これが最大化されるように関数 $p$ を調整してやればいいということです（通常はELBOと呼ばれる指標が使われるが、これは後述する）。

兎にも角にも、$\theta$の発生源となる関数 $p$ を指標 $M(p)$ を使って決めてやろうとしているのです。決して $\theta$ という値をどうするかを決めようとしているのではなく、 $p$ という関数そのものをどのようなものにしてやろうかを決めようとしていることに注意してください。ちなみに、 $M(p)$ は 関数$p$ を引数に取って、スカラー値を返します。このスカラー値が大きければ大きいほど「良い $p$」を入力したということになります。

では関数 $p$ を手当たり次第に試しに入れまくってやって、$M(p)$ を評価していき、一番大きな値を返した $p$ を選べば解決じゃないか！という話です。ただし、当然 $p$ の選び方なんて無数にあるわけなので、この世に存在するありとあらゆる関数（今回の場合は確率分布）を評価して見るなんて無理なのです。

そこで、変分法が登場します。変分法とは、 $M(p)$ に 関数$p$ を試しに入れまくる際に、入れ方を関数の形をほんの少しだけ変えて入れてみるという方法を取ることにします。すなわち、とある $p$ で $M(p)$ を評価したあとに、今度はちょっとだけ $p$ を変えた $ p + \delta p$ で $M(p+\delta p)$ を評価するということです。

もしも上記のことを試してみて、

$$ M(p+\delta p) - M(p) = 0 $$

ということが起こった場合にはどういうことが考えられるでしょう。 $p$ をほんの少しだけ変えても評価値は変わらなかったので、この辺りの $p$ が一先ず解の候補になりそうです。仮に評価指標 $M(p)$ が単峰ならば、上記の $p$ を答えとしてしまうことができます。上記の $M(p+\delta p) - M(p) $ を変分とよんでみることにしましょう。この関数を引数に取ってスカラーの値を返すようなもの（今回で言えば評価指標 $M(p)$ )を「汎関数」と言います。そして、上記の汎関数に入れる関数を少しだけ変えたらどうなるんだろうか？という考えに基づいて変分を考え、関数を決めてやろうという試みを変分法と言います。

なんだか微分と似たような話です！ただし、 $p$ は値ではなく関数なのでもうちょっと話は複雑ですが…。

何が複雑なのかというと、微分であれば、$ f(x) $ の微分を考えるために、$x + \Delta x$ というものを考える必要があります。そのときに $x=5$ に対して $\Delta x = 0.00001$ としてちょっとだけ変えていれてやるか…、みたいなことが出来るわけです。では 関数 $p$ をちょっとだけ変えるために $\delta p$ をどのように選べば良いでしょうか？ 　

　

　 端点を固定しつつ、ちょっとだけ関数の形状をずらす

　 　

　 関数の形状のある一部分のみをずらす 　 上の方が自然な方針に見えますか？それとも下のほうが自然な方針に見えますか？実は私達が求めたいと思っている $p$ の形状の本質を知らない限り、どちらが良いとは言えません（$M(p)$ の評価値を良くするような変更の仕方としか言いようが無いから！）。これが微分に対して変分の方が複雑ですという意味です。（上の方が数学的には扱いやすそうな気はしますよね。物理の１分野である解析力学で学ぶ変分法で最初に学ぶのはおそらくこの形になるでしょう）。

変分ベイズ法の戦略

さて、変分法なるものの存在を知ったところで、変分ベイズ法の話に突入します。まず、私達が求めようとしている関数 $p$ はパラメータ $\theta$ の事後分布です。いま評価指標である $M$ はすでに決まっているとします。上記で見たように、基本的に変分をどのように扱えば良いかの方針というのは簡単には決まりません。$p$ は本来どんな形状であるのかは我々は知りませんし、変分としてどのような変え方をして関数を探していけば良さそうなのかも我々は知りません。

それでも、 $M$ を片っ端から評価していって良さげな $p$ を見つけなければなりません。これを実施するために $p = q(\eta)$ と置いてしまうのです。これは一体何をしたのかというと、 $q(\eta)$ は我々にとって計算しやすく、 $\eta$ を少し変えてやれば $q(\eta)$ の形も少し変わるような良い関数を選んだのです。そうすることで、$p$ の考えうる形状と、そして変分としてどのように変更を加えていくかの方針を $\eta$ に託すことにしたのです。

1. $M(p)$ に基づいて一番良い $p(\theta | D)$ を探したい。
2. $p1(\theta | D), p2(\theta | D),...$ と片っ端から試すのは苦労する。
3. $p$ をほんの少し $\delta p$ だけ変えて、$M(p + \delta p) - M(p)$ を評価していくことに。
4. 関数の探す範囲と変え方 $\delta p$ の方針が分からない。
5. 関数 $p(\theta | D) = q(\eta)$ と置いてしまって、$\eta$ を変えることで変え方と探索する範囲を決定

ということです。そうして得られた $q(\eta)$ をあたかも $\theta$ の事後分布として扱って以後予測分布などに使っていくということになります。この方針で優れている点は、最終的には

$$ M(p(\theta | D)) = M(q(\eta)) $$

と置いてしまったので、変分法で $M(p+\delta p) - M(p) $ を考えることが

$$ \begin{align} & M(p+\delta p) - M(p) \\ \leftrightarrow & M(q(\eta)+\delta q(\eta)) - M(q(\eta)) \\ \leftrightarrow & M(q(\eta + \Delta \eta)) - M(q(\eta)) \\ \leftrightarrow & M(\eta + \Delta \eta) - M(\eta) \end{align} $$

のように書き換えられてしまうことです（$q(\eta)$ は私達が計算しやすく仮定した関数なので、$M(q(\eta))$ という入れ後になった関数は簡単に $M(\eta)$ と入れ子を意識しなくても良い関数に出来てしまうのです）。私達が考えているのは、結局$\eta$ を少し変えることだけなのです。これは従来ニューラルネットワークなどでも使ってきた勾配法が使えます。

そうして $M(\eta)$ を評価して求まった $\eta$ を使って、元々仮定した $q(\eta)$ を $\theta$ の事後分布の形状として採用してやればいいということになります。

TensorFlow2.0 Preview版が出ました！

TensorFlow 2.0発表！

ついに動きがありましたね。APIは下記で見ることが出来ます。名前空間がスッキリしていることに気づくはずです。

www.tensorflow.org

v1.12.0からv2.0へコードを書き換えるためのツールも整備されていく模様です。

tensorflow/tensorflow/tools/compatibility at master · tensorflow/tensorflow · GitHub

また、2.0の発表して間もなく、githubにはチュートリアルのリポジトリが出現しました。さすがは注目度が高いですね。

github.com

コード周辺の変更

さて、TensorFlow2.0でどのように書き方が変わったのかというと、以前からお伝えしてきたとおり、EagerをデフォルトとしたKerasAPIを全面に出した書き方が中心となっていくようです。それに伴って、多くの方が今まで使っていたモジュールにあったはずの関数やクラスがKerasに統廃合されている可能性がかなり高くなっています。

例tf.train.GradientDescentOptimizer()などを始め、tf.train クラスには最適化手法等の実装は一切なくなっておりました（checkpoint等の学習を管理する系のモジュールとなったようです）。代わりに最適化手法は tf.keras.optimizers モジュールに統合されていました。また、tf.lossesなども多くの方が利用していたと思われますが、損失関数などは tf.keras.losses に移動されています。

tf.train → tf.keras.optimizers

tf.losses → tf.keras.losses

細かい例にはなりますが、tf.losses.sparce_softmax_cross_entropyと同様の動きをするものは、tf.keras.losses.sparse_categorical_crossentropyになっているなど、命名方法が若干異なったりするので注意が必要です。また、accuracyやrecallなどを求める便利なmetricsは tf.keras.metricsにまとまっています。

例えばロジスティック回帰を行う場合には下記のようなコードを書くことになるでしょう。

x_train_ = tf.convert_to_tensor(train_inputs_of_numpy)
y_train_ = tf.convert_to_tensor(train_labels_of_numpy)

# Logistic regression model
model = tf.keras.layers.Dense(output_size, activation="sigmoid")

# loss function
def loss_fn(model, x, y):
    predict_y = model(x)

    # return shape is (x.shape[0], ) ; each element is each data's loss.
    return tf.keras.losses.binary_crossentropy(y, predict_y)

def accuracy_fn(model, x, y):
    predict_y = model(x)

    # return shape is (x.shape[0], ) ; each element is 1 or 0.(If y[i] == y_pre[i], the i th element is 1). 
    return tf.keras.metrics.binary_accuracy(y, predict_y)

# optimizer
optimizer = tf.keras.optimizers.SGD(learning_rate=learning_rate)


for epoch in range(num_epochs):
    with tf.GradientTape() as tape:
        loss = loss_fn(model, x_train_, y_train_)
    grads = tape.gradient(loss, model.variables)
    
    accuracy = accuracy_fn(model=model, x=x_train_, y=y_train_)
    
    if (epoch+1) % 5 == 0:
        print(
            "loss: {:0.3f},  acc: {}".format(
                tf.reduce_sum(loss).numpy(),   # using TF function or numpy method
                accuracy.numpy().mean()         # both of ok.
            )  
        ) 
    # update prameters using grads
    optimizer.apply_gradients(zip(grads, model.variables))

tf.layers は完全に廃止されており、tf.keras.layersを利用していくことになるようです。また、地味な変更ですが tf.randomモジュールが出来たので、今まで tf.random_normalなどとしてきたものは tf.random.normalと書くことになります。

学習周りに関してはeager executionを触っていた人にとっては、特に違和感のない変更となっており受け入れやすいものになっています。 一方で、ecosystem周り（TF hubやTF probabilityなど）は2.0向けに整備がまだ完全にはなされていません。要は今回廃止されたAPIを内部に含んでいる場合には動かすことが出来ないというわけです。ここらへんの整備と本体のデバッグ含めてしばらくpreview版が続きそうな雰囲気です。

tutorial

Eagerで書いていたコードを、ブランチ切ってmaster側を2.0に変えていくことにしました。

github.com

また冒頭で述べた下記のリポジトリも整備が進みそうなので注目しておいて良いと思います。

github.com

なにはともわれ、まだ実用段階ではなさそうなので、気長に見守りましょう。

ベイズモデリング勉強の外観

はじめに

今回は下記のツイートが割と評判が合ったので、少し補足のための記事を書きたいと思います。

ベイズやるなら

●モデリング、プログラミング
アヒル本

●推論、計算
須山ベイズ

●理論
渡辺ベイズ

3冊で良いです。

— HELLO CYBERNETICS (@ML_deep) 2019年1月6日

なぜかすごく「いいね」がついているので補足

・アヒル本
数学的難易度低め、データに対するモデリングに注力

・須山ベイズ
数学的難易度は中、腕力がいる、モデリング→推論の流れの数式操作に注力

・渡辺ベイズ
数学的難易度は高、予測分布等の数学的性質（汎化性能との関係）などの理論に注力 https://t.co/28No42C12R

— HELLO CYBERNETICS (@ML_deep) 2019年1月7日

上記発言の意図

まずツイートの意図としては、「アヒル本」、「須山ベイズ」、「渡辺ベイズ」は指導の方向性が異なっており、それがベイズの全体像を埋めていくための知識を補い合っているということを伝えたかったのです。一応それぞれがどのような立ち位置かと言いますと

アヒル本

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

問題設定からデータの把握、そしてモデリング、結果の解釈までをStanの豊富なコード例と共に進めていく極めて実践的な本です。数式は基本的にモデリングをわかりやすく端的に表すためのツールとして用いられており、それほど複雑怪奇な式が出てくることはありません。基本的な確率統計の言葉と、線形代数の計算ができれば読み進めることができる本になっています。

一方で、データ解析を学ぶ本であり、データ解析手法としてベイズモデリングを（好んで）使っているというイメージが合うような気がします。ベイズ統計や周辺のアルゴリズムについて懇切丁寧に解説をしているわけではないということには注意が必要です。

須山ベイズ

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

様々な確率モデルについて「モデリング→推論」の流れを数式で追っていきます。紙面の半分は日本語で残りは計算式だと思ってください（言い過ぎ）。結構腕力がいると思われます。この本を追える、あるいはこの本のように確率モデルを作れるということは、（何らかの計算に向いた言語を知っていれば）そのままプログラムがかけるというわけです（著者の須山氏はJulia langで実装を公開しています）。

問題を解いていくことで理解が深まり身についていくあの感覚そのままです。本書は問題形式にはなっていませんが、答えを見ながらでもちゃんと式をフォローしていくことで勉強になる感覚を覚えているはずです。割とベイズで頻出する式変形を体に慣れさせる感覚の本です。

渡辺ベイズ

ベイズ統計の理論と方法

こちらはベイズモデリングの個々の手法やデータ分析にまつわる具体的な話を扱っている本ではありません。ベイズで扱われる予測分布等の数学的な性質について言及している本であり、理論の話がメインです（一応「方法」とタイトルに書いてありますが、アヒル本のような具体的なデータ解析を想像すると痛い目みます）。

難易度の高い数学が使われており、この本をスラスラ読めるのであればそもそもこんなブログで参考にしないでしょう。 ただスラスラとは読めなくとも、本で紹介されている結論などは知っておいて損はないでしょう。頻度主義とかベイズ主義とか、そんな区別するものでもないことが（結論を認めればですが）わかるはずです。

続いて、ベイズモデリングに関する概観をザッと見ていきます。

確率モデリング

確率モデリングとはデータ $x$ が確率分布 $\hat p(x)$ から生起していると考え、その確率分布 $\hat p(x)$ を手持ちの有限データから探ろうという試みです。そうして作られた確率モデルは新しいデータ$x_{new}$ の予測に使われたり、あるいは確率モデルの形状（あるいはパラメータ）から何らかの解釈に落とし込むことに用いられます。

確率モデリングの概要

通常は手元に有限個のデータ $D = \{x_1, \cdots, x_N\}$ を集めておいて、これらに対する物理的な知見や、ヒストグラム等の可視化から定性的（あるいは定量的）な情報を集めることから始まります。

その後、$x$ が従う確率分布 $\hat p(x)$ は一体どんなものであろうかと思いを馳せながら、既知の確率分布を組み合わせてモデルを構築します。例えば、データはある値を中心に大小に均等にぶれており、中心に近いほどデータの数も多いようなケースは正規分布を使ってみようなどと考えるわけです。

ところがデータが実際に従っている真の分布$\hat p(x)$ が勝手に我々が考えたモデル（例えば正規分布）で表せる保証などありません。単に、そんな気がするのでそう仮定するというだけなのです。そして、そう仮定した中で一番真っ当な形はどんなものであるかを決めることしか出来ません。この時点で、真の分布 $\hat p(x)$ を知ることを我々は諦めていることに注意してください。これから得ようとしている確率モデルは、真の分布と完全に一致していないという意味で間違っています。

それでも「間違っている度合いがどれだけマシか」を測る方法はいろいろあります。実直に「真の分布に比べてどれくらい隔たりがあるのか」ということを近似的に見ることもできます。あるいは「データの予測に対する誤差がどの程度で済むのか」ということを見てやることもできます。これは目的に応じてでしょう。いずれもどれくらいマシなのかを見ることは出来ても、真の分布それ自体を知ることは出来ないことに注意してください。

多くの場合は構築されたモデルから解釈を行ったり、あるいは新しいデータの予測を行ったりすることに使われるはずなので、その目的を達成することが確率モデリングの役割です（決して、集めたデータの真実をピタリと知れるなどと思わないように…）。

確率モデリング手順

ここでデータ $D = \{x_1, \cdots, x_N\}$ を持っているとしましょう。これに対して何らかの考察を行って、確率分布 $p(x | \theta)$ を仮定します。ここで $\theta$ は分布の形状を決定するパラメータです（例えば正規分布ならば平均と分散）。これを確率モデルと呼びます。確率モデルのパラメータ $\theta$ 自体もばらつきを持っており確率変数だと考えます。

例えば「夏場の気温 $x$」に対して正規分布$\mathcal N(x | \mu, \sigma)$を仮定し、「平均気温 $\mu$」なるパラメータを考えることが出来ますが、平均気温自体も年や国によってばらついているはずです（そもそも、いつ、どこの夏場の気温が知りたいのかによりますが）。なのでそのばらつきを反映するような確率分布 $p(\mu)$ を仮定してやることもできるわけです。

話を一般的なところに戻すと、確率モデル $p(x | \theta)$ に対して $\theta$ の確率分布 $p(\theta)$ を更に考えてやることができるというわけです。これを事前分布と呼びます。確率モデル$p(x |\theta)$ と 事前分布 $p(\theta)$ を合わせてベイズ確率モデル $p(x, \theta) = (x |\theta)p(\theta)$ を作ることが出来ます（式自体は条件付き分布の決まり事である）。これは要するにデータとパラメータの同時分布です。ベイズ確率モデリングというのは言わば、データとパラメータの同時分布を考えてやることに等しいのです（言い過ぎか？）。

ところで同時分布と条件付き分布の関係式は

$$ p(x, \theta) = p(x |\theta)p(\theta) = p(\theta | x)p(x) $$

と最右辺の形にできることも注意しておきましょう（余談ですが、同時分布自体はそれぞれの確率変数を完全に対等に扱っており、グラフィカルモデルでは無向グラフになります。一方、条件付き分布を使った表現は有向グラフになっており、ある同時分布を表現する条件付き分布での形式が複数あるというのは、無向グラフを有向グラフに変換する方法は1通りではないということです）。

この関係式をジッと眺めつつ、今、手元にある確率変数の実現値が$ D = \{x_1, \cdots, x_N\}$ であることを思い浮かべれば、

$$ \begin{align} p(\theta | x_1) &= \frac{p(x_1 |\theta)p(\theta)}{p(x_1)}\\ p(\theta | x_2) &= \frac{p(x_2 |\theta)p(\theta)}{p(x_2)}\\ & \vdots \\ p(\theta | x_N) &= \frac{p(x_N |\theta)p(\theta)}{p(x_N)} \end{align} $$

という条件付き分布たちを考えたくなります（条件付き分布というのは $ | $ の右側が確定した場合の 左側の確率変数の確率分布であった）。 今、 $x_i$ が全て独立に同じ分布から生成されてきたものであるとすれば

$$ \begin{align} \displaystyle p(\theta | x_1, \cdots, x_N) &= \frac{p(x_1, \cdots, x_N | \theta)}{p(x_1, \cdots, x_N)} \\ \displaystyle p(\theta | D) & = \frac{\prod_i^N p(x_i | \theta)}{p(D)} \end{align} $$

と表現できます。これを通常は事後分布と呼びます。基本的に、確率モデル $p(x | \theta)$ と 事前分布 $p(\theta)$ を決めた時点でここまではオートマチックに出てくることなります（階層モデルにしたり、複雑にパラメータ・潜在変数が関係し合うモデルでも数式変形あるいはグラフィカルモデルの操作を機械的に行うことで求めたい事後分布を導出することが出来ます）。

事後分布をデータから求めることが推論あるいは学習と呼ばれるフェイズです。

予測モデル

予測モデルの構築方法はいくつか考えられます。まず事後分布 $p(\theta | D)$ は 確率モデル $p(x|\theta)$ と事前分布 $p(\theta)$ を使う上で、$\theta$ がどんな値だと真っ当だと言えるかを確率で表現しているものです。

MAP推定値

そこで事後分布が最大となっている $\theta$ は一番真っ当な $\theta$ であると考えられます。このような $ \theta_{MAP} = {\rm argmax}_{\theta} (p(\theta | D)) $ を決めることをMAP推定と言います（最初からMAP推定すると決めていれば、学習のフェイズで事後分布の推論をせずに最大化問題で楽ができる）。そうして予測時には確率モデルにMAP推定値を代入して $p(x | \theta_{MAP}) $ を予測モデルとして利用することができます。

EAP推定値

次に、事後分布 $p(\theta | D)$ で $\theta$ の期待値 $\theta_{expected} = \int_{\theta} \theta p(\theta |D)d \theta$ を計算して $p(x|\theta_{expected})$ を予測モデルとして使う方法もあります。これはMAP推定に似ているようですが、事後分布が歪んでいたりすると大きく結果が異なって来ることに注意してください（どちらが良いとは言えませんが）。

ベイズ予測分布

最後にベイズ予測分布を紹介します。上記2つは、事後分布からある1つの $\theta$ を計算して使っていました。しかし、ベイズ予測分布は考えうる全ての $\theta$ を考慮してモデルを構築します。ベイズ予測分布は下記の形をしています。

$$ p(x | D) = \int_\theta p(x|\theta)p(\theta|D)d\theta $$

これは確率モデル $p(x|\theta)$ にどんな $\theta $を使えば良いか迷うなら、 $p(\theta | D)$ の重みをつけて色々な $\theta $を混ぜ込んでしまえという心があります。例えば、事後分布が単峰で裾も薄いのであればMAP推定値を使ってしまってもいいかもしれません。しかし仮に事後分布の裾が厚かったり、多峰であったりしたならば1つの値しか使わないなんて明らかに色々な可能性を捨ててもったいないと思うはずです。可能な限りはこのベイズ予測分布を計算してみることをおすすめします。

そして「ベイズ統計」は、このベイズ予測分布が元々考えてたデータ $x$ の真の分布 $\hat p(x)$ を（一致はしないものの上手に）表してくれるであろうと考えて理論が構築されています。したがって、ベイズ理論の枠組みで真の分布との隔たりや、汎化性能を評価しようとする場合にはベイズ予測分布を求めるのが基本ということになります。

ベイズモデリングのまとめ

1. データ $D$ の概要を掴む
2. 確率モデル $p(x | \theta)$ と事前分布 $p(\theta)$ でモデリングする
3. 何らかの推論アルゴリズムで事後分布 $p(\theta | D)$ を得る
4. ベイズ予測分布 $p(x | D)$ を得る
5. モデリングそのものが良かったのかを何らかの方法で評価する

というのが理想的な流れです。

アヒル本は事後分布の計算結果と予測分布の計算結果も掲載されていますが、特に「1.~2.」とその後の結果の解釈「5.」に注力しております。 須山ベイズは「2.~4.」の数式パートを詳細に取り扱っています（github上にはjulia langのコードもあります）。 渡辺ベイズは「2.」と「4.~5.」に関して理論的な立場から言及しており、特に事後分布やベイズ予測分布などベイズならではの推論が真の分布とどのように違うのか、新しいデータに対する予測としてどのような性能を持ちうるのか、そしてそれを評価するための理論的な方法について言及しています。

もちろん、推論アルゴリズムの主力となっているMCMCについて詳しく知りたい場合は、そのための本が別途必要になるでしょうし、確率分布の基本から学びたければそのための本が必要でしょう。また、時系列データであるとか、それに伴い確率過程を扱うような場合には更に高度な数学が（ベイズ推論するかは別として）要求されたりする場合もあります。

上記３冊を読めば「完璧！」と言えるかは分かりませんが、それでもベイズの大事な部分を網羅しつつ、ど真ん中の最重要なところを３冊で重ね塗りできるようなイメージなのでおすすめいたしました。

Optunaでハイパーパラメータチューニング

【簡易速度比較】TensorFlow vs PyTorch

【TF Advent Calendar 2018】TensorFlow Eager Executionの使い方 step by step

TensorFlow Eager Execution Tutorials

TensorFlow Eager_Execution Tutorials 始めました。

PyTorchのTutorialsの充実具合に影響されて始めました。githubにあるPyTorchのtutorialリポジトリを参考に、TensorFlow Eagerへ焼き直し、あるいは適宜内容を変更し、今後も追加していきます。コードは今のところ全てgoogle colabで書かれています。

英語の読み書きは苦手なので雰囲気だけで…。

github.com

フーリエ変換と線形システム【振動解析への応用】

TensorFlow Probability コードメモ ① 確率分布周りの基本と最尤推定

【主成分分析（PCA）まとめ】分散最大化・確率的主成分分析・ベイズ主成分分析（MAP推定）まで

【確率的プログラミング】Edward2, Pyro, PyStanのベイズ線形回帰コードメモ

TensorFlow Eager Execution + Keras API の基本