作って遊ぶ機械学習。

～基礎的な確率モデルから最新の機械学習技術まで～

第一線のAI研究者が注目する最新機械学習技術6選（NIPS2015招待講演より）

まとめ

ちょっと前になりますが、昨年12月に行われた機械学習のトップカンファレンスであるNIPS2015の講演ビデオが上がっているようなのでチェックしてみました。今回ご紹介するのはケンブリッジ大学のZoubin Ghahramani教授の研究です。

datasciencereport.com

同教授は今後excitingな機械学習の基礎・応用に関する取り組みとして次のような６テーマを紹介しています。

 

・Bayesian Nonparametrics

ベイジアンノンパラメトリクスでは、無限次元を持つモデルを仮定することにより、データ量に応じて適切なモデルを学習することができます。こういったモデルは関数上の確率分布（確率過程）を考慮することによって実現できます。代表的な例はガウス過程、中華料理店過程、インド料理過程などです（ご飯ばっかりですね）。

ガウス過程は、回帰や識別、ランキングや次元削減に使われています。またある種のニューラルネットワークがガウス過程の一例であることが1994年に示されているそうです。

・Probabilistic Programming

確率モデルを構築し、推論アルゴリズムを導き、実装をする・・・といったプロセスはものすごく時間がかかって大変です。Probabilistic programmingは確率モデル（あるいは確率的にデータを生成するプロセス）を記述するための言語で、推論はサンプリングや変分近似で勝手にやってくれます。代表的なのはBugsやStanです。

・Bayesian Optimization

日本語だと「ベイズ的最適化」の名前でコンセンサスが取れていたと思います。ガウス過程を使い未知の関数の最適値を探す手法です。未知の関数の評価にとても時間orコストがかかるような状況を仮定しており、いかに少ない評価回数で最適値を探せるかがポイントになります。

・Data Compression

シャノンの定理によると、すべての圧縮アルゴリズムは確率モデルに基づくそうです。ベイジアンノンパラメトリクスを使ったアダプティブな圧縮アルゴリズムが近年非常によい圧縮性能を出しているようです。

・Automatic Statistician

データを食わせるだけで何かしらレポートを返してくれるような野心的なアルゴリズムの話です。アイデアは非常にシンプルで、いくつかのカーネルを組み合わせて周辺尤度を評価することによって、データをうまく説明するような適切なモデルを探索します。

・Rational Allocation of Computational Resources

大量のデータと複数のモデルがある中、リソース（時間、CPU、メモリ、ディスクサイズ）が限られた状況下でいかに最適な機械学習システムを構築するか、という課題設定です。これを不確実性がある中での逐次的な意思決定問題として扱います。

 

・個人的な雑感

Bayesian nonparametricsは従来のベイズ学習に対する正当な拡張であるため、今後も基礎的な研究は発展していくかと思います。そろそろ大きなアプリケーションがほしいところですが、ディープラーニングなどの成功しているモデルもかなりの部分はこれに置き換わり、より洗練されたアルゴリズムになるのではないかと予測しています。

Probabilistic programmingは基本的に推論性能はあまりよくないようです。僕は使ったことがないのですが、とりあえず複数のモデルを素早く試したいときなんかには便利だと思います。普通にモデル構築→アルゴリズム導出→実装ってやってるとそれだけで１日使いますからね・・・。

ベイズ的最適化は非常に概念もわかりやすく応用先も広いのでぜひ多くの人に使っていただきたいと思います。日本の京などのでっかいコンピュータで何かシミュレーションを行う場合は、Bayesian optimizationで全体を包んであげると解探索が早くなるかもしれません。あまりいないと思いますが、例えば遺伝的アルゴリズムとかを使っている人は早々にこちらへの置き換えを検討した方が良いと思います。

データ圧縮も面白い領域です。結局のところ、jpegに代表されるように、非可逆圧縮というのはデータの中の捨ててもいい情報と残さなきゃいけない情報をうまく見つけることが基本になっています。これは機械学習における特徴量抽出といっしょです。画像だったら、人間が認知に必要でない情報は捨ててもいい可能性が高いと言えるわけですね。データ圧縮もかなりドメイン依存が強いかと思いますが、汎用的なアルゴリズムでどこまでいけるのかは興味深いですね。 （追記：すみません、招待講演の内容では「非可逆圧縮」ではなく、情報損失のない「可逆圧縮」の話だったようです。データに関する分布をベイジアンノンパラメトリクスを使って推定することにより、データを表現するための符号長を統計的に短くすることができるという話です。研究のポイントとなる点は、zipなどで使われている性能の良い可逆圧縮アルゴリズムと同等のものがベイジアンノンパラメトリクスに基づいたアルゴリズムで解釈できることを示し、かつ性能も向上させたことです。すごいですね。）

Automatic statistiicianはデータサイエンスにおける究極的な課題の１つなんじゃないかと思います。ただ、質問でもありましたが、決してデータサイエンティストの職業を奪うようなものではなく、あくまで便利なツールであるという位置づけだそうです。

リソースアロケーションの話は実応用で非常に重要です。まぁそれ以前に日本だったらちゃんと複数の確率モデルを作れる人がいないと話にならないですが・・・。

 

以上です。

多峰性の事前分布ってどうやって作るの？

グラフィカルモデル

ベイズ学習では、複雑なモデルにおけるパラメータの学習を効率的に行うために、しばしば観測モデルに対する共役な事前分布を仮定します。例えばベルヌーイ分布のパラメータの事後分布を推定するために、事前分布をベータ分布にしたりします。ウェブを検索すると様々な確率分布のパラメータに対応する共役事前分布のリストを見つけることができますが、そもそも、データを表現する観測モデルが決まってしまうと自動的にそのパラメータの事前分布が決まってしまうのは何だか窮屈な感じがしてしまいます。

今回はグラフィカルモデルを工夫することによって、もっと自由な、例えば多峰性の事前分布などを簡単に作れることをお話ししたいと思います。

 

[必要な知識]

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

 

さて、改めてベイズ学習におけるパラメータの学習を書いてみると次のようになります。

{ p(\theta| x) \propto p(x|\theta)p(\theta) }

ここで{x}は観測データ、{\theta}はモデルのパラメータです。さらに{p(\theta|x)}はパラメータの事後分布、{p(x|\theta)}は尤度関数、{p(\theta)}はパラメータの事前分布と呼ばれています。ベイズ学習を言葉に翻訳すると、

「パラメータに関する我々の事前知識{p(\theta)}は、観測データに対するモデルの尤もらしさ{p(x|\theta)}を通して、事後の知識{p(\theta|x)}に変換される」

というプロセスになります。

で、今回は尤度関数を計算するモデルとして次のようにガウス分布を仮定します。

{ p(x| \theta) = \mathcal{N}(x| \mu, \Sigma_x ) }

ここで今回推定したいパラメータは平均値{\theta=\mu}であり、分散{\Sigma_x}に関しては既知で固定であるものとします。

 

・単峰性事前分布によるガウス分布の推定

まず始めに、ガウス分布の平均値パラメータに対する一番シンプルな事前分布の選び方は、共役事前分布であるガウス分布を使うことです。*1

{ p(\mu) = \mathcal{N}(\mu | m, \Sigma_{\mu}) }

ここで{m}と{\Sigma_{\mu}}はそれぞれ平均値と分散パラメータであり、{\Sigma_x}と同様固定値であるものとします。

このモデルに対して{N}個のデータを観測した後の事後分布は、グラフィカルモデルで表現すると次のようになります。

f:id:sammy-suyama:20160220183059p:plain

さて、さっそく事後分布を計算してみましょう。事後分布はベイズの定理から次のように計算できます。

{ p(\mu|x) \propto p(x|\mu)p(\mu) }

両辺対数を取って計算を進めることにします。{\mu}に関する分布を求めたいので、無関係な項はすべて定数項{c}に吸収させてしまいます。

f:id:sammy-suyama:20160220215841p:plain

結果はやはりガウス分布になり、新たに文字{\hat{m}}と{\hat{\Sigma}_{\mu}}を導入すれば次のように書けます。

f:id:sammy-suyama:20160220220026p:plain

・多峰性事前分布によるガウス分布の推定

さて、今度は先ほどの単純な共役事前分布と違い、{\mu}に対して{K}個の異なる平均値の分布からなる、多峰性の事前分布を構築してみましょう。感覚的な例をいうと、推理小説で「犯人が一人なのはわかっているけど、容疑者が複数いる」というような状況設定です。

これは次のグラフィカルモデルで示すように、新たな潜在変数{z}を導入することによって実現できます。

f:id:sammy-suyama:20160220183117p:plain

まず始めに新しく登場した{z}に対する確率分布を定義してみましょう。ここではもっとも単純な離散変数の分布であるカテゴリカル分布を使います。

{\begin{eqnarray} p(z) = Cat (z|\pi) = \prod_{k=1}^{K} \{ \pi^{(k)} \}^{z^{(k)}} \end{eqnarray} }

ただし{\sum_{k=1}^{K} \pi^{(k)}=1 }です。さらに、{\mu}の分布は次のように{z}の条件付き分布になります。

{ \begin{eqnarray} p(\mu|z) = \prod_{k=1}^{K} {\mathcal{N}(\mu| m^{(k)}, \Sigma_{\mu}^{(k)})} ^{z^{(k)}} \end{eqnarray} }

このような{z}に対する記法は1 of K表記というように呼ばれたりします。例えば{z=(0,0,1,0)}のように、{z}の{K}個ある要素のうち一つだけが１を取ることによって、どの分布が使われるのかを指示（indicate）することができます。

それでは事後分布を計算してみることにしましょう。

{p(\mu, z| x) \propto p(x,\mu,z) = p(x|\mu)p(\mu|z)p(z) }

今回は{\mu}と{z}の2つの事後分布を求めなければなりませんが、これはそれぞれ別々に解析的に計算できます。まず上の数式から{\mu}に関わる項のみを取り出してみましょう。例によって対数計算を行うと、

{ \ln p(x,\mu,z) = \ln p(x|\mu) + \ln p(\mu|z) + c }

となります。先ほどの単峰性の例と同じように{\mu}のみに注目して計算を行うと、結果として{\mu}の事後分布は{z}に依存する形となり、次のような（条件付き）ガウス分布になります。

f:id:sammy-suyama:20160220220259p:plain

さて、今度は{z}の事後分布も求めてみましょう。これは{\ln p(z|x) = \ln p(\mu|z) - \ln p(\mu|z,x)}を考慮すれば簡単に計算できます。

{ \begin{eqnarray} \ln p(z|x) = \sum_{k=1}^{K} z^{(k)} (-\frac{1}{2} {\mu^{(k)}}' \Sigma_x^{-1} \mu^{(k)} -\frac{1}{2} \ln |\Sigma_{\mu}^{(k)}| \\ +\frac{1}{2} {\hat{\mu}^{(k)}}' \hat{\Sigma}_x^{-1} \hat{\mu}^{(k)} +\frac{1}{2} \ln |\hat{\Sigma}_{\mu}^{(k)}| +\ln \pi^{(k)})+c \end{eqnarray} }

{z}に掛け算されている部分が事後分布のパラメータになります。書き直すと次のようになります。

{ \begin{eqnarray} p(z|x) = Cat (z|\hat{\pi}) = \prod_{k=1}^{K} \{ \hat{\pi}^{(k)} \}^{z^{(k)}} \end{eqnarray} }

{\hat{\pi}}は次のように計算されます。

{ \begin{eqnarray} \hat{\pi}^{(k)} \propto exp\{ -\frac{1}{2} {\mu^{(k)}}' \Sigma_x^{-1} \mu^{(k)} -\frac{1}{2} \ln |\Sigma_{\mu}^{(k)}| \\ + \frac{1}{2} {\hat{\mu}^{(k)}}' \hat{\Sigma}_x^{-1} \hat{\mu}^{(k)} +\frac{1}{2} \ln |\hat{\Sigma}_{\mu}^{(k)}| +\ln \pi^{(k)} \} \end{eqnarray} }

実装上は{\hat{\pi}}が足して1になるように正規化してください。

以上、事後分布{p(\mu,z|x)}を見てみると、事前分布{p(\mu,z)}と同じ形状を持っていることになります。ガウス分布の平均値だからといって単純にガウス事前分布を使う必要はなく、今回のような多峰性のガウス分布を事前分布にしても、最適化のような計算をせずに解析的に事前分布と同じ関数形の事後分布が求まるんですね。*2

 

さて、ちょっと思った以上に計算式ばかりの記事になってしまったので今日は一旦ここで切ります。今回はガウス分布の平均値パラメータの事後分布を推定する例を使って、多峰性の事前分布をどのように構築したらよいのかを説明しました。次回は今回紹介した事前分布を使った簡単な実験を行い、ベイズ学習における興味深い特性に関してお話したいと思います。

 

＝＞続きの記事はこちらです。

アニメでわかるベイズ推論によるパラメータ学習 - 作って遊ぶ機械学習。

 

*1:ガウス分布の平均値の共役事前分布は「たまたま」ガウス分布になります。ガウス分布を観測、平均値ともに使っているので気を付けてください。

*2:ちなみにこの多峰性事前分布によるモデルは、混合ガウス分布とは異なります。混合ガウス分布の場合は潜在変数がデータ{x_n}ごとに割り当てられます。今回の例と違い、解析的に求めることは不可能で、変分近似のような近似推論法を使う必要があります。

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

グラフィカルモデル ベイズで解釈

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

 

・教師あり学習（今回）

・半教師あり学習

・共変量シフト

・転移学習

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

 

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

  

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

machine-learning.hatenablog.com

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

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

 

＜有向分離（D-separation）＞

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

f:id:sammy-suyama:20180813144416p:plain

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

 

・具体例）教師あり学習

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

f:id:sammy-suyama:20160213234013p:plain

$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$が観測されます。描き直してみましょう。

f:id:sammy-suyama:20160213230120p:plain

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

\[ 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。

f:id:sammy-suyama:20160213230508p:plain

さて、$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)$の項をダイエット

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

f:id:sammy-suyama:20160213230527p:plain

この項では$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

f:id:sammy-suyama:20160209232422p:plain

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

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

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

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

f:id:sammy-suyama:20160209233138p:plain

はい、こんな感じです。このとき、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

 

f:id:sammy-suyama:20160209233550p:plain

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

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

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

f:id:sammy-suyama:20160209233607p:plain

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

\[ 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

f:id:sammy-suyama:20160209234931p:plain

最後ですが、この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を条件付けしてみましょう。

f:id:sammy-suyama:20160209234943p:plain

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

\[ 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 動く

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

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

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

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

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

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

 

[必要な知識]

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

 

 

ギブスサンプリングでは、推定したい多変数の事後分布{p(z_1, z_2|x)}に対して、下記のような手続きで{n+1}個目のサンプル値{z_1^{(n+1)}}と{z_2^{(n+1)}}を拾ってきます。

f:id:sammy-suyama:20180716154410p:plain

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

 

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

<変分近似>

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

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

<ギブスサンプリング>

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

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

 

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

 

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

f:id:sammy-suyama:20180716154455p:plain

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

では、実際に条件付き分布を求めてみましょう。{x_2^{(n)}}を固定されたサンプル値であるとして、{x_1}の確率分布を求めてみます。f:id:sammy-suyama:20180716154513p:plain

これを正規化することにより1次元のガウス分布が得られます。f:id:sammy-suyama:20180716154530p:plain

ただし、f:id:sammy-suyama:20180716154549p:plain

です。

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

{x_2^{(n+1)}}のサンプルも全く同様なので省略します。

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

f:id:sammy-suyama:20180716154615p:plain

 

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

f:id:sammy-suyama:20160203200821g:plain 

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

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

f:id:sammy-suyama:20160203201223p:plain

 本当は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}$のまま更新されず、平均値だけ更新されていくようなアルゴリズムになりましたね。

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

 

f:id:sammy-suyama:20160131135817g:plain

f:id:sammy-suyama:20160131135416p:plain

上図では、青い楕円が推定したい真のガウス分布で、赤い楕円が推定中の近似分布です（σ=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$が解析に求まってしまいます。今回は実際の多くの応用のように、あえて繰り返しアルゴリズムを使って解を求めます。