統計学における平均場近似

Stat

統計学における変分法に興味を持ってしまったため、変分ベイズの前に必要な平均場近似について勉強したことをまとめます。参考書は主に渡辺澄夫先生の解説文で、残りの部分は『計算統計 I 』や樺島先生の論文を参考にしました。『パターン認識と機械学習（下）』にも説明があります。


『物理学者でない人のための統計力学（渡辺澄夫先生）』
（参考にさせて頂いた資料は直接リンクできなかったので、代わりに同じ内容とおもわれるもののリンクをのせています。検索するとより詳細な資料が見つかります。）

『グラフィカルモデルと平均場近似（樺島祥介先生）』
（計算統計Iの平均場近似の箇所とほとんど同じです。）

計算統計 I―確率計算の新しい手法 (統計科学のフロンティア 11)

計算統計 I―確率計算の新しい手法 (統計科学のフロンティア 11)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

機械学習とか統計学について素人なので嘘八百書いていたらごめんなさい。

0.　何をしようとしているのか？

統計モデルにおける計算問題には2種類ある。平均等統計量の算出（順問題）とパラメータの推定（逆問題）。不完全データが存在するモデルでは、この両方が問題になることがある。例えば、EMアルゴリズムによる推定など。

統計モデルの説明変数が非常に多い場合、順問題である平均の算出すら計算量的に大変であることが多い。このような時には近似的な計算法が必要になる。近似計算法にも、確率的近似法、決定論的近似法という2つの系統がある。前者ではマルコフ連鎖モンテカルロ法が知られ、後者では平均場近似や変分ベイズ法が知られている。

今回はこの平均場近似を取り上げる。平均場近似というのは、物理学の一分野である統計力学で生まれ育ったものであるため、統計力学の世界を振り返る。

1.　物理における平均場近似

僕も物理専攻でなかったので統計物理とかよくわからないし、今のところそれを学ぶ意欲がなかったりするので、数学の部分だけ渡辺先生の資料を読んで勉強します。

1.1.　無限次元空間

統計物理では気体中の分子運動等を対象としているため、対象の次元が非常に大きい（分子の個数をNとすると、位置で3N、運動量で3N、計6N次元。Nは大体6.02×10^23）。もうこれはある種の無限次元と思いましょう、というところが出発点。ということで、次のような離散無限（可分なヒルベルト空間）の世界を考える。

　[tex:R^\infty=\{(x_1, x_2, ... , x_n, ...); -\infty

1.2.　無限次元空間における確率分布 (1) - 標準正規分布

上記で定義されている無限次元空間において、確率分布を定義したい。
まず用語の定義。確率分布のサポートとは。（密度関数が存在すれば、その関数の台（サポート）のこと。）

　P(A) = 1となる最小の閉集合Aを確率分布Pのサポートと呼ぶ。

次のようなR^n上の標準正規分布をまず取る。

　p_n(x)=p(x_1, x_2, ... , x_n) \propto \exp(-||x||^2/2)

これを使って、R^{\infty}上の分布で、n次元までは正規分布、その他の変数は適当な分布関数に従う、次のようなものが構成できる。

　p_n^{*}(x)=p_n(x_1, x_2... , x_n) \phi(x_{n+1}, x_{n+2}, ...)

このような確率分布の列

　p_1^{*}(x), p_2^{*}(x), p_3^{*}(x), ...

を取ると、この場合は\phi(y)の取り方によらずある確率分布に収束する。

　\lim_{n \to \infty} p_n^{*}(x) = p^{*}(x) \propto \exp(-||x||^2/2)

これを無限次元空間における標準正規分布と定義する。重複対数の法則より、この分布のサポートはH(0)に含まれている。

1.3.　無限次元空間における確率分布 (2) - 混合正規分布

正規分布の場合は上手くいった、有限次元分布の極限としての無限次元空間上の分布だが、混合正規分布では収束先は複雑になる。正規分布の場合と同様に、有限次元の混合正規分布をまず取る。

　1=(1,1,...,1) \in R^n
　q_n(x) \propto \exp(-||x||^2/2)+\exp(-||x-1||^2/2)

これを使って、R^{\infty}上の分布で、n次元までは有限次元混合正規分布、その他の変数は適当な分布関数に従う、次のようなものが構成できる。

　q_n^{*}(x)=q_n(x_1, x_2... , x_n) \phi(x_{n+1}, x_{n+2}, ...)

このような確率分布の列

　q_1^{*}(x), q_2^{*}(x), q_3^{*}(x), ...

を取ると、この場合はq_n^{*}(x)のサポートによって収束先の確率分布q^{*}(x)が違う。

  1. q_n^{*}(x)のサポートがR^{\infty}全体の場合、q^{*}(x) \propto (1-a) \exp(-||x||^2/2) + a \exp(-||x-1||^2/2)
  2. q_n^{*}(x)のサポートが常にH(0)に入っているとき、q^{*}(x) \propto \exp(-||x||^2/2)
  3. q_n^{*}(x)のサポートが常にH(1)に入っているとき、q^{*}(x) \propto \exp(-||x-1||^2/2)
1.4.　分配関数と自由エネルギー

再び混合正規分布が対象だが、今度は混ぜる正規分布の分散および混合比率が違うモデルを考える。

　q_n^{*}(x) = Z(n)^{-1} \{ n \exp(-||x||^2/2) + \exp(-||x-1||^2/4) \}
　Z(n)=\int \{ n \exp(-||x||^2/2) + \exp(-||x-1||^2/4) \} dx

このZ(n)は分布の全積分が１になるような正規化定数で、分配関数という名称がついている。q_n^{*}(x)の収束先のサポートはどうなるだろうか？　これが分配関数の漸近挙動でわかるというのが本節の結論。この例の場合、分配関数は簡単に計算できて、次のようになる。

　Z(n) = n \( 2 \pi \) ^ {n/2} + \( 4 \pi \) ^ {n/2} \simeq \( 4 \pi \) ^ {n/2}

したがって、H(1)より外では、Z(n)の寄与によって無限大に飛ばしたときに消えてしまうことがわかり、収束先のサポートがH(1)に入っていることがわかる。

　q^{*}(x) \propto \exp(-||x-1||^2/4)

同様のことは、自由エネルギーを用いた議論でも確認できる。F(n)=-logZ(n)を自由エネルギーといい、その部分的な自由エネルギーを以下で定義する。

　f_1(n)=-\log{ \{ \int n \exp{(-||x||^2/2)} dx \} }
　f_2(n)=-\log{ \{ \int n \exp{(-||x-1||^2/4)} dx \} }

すると次の式が成り立ち、

　F(n)=-\log{ \{ \e{ \{ -f_1(x) \} } + \e{ \{ -f_2(x) \} } \} }

Laplace近似より次がわかる。

　F(n)=\min{ \( f_1(n), f_2(n) \)} \simeq f_2(n)

このことから、「無限大への極限の分布は、自由エネルギーを最小にする状態を持つ」という原理が確認できた。

1.5.　平均場近似

分配関数が計算できれば（その漸近挙動から）無限次元空間上の確率分布の解析ができるということだが、厳密に分配関数がわかるモデルは少ない（可解模型と呼ばれるらしい。無限次元Lie代数とかなんとか）。そのような場合に、結果が正しいかわからないがもしかしたら何かわかるかもしれない近似法として、平均場近似と呼ばれるものがある。

ある確率分布p(x)を持つモデルがあるが、変数の次元が高すぎるため平均等の計算が容易にできない場合に、計算が容易な確率分布の族Q= \{ q(x) \} で元の確率分布p(x)を十分よく近似できたとする。そのとき、このq(x)を用いて平均等の計算を行うということを「平均場近似」と呼んでいる。

分布の近さをはかるために、次の相対エントロピーを最小にするものを選ぶことがよく行われる。この相対エントロピーは統計学ではカルバック・ライブラー情報量と呼ばれていることが重要である。

　\int q(x) \log{\frac{q(x)} {p(x)}} dx

2.　統計学における平均場近似

原理的な部分は今までの説明で終わっている。残りの部分で、グラフィカルモデルの場合への応用を見る。

2.1.　グラフィカルモデル

グラフィカルモデルの説明については省きます。余裕があれば自分の復習のためにいつか書くかもしれません。ガウシアングラフィカルモデルを想像しながら以下を読んでください。

N次元の状態変数xに対して、ガウシアンの時の２次形式に相当する、次のポテンシャル関数を用いて確率分布が定まるモデルを考える。

　\phi(x) = \sum_{\mu}{ \phi_{\mu}(x_{\mu})}
　P(x)=\frac{ \exp [ \phi(x) ]} {\sum_{x^{\prime}} {\exp [ \phi(x^{\prime}) ] }}

ポテンシャル関数を構成する各関数の引数に現れる添字集合をクリークという。変数間の依存関係を、このクリークの（位相）幾何的な構造で表しているのがグラフィカルモデル、と理解してます。

グラフィカルモデルの場合、変数の依存関係によって多重積分が何回も登場し、計算量が指数関数的に増大する。これを平均場近似を用いて少ない計算量に落とそうというのがこれからの目的。

2.2.　各変数が独立なグラフィカルモデル

最も簡単なグラフィカルモデルである、各変数が独立なモデルを考える。このとき分布関数は次のように表される。

　P(x)=\prod_{i=1}^{N} \frac{ \exp [ \phi_i(x_i) ]} {\sum_{x_i^{\prime}} {\exp [ \phi_i(x_i^{\prime}) ] }}

この場合、積分は１重積分しか登場しないため、計算量はO(N)に過ぎない。このように高速計算が可能なグラフィカルモデルは、他にも１次元鎖の場合やジャンクションツリーの場合が知られている。

2.3.　ナイーブ平均場近似

高速計算が可能なモデルの族（独立、１次元鎖、ジャンクションツリーなどの合成）から１つモデルを取り出し、その分布関数をQ(x)とする。考えているモデルP(x)とこのQ(x)の近さを測る量として、次のカルバック・ライブラー情報量が有名である。

　D(Q|P)=\sum_x {Q(x) \log \frac {Q(x)} {P(x)}}

このカルバック・ライブラー情報量を利用して、モデル族として変数が独立な場合を取った場合に、最もよい近似をしているモデルを求める計算アルゴリズムは次のようになる。

　Q(x)=\prod_{l=1}^N Q_l(x_l)

　D(Q|P)=\sum_x Q(x) \log{\frac{Q(x)}{P(x)}}
　=-\< \phi(x) \>_Q + \sum_{l=1}^N \sum_{x_l} {Q_l(x_l) \log{Q_l(x_l)} } + \log \sum_x \exp [ \phi(x) ]

これを最小にするようにQ_l(x_l), \( l=1, 2, ... , N \)を決める。ただし、平均操作の記号として以下を使っている。

　\< \cdot \>_Q = \sum_x \prod_{l=1}^N Q_l(x_l) \( \cdot \)
　\< \cdot \>_{Q \backslash l} = \sum_{x \backslash l} \prod_{j \neq l}^N Q_j(x_j) \( \cdot \)

カルバック・ライブラー情報量を最小にするのは、次の分布である。（ここよくわからず…）

　Q_l(x_l)=\frac{\exp [ \< \phi(x) \>_{Q \backslash l} ] } {\sum_{x_l} \exp [ \< \phi(x) \>_{Q \backslash l} ] }

これは反復法によって解くことができる。

3. まとめ

基本的な原理は、高速に計算可能なモデルがあれば、考えているモデルとカルバック・ライブラー情報量を最小にするモデルをとって平均計算を代用すればよい。ただし、答えが近いという保証はできない。研究上は、実際に可能なモデルと平均場近似の場合の結果を対比させることで手法の是非を論じているとのこと。

グラフィカルモデル以外にもいろいろな応用があるらしい。様々な文献がある。