グラフィカルモデルを使いこなす！～有向分離の導入と教師あり学習～

さて、前回はグラフィカルモデルの描き方と簡単な事後確率の推論をやってみました。今回以降は、下記のようなもう少し現実的な確率モデルに対する推論をグラフィカルモデル上でやってみる予定です。

・教師あり学習（今回）

・半教師あり学習

・共変量シフト

・転移学習

・潜在変数モデル（EMアルゴリズム）

今回は導入として、有向分離（D-separation）と呼ばれる、より複雑なモデルに対する確率変数の独立性をチェックするための手法を紹介します。これを使って、教師あり学習である回帰モデルや識別モデル（2つともグラフィカルモデル上の区別はないです）に対する推論結果がどうなるかを見てみたいと思います。

今回やることも基本的には前回の３つのノードを使った単純なグラフィカルモデルと同じです。

machine-learning.hatenablog.com

あるグラフィカルモデルが与えられ、さらに一部のノードが観測されたとき、残りの観測されていないノードの確率分布（事後分布）を計算します。このとき事後分布は一般的には複数の確率変数が絡みあった複雑な形状をしてしまっています。グラフィカルモデルから変数間の独立性を読み取り、事後分布をよりシンプルな積に分解して表すのが今回の課題になります。

というわけで、その独立性を発見するためのシステマチックな手法である有向分離をさっそく紹介します。

＜有向分離（D-separation）＞

与えられたグラフィカルモデル上でノードＡとノードＢが独立であるか判断したいとします。AとB間の間のすべての経路がブロックされていればＡとＢは独立になります。あるノードCが経路をブロックしているかどうかを判定するには次のフローチャートを使います。

むむ、なんだかでかい図が出てきてちょっと嫌な感じです。慣れないうちは非常にめんどくさそうな手続きに見えてしまいますが、慣れてもめんどくさいので諦めてください。この図は次の具体的なモデルで使っていきましょう。

・具体例）教師あり学習

では、今回は教師あり学習をグラフィカルモデルで表現し、それに有向分離を適用してみましょう。一番単純な教師あり学習は回帰・識別モデルです。これはグラフィカルモデルでは次のように書くことができます。

$x_1$と$y_1$は教師あり学習のための訓練データということにしてください。そして$x_2$と$y_2$はテストデータです。$\theta$は$x$から$y$を推定するための未知のパラメータです。$y$が連続変数を取る場合は回帰で、離散値を取る場合は識別（分類）になります。このグラフに対応する式も書いておきましょう。

\[ p(y_2, x_2, y_1, x_1, \theta) = p(y_2|x_2, \theta)p(x_2)p(y_1|x_1,\theta)p(x_1)p(\theta)\]

グラフと式をよく見比べて対応が取れていることを確認してみてください。モデル作りはこのようにまだデータが1つも観測されていない状態からスタートします*1。

実際には、教師あり学習では、$x_1$と$y_1$が訓練データとして観測され、さらに予測したい$y_2$に対する入力値$x_2$が観測されます。描き直してみましょう。

手元に持っている観測データが条件付けられましたね。ついでに式も書いておくと次のようになります。

\[ p(y_2, \theta | x_1, y_1, x_2) \]

さて、この分布を単純化してみることにしましょう。まず始めに、確率の乗法定理を使って事後分布を積の形で書いてあげます。

\[p(y_2, \theta | x_1, y_1, x_2) =p(y_2| \theta, x_1, y_1, x_2)p(\theta | x_1, y_1, x_2)\]

単純に$y$と$\theta$を2つの項に分けて書いてみただけです。まだちょっと式が長いので、ここでいよいよ有向分離を導入して2つの項をそれぞれダイエットしてみることにしましょう。

・$p(y_2| \theta,  x_1, y_1, x_2)$の項をダイエット

さて、この項をグラフィカルモデルで描いてみましょう。この項では$\theta$が条件付けられているので黒丸になります*2。

さて、$y_2$と他の変数たちの独立性を見ていきましょう。まず始めに、$y_2$と$\theta$は独立でしょうか？そんなわけないですよね。なぜなら2つのノードは隣り合っているので、どうあがいても依存してしまいます。同様の理由で$x_2$も隣にいるので依存してしまいます。

では、$y_1$はどうでしょうか？$y_2$と$y_1$の間の経路は、$\theta$を経由する以外にはありません。したがって、$\theta$が2つのノードをブロックしているかどうかを確かめればOKです。ということで、先ほどのフローチャートを使ってみましょう。

1. $\theta$は●？  Yes!
2. $\theta$は→●←？ No!   =>   ブロックする（B1）

はい、どうでしょうか。ちゃんとB1の結果にたどり着けたでしょうか？$\theta$が$y_2$と$y_1$を結ぶための唯一の経路をブロックしてしまっているので、２つの変数は独立であることが分かります。同じ理由で$y_2$と$x_1$の経路も$\theta$によってブロックされますね。

というわけで、$y_2$に対して、$\theta$と$x_2$には依存関係があり、$y_1$と$x_1$に対しては独立であることがわかりました。これで式が次のようにダイエットできます。

\[ p(y_2| \theta,  x_1, y_1, x_2)=p(y_2| \theta, x_2) \]

・$p(\theta | x_1, y_1, x_2)$の項をダイエット

さて、この項も改めてグラフィカルモデルを描いてみましょう。

この項では$y_2$が登場していないので点線で描いてみました。

さて、ここから$\theta$と他の確率変数との依存性を見ていきます。まず見た瞬間すぐにわかるのは、$x_2$とは独立であるということです。経路そのものが存在しませんよね。

一方で$y_1$は隣り合っているので依存してしまいますね。

では$x_1$はどうでしょうか？先ほどと同様、$\theta$と$x_1$の間の経路を見てみると、唯一、$y_1$が間にいます。これが経路をブロックしているかどうかをフローチャートを使って調べればOKですね。

1. $y_1$は●？  Yes!
2. $y_1$は→●←？ Yes!   =>   ブロックしない（UB1）

おっと、さっきと違って経路上のノードがブロックしませんね*3。この場合は、$\theta$と$x_1$は依存することがわかりました。したがってこの項では$x_2$のみが消えて

\[ p(\theta | x_1, y_1, x_2) =p(\theta | x_1, y_1)\]

となることがわかります。

さて、ちょっと長かったですが、以上、教師あり学習のモデルに対する事後分布を考えると次のような形になることがわかりました。

\[p(y_2, \theta | x_1, y_1, x_2)=p(y_2| \theta, x_2)p(\theta | x_1, y_1)\]

この式の意味するところを解釈してみましょう。

まず右側の$\theta$の分布を見てください。この分布は、パラメータ$\theta$の分布を学習するには訓練データである$x_1$と$y_1$だけ必要だよ、と言っています。テスト入力$x_2$は学習には影響しないようですね。次に左側の$y_2$の分布を見てみると、未観測である$y_2$の分布を推定するためにはパラメータ$\theta$と入力値$x_2$のみが必要で、過去の学習データである$x_1$は$y_1$がどうだったかなんて知らないよ、と言っています。これは教師あり学習のひとつのシンプルさであり、また制限でもあると言えます*4。

ちなみに、実際の応用の場面では学習後の$\theta$なんかどうでも良くて$y_2$の予測だけが知りたいということが多いと思います。この場合は次のように確率の加法定理を使って$\theta$を消してあげる操作が必要になります。これは周辺化（marginalization）と呼ばれています。

\[p(y_2 | x_1, y_1, x_2)　= \int p(y_2| \theta, x_2) p(\theta | x_1, y_1) d\theta \]

この式が簡単に計算出来るかどうかは具体的な確率分布（ガウス分布など）の設定の仕方に依ります。*5

さて、ちょっと長くなってしまったので「めんどくさ」って思われている方もいるかもしれません。しかし、今回の教師あり学習の例をよく振り返って見ると、ベイズ学習で行っていることは確率モデルを設計してその事後確率を推定しただけです。ベイズ推定が最尤推定を発展させたものとして説明されているのをよく教科書とかで見かけますが、そうではなく、単に確率の加法定理と乗法定理を使って条件付き確率分布を求めているだけと捉えると、ベイズの考え方のシンプルさがわかっていただけるかと思います。

次回以降はもっと複雑なモデルに対して有向分離を適用し、事後確率を求めていきたいと思います。

*1:グラフ上ではまるで訓練とテストでデータが１つずつしかないように見えますが、今回の例ではデータが$N$個ある場合でも同じ議論になるので省略しました。

*2:黒丸が「観測データ」だと思い込んでいるとちょっと混乱するかもしれません。この項は$\theta$がある値で条件付けられた場合の$y_2$の確率分布を表しています。

*3:ブロックしないということは，「２つのノードの間の独立性をグラフからでは示せない」ということを意味しています．グラフ上では独立でないように見えても，具体的な数式によるモデル化次第では独立にもなり得ることに注意してください．一方，経路がすべてブロックされる場合は独立性は成り立ちます．

*4:発展的話題ですが、テストデータの入力$x_2$もパラメータの学習に含めたい場合は、共変量シフトと呼ばれる入力$x$の分布を考慮する学習モデル等を考える必要があります。また学習データそのものによって柔軟に未知変数の分布を推定したい場合には、ガウス過程等のカーネルモデルを使う必要があります。

*5:この例だと、例えば$p(y|x)$に対してガウス分布を設定し、$p(\theta)$に対して共役事前分布であるガウス・ウィシャート分布などを設定すれば解析的に計算できます。また別の機会で具体的な計算方法を説明したいと思います。

グラフィカルモデルによる確率モデル設計の基本

今回から数回にわたって、グラフィカルモデルを利用した確率モデルの設計についてお話しします。従来の統計モデルと比べ、機械学習を機械学習たらしめているものの一つは、扱う現象の複雑さにあると言えます。複雑な現象を解析するためにはそれに見合った複雑なモデルが必要で、それを簡潔に記述するための方法としてグラフィカルモデルが開発されました。

「グラフィカルモデルを使って現象をモデル化し、必要に応じて近似推論法を用いて未知の値を推定する」

という一連の流れが身につくと、いろんなデータサイエンスの課題に対してシンプルかつフォーマルに取り組めるようになります。

それではまず始めに、超超超重要な確率の加法定理と乗法定理の確認をしてみましょう。 

・加法定理（sum rule）*1 

\[ p(x) = \sum_y p(x,y) \]

・乗法定理（product rule）

\[ p(x,y) = p(x|y)p(y) \]

$p(x,y)$は同時分布（joint distribution）と呼ばれ、対称性$p(x,y) = p(y,x)$が成り立ちます。$p(x|y)$は条件付き分布（conditional distribution）と呼ばれ、$y$が与えられた時の$x$の分布です（$y$の分布ではありません）。

ベイズ学習と呼ばれる機械学習の技術は、実は上の2つのルールをただ黙々と使っているに過ぎません。残りは、これらを具体的に計算するための部品である確率分布たちと、（必要に応じて使う）近似推論くらいしかありません。

また2つのルールから有名なベイズの定理が導けますので、重要な結果として覚えておいて損はないです。

・ベイズの定理（Bayes' theorem） 

\[ p(x|y) = \frac{p(y|x)p(x)}{\sum_x p(x,y)} \]

さらに用語として知っておいてほしいのが独立です。次の式が成り立つときに限り、$p(x)$と$p(y)$は独立であると言います。

\[ p(x,y) = p(x)p(y)\]

ちなみに両辺を$p(y)$で割ってあげると、分布の独立性は次のようにも書けます。

\[ p(x|y) = p(x)\]

これも覚えておいて損はないです。

さて、ここからが今日の本題なのですが、今回はDAG（Directed Acyclic Graph）と呼ばれる、ループ構造のないもっともポピュラーなグラフィカルモデルを扱うことにします*2。まず上記の表記を使った確率モデルと、ノードと矢印を使ったDAGとがどのように対応付けられるのかを確認したいと思います。こういうのは細かい定義をうだうだ言う前に例示した方が早いかと思います。 

モデル１） head to tail

3つの白丸ノードA,B,Cが右向きの矢印のみでつながっています*3。白丸が3つなので確率モデルの式も次のように3つの分布の積で書けます。

\[ p(A,B,C)=p(A)p(B|A)p(C|B) \tag{1} \]

ノードそれぞれに対して確率分布を$p(\cdot)$を置き、矢印の元になっているノードを「条件」として右側に書けばいいんですね。

さて、もしこのモデルでノードBが「観測」された場合はどうなるのでしょうか。「観測」とは、確率変数に具体的な数値が与えられるという意味です（条件付けられる、の方がより正確です）。グラフィカルモデルでは観測されたノードは黒丸で表現されます。

はい、こんな感じです。このとき、AやCの確率分布はどうなるのでしょうか。Bは黒丸なのでもはや分布を持たず、白丸だけの分布を考えればOKですね。

\[p(A,C|B)\]

ということで、これを式(1)と、確率の加法定理と乗法定理を使ってちょこっと計算してみましょう。

\[p(A,C|B) = \frac{p(A,B,C)}{p(B)} = \frac{p(A)p(B|A)p(C|B)}{p(B)} = p(A|B)(C|B) \]

というわけで、ノードAとCは、Bが観測された状態では「独立になる」ことがわかりました。これは条件付き独立って呼ばれています。あるいは「BはAとCをブロックする」という言い方もします。

モデル１におけるこの結果はとても重要です。一度Bが観測されてしまえばAとCの値は相関を持たなくなるんですね。

モデル2） tail to tail

さて、AB間の矢印の向きが先ほどと違い、今度はノードBがA,Cの親になっているような状態です。式で書くとこんな感じです。

\[ p(A,B,C) = p(B)p(A|B)p(C|B) \tag{2} \]

先ほどと同様、Bを観測してみることにしましょう。

この事後分布が先ほどの例のように条件付き独立になるかどうか計算してみましょう。

\[ p(A,C|B) = \frac{p(A,B,C)}{p(B)} = \frac{p(B)p(A|B)p(C|B)}{p(B)} =p(A|B)p(C|B) \]

はい、というわけで今度もまたAとCはBが与えられたとき独立になるようです。ノードBはAとCをブロックする、と言うんでしたね。

モデル3） head to head

最後ですが、このhead to headと呼ばれる3つ目のモデルが一番のやっかいものです。これを説明したいがために前の2つのモデルを紹介したと言っても過言ではありません。

対応する式は次のように書けます。

\[ p(A,B,C) = p(A)p(C)p(B|A,C) \tag{3} \]

今までのモデルと違い、2つの変数が条件付けられているような項$p(B|A,C)$が出てきましたね。ちなみにこのとき、グラフから読み取っても直感的かと思いますが、AとCは独立です。念のため加法定理を使って調べてみましょう。

\[ p(A,C) = \sum_B p(A,B,C) =\sum_B p(A)p(C)p(B|A,C) = p(A)p(C) \]

$\sum_B p(B|A,C)=1$を使ったのは大丈夫でしょうか。Bの確率分布なので足し合わせれば絶対に1ですね。というわけで、AとCは独立です。これは観測されていないノードBがAとCをブロックしている、とも言えます。

さて、次にBを条件付けしてみましょう。

今までどおり、対応する事後分布が条件付き独立になるか調べてみましょう。

\[ p(A,C|B) = \frac{p(A,B,C)}{p(B)} = \frac{p(A)p(C)p(B|A,C)}{p(B)} \]

これは他の2例と違って一般的に$p(A|B)p(C|B)$に分解することは出来ません。これ以上どんなに式をコネくり回してもダメです。もともとは独立だったAとCが、Bが観測されることによって互いに依存するようになってしまいました。

このように、モデル3ではBを観測することによって、AとCの分布がより複雑なものになってしまいました。実はこの現象が、ベイズ学習において近似推論手法（変分近似*4やMCMC*5）が必要とされる理由なんですね。

さて、今回はグラフィカルモデルの導入と、簡単なモデルにおける事後分布の推論を見てみました。次回は有向分離と呼ばれる、もっともっと複雑なグラフィカルモデルに対するノードの独立性の判定アルゴリズムを紹介します。さらに、具体的ないくつかの例（教師なし学習、回帰・識別、半教師あり学習、共変量シフト、転移学習、など）をグラフィカルモデルで表現し、それらに対して有向分離を適用してみたいと思います。

[続き・関連]

machine-learning.hatenablog.com

今回の記事がよくわからん！という方には，次のような入門書もあります．

books.rakuten.co.jp

*1:$y$が連続値を取る場合は和$\sum_y$の代わりに積分$\int dy$を使えばOKです。

*2:このような確率モデルの表記の仕方はベイジアンネットワークとも呼ばれます。ちなみに考案者である計算機科学者のJudea Pearlは2011年にこの功績によってチューリング賞を受賞しています。

*3:ちなみに矢印の根元にあるノードは先にあるノードに対して親であると言います。逆は子です。この場合は例えばAがBの親であるとか、CがBの子であるとかって言ったりします。祖先や子孫も直感どおり定義されます。

*4:変分近似の基礎 

http://machine-learning.hatenablog.com/entry/2016/01/23/123033

*5:MCMCと変分近似の比較 

http://machine-learning.hatenablog.com/entry/2016/02/04/201945

もうバグに悩まされることもない？MITの研究者が人工知能を利用した自動デバッグシステムを開発

こんにちは。

MITの研究者が従来法よりも10倍多くのバグを修正できるアルゴリズムを開発したそうです。

news.mit.edu

・従来法の問題点

ソースコードのバグを自動修正するという研究は従来からありました。一般的な方法としては、修正したいソースコードに対する正解セット（入力値とそれに対する正しい出力値）をいくつか人間が用意してあげて、その正解セットに合致するような修正提案をシステムが一生懸命探すというものです。このような手法は論理的な計算を多く必要とするので時間がかかるのと、結局修正コードを作ってもらっても与えたセット以外の入力に対してめちゃくちゃな値を出して使えないとかいった問題がありました。

・大量のパッチからプログラムの修正パターンを学習

今回のMITの研究者が開発したアルゴリズムの革新的なところは、ソースコードの修正に機械学習技術*1を利用した点です。研究者らは、GitHubから大量のオープンソースプログラムの修正パッチを取得し、それを「学習データ」としてコードを自動修正するための汎用的な規則を抽出しました。機械学習アルゴリズムを使うには、生データ（今回はソースコード）に対してどのような特徴量を設計したらいいのかがひとつの重要なポイントになります。この研究ではソースコード上の変数にいくつかの特性（変化する値か定数か、グローバルかローカルか、など30種類）を与えることによって特徴量を設計しているようです。この特徴量の空間でバグの入ったコードと修正パッチの規則性を学習します。従来技術では1個か2個ほどのバグしか修正できなかったのに対して、今回の手法は15から18個ものバグを修正できたそうです。

・これからの課題と雑感

今のところはソースコード全体に影響を及ぼすような複雑なバグは修正できず、ローカルな一対一対応のような小さなバグを修正するのみとなっているようです。とはいえ、それでも人間のかなりの作業量を減らせることに間違いはありません。また、アルゴリズムの学習データとなるソースコードは日々オープンに蓄積されていきますし、今回の技術もまだまだ洗練されているとは言えないので、これからの技術改善に注目です。

*1:「人工知能」は人間などの「知能」をコンピュータ上に実現するための技術の総称です。それに対して「機械学習」はもうちょっと工学寄りで、データから規則的な構造を自動発見して未来の予測や判断を行うための技術の総称です。

MCMCと変分近似

今回は代表的な２つの確率分布の近似推定手法であるMCMCと変分近似を比較します。変分近似に関しては複数回にわけて記事にしているのでそちらを参照されるとよいです。

変分近似（Variational Approximation）の基本（１）

変分近似（Variational Approximation）の基本（２）

変分近似（Variational Approximation）の基本（３）

さて、MCMC（Markov Chain Monte Carlo）は、サンプリング手法の一種です。サンプリングでは、解析的に計算できない事後分布の統計量などを、データをサンプリングすることによって近似的に求めます。今回はMCMCの中でも一番シンプルで便利なギブスサンプリング（Gibbs sampling）と呼ばれる手法を紹介します。

前回取り上げた2次元ガウス分布の近似問題をまた例に取り上げます。今回はこの分布を近似推定するギブスサンプリングのアルゴリズムを導き、過去に導出した変分近似による推定アルゴリズムと簡単な比較してみたいと思います。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 前回の変分近似の内容
• 多次元ガウス分布
• KL divergence

ギブスサンプリングでは、推定したい多変数の事後分布に対して、下記のような手続きで個目のサンプル値とを拾ってきます。

他の変数のサンプル値で条件付けされた分布から新たにサンプル取る、というような手続きですね。これ以上の多変数の場合も同様です。事後分布から一度にすべての変数のサンプルを取るのが難しいのであれば1つずつ変数をサンプルしてしまえばいい、という発想です。この方法が真の事後分布からサンプルしていることを証明することができるのですが、今回は割愛します。*1

さて、変分近似の更新則と比較してみると面白いです。 

<変分近似>

- より単純な事後分布の関数形を仮定し、積分を解析的に行えるようにする。

- 注目している確率分布以外の確率分布に関して期待値を取ることにより、を得る。

<ギブスサンプリング>

- より単純な条件付き分布からサンプルし、積分をサンプルで近似する。

- 注目している確率変数以外の確率変数をサンプルされた値で固定し、をサンプルする。

こうして比較してみると、ギブスサンプリングと変分近似が兄弟のような関係になっていることがわかります。あとで見るように、実際にアルゴリズムを導いてみると同じような式が出てきます。

さて、前回と同様、今回ギブスサンプリングで推定したいのは下記のような2次元ガウス分布です。

前にも触れましたが、推定したい分布が2次元ガウス分布だとわかっているので、実はあんまり意味のない問題を解いていることになります。しかしこの単純な例を使うと、真の分布との間の近似誤差を解析的に計算できるので、２つの手法の近似精度を定量的に比較できるというメリットがあります。というわけで、今回は「1次元のガウス分布は簡単にサンプルできるけど、2次元になるとサンプルできない」という架空の想定で進めたいと思います。

では、実際に条件付き分布を求めてみましょう。を固定されたサンプル値であるとして、の確率分布を求めてみます。

これを正規化することにより1次元のガウス分布が得られます。

ただし、

です。

1次元のガウス分布によるサンプリングは「簡単」なんでしたね。変分近似と似ていますが、分布を細かくした結果が十分に簡単にサンプリングできるような形になっている必要があります*2。また実装でこの式からをサンプルするには、お使いのプログラミング言語にたぶんあるrandn()のような1次元ガウス分布からサンプルを取るライブラリーを使えばOKです。 

のサンプルも全く同様なので省略します。

さて、目的は2次元のガウス分布を近似したかったんでした。上の手続きで個のサンプルを得た後、求めたい分布の統計量（平均、分散）は次のように推定することができます。

では実際に動作させてみましょう。 

 

 青い楕円が推定したい真の2次元ガウス分布で、赤い楕円が推定中の分布です。ピンク色の点が実際のサンプルです。最終的に形状（共分散）も含めてきちんと正解の分布に収束していっていることが分かります。これは分解してしまった変数間の相関を補足できない変分近似とは対照的です。（前回の記事参照）

次にギブスサンプリングと前回導いた変分近似の近似能力を比較してみたいと思います。評価は、真の分布（2次元ガウス分布）とのKL divergenceを直接計算することにより定量的に求めたいと思います。繰り返しになりますが、今回は既知の2次元ガウス分布を使っているためこのような比較ができるのですが、実際は真の確率分布に対するKL divergenceは解析的には求まりません。

 本当は2つの手法を比較するためには横軸は実計算時間を取らないとだめなのですが、どっちもこの例では高速すぎるので繰り返し回数にしました。一般的には変分近似の方がギブスサンプリングよりも早く収束するのですが、この例だと分布が十分に簡単なので、ギブスサンプリングの方が少ない繰り返し数であっという間に変分近似（最適値.46）を抜き去っています。

以上、今回の記事からわかる範囲と個人的な見解にもとづいて変分近似とギブスサンプリングの利点欠点をまとめます。

<変分近似>

- 良い点

計算が早い。収束判定がしやすく、バグも発見しやすい（目的関数があるので）。

- 悪い点

導出がちょっと大変。近似の仕方によっては重要な相関が取れなくなる。

<ギブスサンプリング>

- 良い点

理論的には真の事後分布からサンプルしていることになる。導出が比較的簡単。並列化しやすい。

- 悪い点

大規模なモデルでは収束が遅い。またどれくらいサンプルすればいいのか判断がしにくい。

一般的には小規模問題ならギブスサンプリング、大規模問題なら変分近似、という風に解釈もできますが、計算環境や問題によっても変わってくるので出来れば両方試すのが理想です。

発展的な話題として、2つの技術を組み合わせたような手法を作ることも可能です。他にも、崩壊ギブスサンプリング（collapsed Gibbs sampling）という、日本語だと少し中2病こじらせたような名前の手法があり、これに関してもまた別の機会で書きたいと思います。

今回の記事がよくわからん！という方には，次のような入門書もあります．

books.rakuten.co.jp

*1:サンプリングの手続きが不変（invariant）であることと、エルゴディック（ergodic）であることを示す必要があります。ギブスサンプリングは証明済みであるので安心して使えます。

*2:今回の例のように必ずしも正規化できる必要はありません。よくわからない確率分布が出てきても、例えば棄却サンプリングのような方法を使えばうまく点がサンプルできる可能性があります。

変分近似（Variational Approximation）の基本（３）

「作って遊ぶ」を題目として掲げておきながらまだ作っても遊んでもいなかったので、今回はそろそろ何か動くものを載せたいと思います。

さて、前回得られた変分近似のアルゴリズムを導出するための手引きを使って、今回は世界で一番簡単だと思われる２次元ガウス分布に対して近似推定をやってみたいと思います。*1

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 前回の記事の内容
• 多次元ガウス分布

今回は２次元のガウス分布の近似推定を例として行いますが、実を言うと、多次元ガウス分布は積分も解析的にできますしサンプリングも簡単にできるような単純な分布なので、近似分布をわざわざ求める意味は皆無です。しかし、この例は計算がとてもシンプルで変分近似の導出手順を説明しやすいのと、変分近似が「近似してしまっているもの」が何なのか明確化することができるので、基本を説明するには十分な例だと思います。

では、次のような２次元のガウス分布を変分近似を使って近似推定してみましょう。

\[ p(x_1,x_2|\mu_1,\mu_2, \Lambda) = \mathcal{N} \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] \bigg| \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right], \left[ \begin{array}{cc} \Lambda_{1,1} & \Lambda_{1,2} \\ \Lambda_{2,1} & \Lambda_{2,2} \end{array} \right]^{-1} \Bigr) \\ \propto exp\Bigl\{ \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] - \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right] \Bigr)' \left[ \begin{array}{cc} \Lambda_{1,1} & \Lambda_{1,2} \\ \Lambda_{2,1} & \Lambda_{2,2} \end{array} \right] \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] - \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right]\Bigr) \Bigr\} \]

$x_1$,$x_2$はガウス分布からサンプルされる確率変数です。$\mu_1$,$\mu_2$は平均値、$\Lambda$は$x$の精度行列*2で、２つとも今回は適当な値で固定されたパラメータです。

さて、これをある近似分布$q(x_1,x_2)$で推定しようと思います。２変数の確率分布を分解して推定しようとしているので、次のように2つの分布に分解するしか今回は選択がないです。

\[q(x_1, x_2) = q(x_1)q(x_2) \]

さて、前回紹介した公式

\[ \ln q(z_1) = \langle \ln p(z_1, z_2| x) \rangle_{q(z_2)} + c \]

を適用してみましょう。*3

\[ \ln q(x_1) = \langle \ln p(x_1, x_2|\mu_1,\mu_2, \Lambda) \rangle_{q(x_2)} + c \\ = -\frac{1}{2} \langle x_1^{2}\Lambda_{1,1} - 2 x_1 (\Lambda_{1,1} \mu_1 - \Lambda_{1,2}(x_2 - \mu_2)) \rangle + c \\ = -\frac{1}{2} \{x_1^{2}\Lambda_{1,1} -2 x_1 (\Lambda_{1,1} \mu_1 - \Lambda_{1,2}(\langle x_2 \rangle - \mu_2)) \} + c \]

求めたいのは$x_1$に関する確率分布です。なので$x_1$にだけ注目し、無関係な項をすべて定数$c$に吸収させてしまうのが計算上のポイントです。さらにブラケット$\langle \cdot \rangle$を使って表現した期待値計算ですが、ここでは$x_2$のみに適用してあげればOKで、$x_2$に無関係な項たちはブラケットをするりと抜け出すことができます。

さて、この式をよく見てみると、$x_1$の「上に凸の2次関数」になっていることがわかります。対数計算の結果が上に凸の2次関数になっているということは、この確率分布は1次元のガウス分布であることを表しています。したがってこの式から、平均と分散を求めてあげれば近似分布が求まります。*4

\[ q(x_1) = \mathcal{N}(x_1| m_1, \Lambda_{1,1}^{-1}) \]

ただし、

\[ m_1 = \langle x_1 \rangle = \mu_1 - \Lambda_{1,1}^{-1}\Lambda_{1,2}(\langle x_2 \rangle - \mu_2) \]

です。

$q(x_2)$の計算も同様に計算でき、下記のようになります。*5。

\[ q(x_2) = \mathcal{N}(x_2| m_2, \Lambda_{2,2}^{-1}) \]

\[ m_2 = \langle x_2 \rangle = \mu_2 - \Lambda_{2,2}^{-1}\Lambda_{2,1}(\langle x_1 \rangle - \mu_1) \]

以上から、得られる疑似コードは次のようになります。

1. $m_2$をランダムに初期化する。
2. $m_1$を更新する。
3. $m_2$を更新する。
4. 以上、2と3を十分な回数まで繰り返す。

結果的に近似分布の精度は$\Lambda_{1,1}$と$\Lambda_{2,2}$のまま更新されず、平均値だけ更新されていくようなアルゴリズムになりましたね。

さて、これを実装して動かした結果が次のものです。

上図では、青い楕円が推定したい真のガウス分布で、赤い楕円が推定中の近似分布です（σ=1にあたるところで楕円を描いています）。繰り返し回数は50回にしています。ランダムな平均値から出発して、だいたい15回目くらいの更新で収束しているように見えますね。下図では対応する真の分布と近似分布との間のKL divergenceを繰り返しごとにプロットしてみました。

要点を少し挙げてみます。

１、近似分布が共分散を表現できていない

上図の赤い楕円（近似分布）は常に軸に平行になっており、絶対に斜め向きにはなってくれません。これは近似分布の独立性（分解）を仮定しているため、x1とx2の相関が表現できなくなっているためです。このように変分近似では、最初に独立性を仮定してしまった変数間の相関は取ることができません。

２、KL divergenceが単調に減少している

下図を見ればわかるように、真の事後分布と近似分布との間のKL divergenceが単調に減少していることが分かります。これは、毎回の更新でKL divergenceを最小化する方向に近似分布を修正しているからです。もしこれが増加する場合はバグなので、更新式やソースコードを見直してみる必要があります。

今回の実験に関するソースコードは時間があるときにGitHubに上げようと思います。

すっかり忘れていましたが、Juliaで実装したものをGitHubに公開しました。

MLBlog/demo_simpleVI.jl at master · sammy-suyama/MLBlog · GitHub

基本的にはJulia上で

julia> include("demo_simpleVI.jl")

と叩けば今回のような図が色々出てくるかと思います。

次回以降は、もうちょっと複雑な、でも現実的なモデルに対して変分近似を適用してみたいと思います。

[続き・関連]

MCMCと変分近似 - 作って遊ぶ機械学習。

*1:まったく同じ例がPRMLにもあります。世界で一番簡単な例なので許してください。

*2:精度行列は共分散行列の逆行列なのですが、今回の例では共分散行列を使うよりもこの精度行列の方が導出がスッキリします。

*3:前回導いた公式では右辺は観測データによって条件づけされていますが、今回はそれに該当するのはガウス分布のパラメータ（$\mu$、$\Lambda$）です。ベイズ学習における推論では、既知の固定パラメータと観測データには数学的な区別はありません。

*4:高校数学で教わった平方完成を使います。2次式を$-\frac{1}{2}(x - m)' \Lambda (x - m)$と置いて展開してあげれば、逆からどのように平均と精度を求めたらいいかがわかるかと思います。

*5:実は2つの更新式を連立方程式として解くと2つの平均値$\langle x_1 \rangle$と$\langle x_2 \rangle$が解析に求まってしまいます。今回は実際の多くの応用のように、あえて繰り返しアルゴリズムを使って解を求めます。

Googleの人工知能が囲碁のトッププレイヤーを撃破

オセロ、チェス、将棋に続き、ついに囲碁までもコンピュータが人間を上回るようになったみたいですね。

www.wired.com

概略は以下の通りです。 

・GoogleのDeepMindの研究員が開発したAlphaGoという囲碁の人工知能システムが、欧州のチャンピオンを倒した。

・深層学習（ディープラーニング）を使って過去の対局データ（棋譜）を大量に学習し、人間のように局面を「直観的に」判断するような技術を実現した。

・他にも強化学習や従来の探索アルゴリズムを組み合わせたり、AI同士を競わせるなどして戦略を向上させた。

んーなるほど。どれだけ過去に同じアプローチで研究が行われていたかはちょっと調べていないのですが、私はもう数年はかかるんじゃないかと思っていました。

・囲碁の難しさ

囲碁はオセロやチェスと比べると格段に難しい問題です。IBMが開発したディープブルーというシステムがチェスのチャンピオンを倒したのは1997年ですが、当時の技術はしらみつぶし探索に近い手法を利用して最善の手を選択していました。

しかし囲碁ではボードのサイズが19x19と非常に大きく、単純計算をするとたった7手先を読むだけで80京通り（800000000000000000通り）の組み合わせを探索しなければならず、現状ではどんなにパワフルなコンピュータであってもしらみつぶし探索は使えません。

・深層学習（ディープラーニング）を利用

今回囲碁チャンピオンを倒したアルゴリズムに関して特筆すべき点は、パターン認識のアルゴリズムである深層学習を用いているところでしょう。これは極めて自然なアプローチだと思います。囲碁のプロ棋士たちはよく「模様」という言葉を使って盤上の戦況を判断します。これは本当に文字通り「模様」で、特に対局の序盤から中盤では「白石はこのあたりが厚いから（集まっているから）有利そうだ」とか「このあたりに石を打っておけば黒の陣地拡大を抑えられる」といった、大局的な「見た目」から次の一手を決めます。終盤になってくると、可能な手数が限られてくるので、しらみつぶし探索のようなロジカルな読みに次第に移行していきます。大局的な（直観的な）局面の判断と、局所的な（精緻な）読みの組み合わせを要するところが囲碁の難しいところであり、面白いところでもあります。今回のアルゴリズムでは、このような人間の棋士が普通に行っているような「直観的な判断」を、画像認識などで使われているパターン認識の技術を用いてコンピュータに学習させたわけです。

・大量の棋譜データで学習

さて、学習の仕組みがあるだけではアルゴリズムは強くなりません。特筆すべきもう1つのポイントは、充実した過去のプロ棋士の対局データ（棋譜）を利活用したことです。これをアルゴリズムに与えることにより、どのような局面でどのように手を打ったら良いのかを学習することができます。このアプローチは将棋のAIでも同じで、電脳戦でプロ棋士に勝ったコンピュータは大量の棋譜データから指し手の方針を学習しています。

・現時点で「最良」のアプローチではあるが「最高」ではない

ここからは私の個人的な考えです。私はこの深層学習を使ったアプローチが現時点では「最良」ではありますが、囲碁を解くための「最高」のアプローチではないと思っています。強い囲碁AIを作るという課題は、深層学習が得意とする画像や音声などの認識技術とは決定的に違う点があります。それは目標の設定の仕方です。

例えば、画像認識は、写真に写っている物体が何なのかを当てるといった問題ですが、これは人間が正しいラベルを与えているので、人間の認識能力を実現することが最終的なゴールになります。したがって深層学習のような人間の認識プロセスを模倣するようなアプローチがおそらく一番素直かつ最適なアプローチであると言えるでしょう。

しかし囲碁はゴールが異なります。囲碁は厳格なルールなもとに作られた、極めて人工的な問題です。人間の指し手を模倣すること自体がゴールではありません。囲碁では、究極的にはしらみつぶし探索が最強であり、人間は過去の長い歴史の経験をもとに「そこそこ効率的な」探索手法を発見し、それを棋譜や定石にまとめたに過ぎません。極端な話になりますが、例えば量子コンピュータなどが実現して極めて困難な探索問題をしらみつぶしに近い手法で高速に解けるようになったとしたら、今のプロ棋士や今回のテクノロジーとは全く異なる指し筋になるでしょうし、比べ物にならないほど強いものになると思われます。

人間のプレイヤーを倒してしまった以上、今後この研究が進むかどうかはわかりません。しかし、人間もAI技術も囲碁という問題に対してある程度の近似解を得たに過ぎず、我々はまだ囲碁を制していないと言えると思います。

最尤推定、MAP推定、ベイズ推論

今回は、最尤推定、MAP推定（事後確率最大化推定、正則化）、ベイズ推論*1の関係性を見ていきたいと思います。結論から言うと、最尤推定とMAP推定はベイズ推論の特殊な近似方法であると見ることができます。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• KL divergence
• 確率の加法定理、乗法定理

\[ \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \]

$x$を観測データ、$\theta$をパラメータとした確率モデル$p(x, \theta)$に対して、それぞれの推定方法は一般的には下記のように認識されているようです。

・最尤推定

\[ \theta_{ML} = \argmax_{\theta} \{ p(x|\theta) \} \tag{1} \]

・MAP推定（正則化）

\[ \theta_{MAP} = \argmax_{\theta} \{ p(x|\theta)p(\theta) \} \tag{2} \]

・ベイズ推論

\[ p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)} \tag{3} \]

ベイズ推論では基本的にベイズの定理を用いて事後分布を求めているだけであり、argmaxのような関数の「最適化」は含まれていません。*2

まずはベイズ推論に比較的「近い」MAP推定から見ていきましょう。MAP推定をベイズ推論の近似であるとすると、事後分布$p(\theta|x)$に次のようなパラメトリックな分布を仮定していることになります。

\[ p(\theta | x) \approx q(\theta|\hat{\theta}) = \delta(\theta - \hat{\theta}) \]

ここで$\delta(\theta)$はデルタ分布であり、次のような性質を持ちます。

\[ \delta(\theta) = \begin{cases} +\infty & (\theta=0) \\ 0 & (otherwise) \end{cases} \]

\[ \int_{-\infty}^{+\infty} \delta(\theta) d\theta = 1 \]

\[ \int_{-\infty}^{+\infty} \delta(\theta)f(\theta) d\theta = \langle f(\theta) \rangle_{\delta(\theta)} = f(0) \]

簡単に言ってしまうと、$\theta=0$の点で無限大に尖った形状を持つような確率分布です。したがって$\delta(\theta - \hat{\theta})$は$\theta = \hat{\theta}$の時に無限大の値を持つような確率分布です。

さて、この仮定のもとで真の事後分布とのKL divergenceを最小化してみましょう。

\[ KL(q(\theta|\hat{\theta})||p(\theta|x)) \\ = \langle \ln q(\theta|\hat{\theta}) \rangle_{q(\theta|\hat{\theta})} - \langle \ln p(\theta|x) \rangle_{q(\theta|\hat{\theta})} \\ = \ln q(\hat{\theta}|\hat{\theta}) - \ln p(x|\hat{\theta}) - \ln p(\hat{\theta}) \\ = - \{ \ln p(x|\hat{\theta}) + \ln p(\hat{\theta}) \} + c \]

$\ln q(\hat{\theta}|\hat{\theta})$は$\hat{\theta}$の値にかかわらず無限大の値を取るため、 定数項$c$に吸収させました。さて、この式を$\hat{\theta}$に関して最小化するので、

\[ \theta_{MAP} = \argmin_{\hat{\theta}} \{ KL(q(\theta|\hat{\theta})||p(\theta|x)) \} \\ = \argmax_{\theta} \{ \ln p(x|\theta) + \ln p(\theta) \} = \argmax_{\theta} \{ p(x|\theta)p(\theta) \} \]

となり、(2)が得られます。

次に最尤推定を見てみましょう。最尤推定をベイズ推論の近似であるとすると、事後分布$p(\theta | x)$と事前分布$p(\theta)$に次のような仮定を置いていることになります。

\[ p(\theta | x) \approx q(\theta|\hat{\theta}) = \delta(\theta - \hat{\theta}) \]

\[ p(\theta) = c \]

事後分布が無限大に尖っている分布に対して、事前分布は無限に平坦な無情報事前分布です*3。さて、MAP推定の場合と同様、真の事後分布とのKL divergenceを最小化してみます。

\[ KL(q(\theta|\hat{\theta})||p(\theta|x)) \\ = \langle \ln q(\theta|\hat{\theta}) \rangle_{q(\theta|\hat{\theta})} - \langle \ln p(\theta|x) \rangle_{q(\theta|\hat{\theta})} \\ = \ln q(\hat{\theta}|\hat{\theta}) - \ln p(x|\hat{\theta}) - \ln p(\hat{\theta}) \\ = - \{ \ln p(x|\hat{\theta}) \} + c \]

ここでも$\hat{\theta}$に無関係な項は$c$にまとめました。この式を$\hat{\theta}$に関して最小化したいので、

\[ \theta_{ML} = \argmin_{\hat{\theta}} \{KL(q(\theta|\hat{\theta})||p(\theta|x)) \} \\ = \argmax_{\theta} \{ \ln p(x|\theta) \} = \argmax_{\theta} \{ p(x|\theta) \} \]

となり、(1)が得られます。

以上のように、ベイズ推論の観点からすると、MAP推定や最尤推定は、事前分布や事後分布に対して非常に極端な分布を仮定して推定していることになります。

・MAP推定：

　事後分布を無限に尖った確率分布であると仮定

・最尤推定：

　事後分布は無限に尖った確率分布であり、事前分布は無限に平坦な確率分布であると仮定

事実、KL divergenceに基づく近似の観点からすると、どちらの推定方法も真の事後分布から「無限大に遠い」ような近似事後分布を計算していることになってしまっています。
ただし2つの手法にも利点はあります。それはデルタ関数を使うことによって難しい積分計算（期待値計算）を回避することができるという点です。積分さえ回避してしまえば、あとはKL divergenceを勾配法などを用いて最小化すればOKです。勾配法を使うには微分を要しますが、一般的に積分をするよりははるかに簡単です。

*1:ここでは推定（estimation）と推論（inference）を分けています。推定はある特定の値を何らかの方法で求めるのに対して、ここでは推論は確率計算によりある確率分布を求めることを意味しています。ちなみにMAP推定のことをベイズ推定と呼ぶ場合もあり、非常にややこしいことになっています。

*2:事後分布が解析的に求められない場合に限り、変分近似やラプラス近似のような最適化手法を使います。

*3:厳密に言うと、連続値を取る確率分布の場合、定数cを$-\infty$から$+\infty$まで積分すると無限大に発散してしまうため、$p(\theta) = c$は正しい確率分布ではないです。ここでは便宜上、$\theta$に事前に偏りがないという意味で定数$c$を置いています。