多変量ベルヌー分布の変数間の相関を考慮した同時確率分布を考える

タイトルがやたら難しい!おそらく、3、4日前の自分がこのタイトルを見たら、記事の中身を決して読まなかっただろうと思う(笑)

はじめに:それはベルヌーイ分布だった

最近、何か色々考えたりプログラムしていて、自分のやろうとしていることを、後で、その理論的枠組みを知ることが多くなった。もともと、中野馨氏のアソシアトロン をやっていて、そこでベイズ定理を用いて概念分類をしているのが面白くて、やってみたら、後で、それが「単純ベイズ分類器 Naive Bayes Classifier」というものであることを知って、色々な知識が増えた。それをやったら、今度は、「単純」ではダメだ。ピクセルの変数間の相互関係を考慮しないとMNISTの正解率を84%以上には上げれないと思って、どうしたらいいんだと色々調べたら、そのために必要なことが、実は多変量ベルヌーイ分布の相関を考慮したモデルの問題なのだということを知った。

そもそも、MNISTデータのテスト正解率をあげるのに、なぜ相関を考慮する必要があるかをまず書いておこう。単純ベイズ分類器の場合、一つのピクセルの事前分布を調べるが、その変数の分布は独立であると想定する。MNSITデータは、一つの手書き数字に、28X28=748ピクセルを用いている。Naiveの場合、その第iピクセルが0であるか1であるか(MNISTは、0から255までの数値になっているが、経験的に、これを0,1に変換しても結果はそれほど変わらない)になっているが、その0か1かの値をとる変数は他のピクセルの変数とは独立に与えられる確率変数であると想定するのである。しかし、実際は、数字が形をなすということは、ピクセルのどこかの値が1になることと、別のピクセルの値が1になることが、連動していることを意味している。連動しているからこそ、数字として人間が認識できるのである。だから、相関を考慮することはどうしても必要なのである。

また、ピクセルが独立しているなどという想定は、脳の機構とは全く違ったものになる。ニューロンが他のニューロンとは独立に発火するというのは、脳的ではないし、面白くない。相関を考慮するということは、ニューロン同士の関係性をある程度シミュレートすることにもつながる。

というわけで、ここにたどり着いたのである。この「多変量ベルヌー分布の変数間の相関を考慮した同時確率分布の導出」という問題は、一つのピクセルの値をベルヌーイの二値の分布ではなく、正規分布と仮定すること、さらには、多くのピクセルの変数を相関のある多変量正規分布と仮定することよりも、問題が複雑なのではないかと思う。

そもそも、1,0の二値の分布が独立でない場合の同時分布とはなんだろうか、そこが全く掴めなかった。ネットでも、わかりやすく解説しているところがなかなか見つからなかった。しかし、自分の理解した範囲で、以下で解説してみよう。

2変数のベルヌーイ分布:共分散がある場合

MNISTは784ピクセルで1文字を表す。が、ここで目一杯単純化して、あるパターン(MNISTの場合は文字)が2ピクセルで表されたとしよう。一つのピクセルが0か1の値を持つとする。このパターンのもとで、二つのピクセルがどのような値を持つかを考えるのである。最大可能状態なパターンは4つしかなく、それをとし、それらが二つのピクセルの値をそれぞれである。一般的に、第1ピクセルの値を確率変数(ベルヌーイ確率変数)としてと表すと考えておいても良い。

今、このパターンのもとで、状態が現れる確率は、の二変数の同時確率であり、それをで表すことにしよう。この同時確率は、一般ベイズ分類器の事前確率になるものであり、私たちが求めたいところのものである。つまり、あるパターンのもとで、状態iが現れる確率がわかれば、状態iが現れた元でその原因となるパターンの確率がベイズ定理によって推計することができる、元のパターンを推定することができるわけである。

今、このパターンのもとで、第1ピクセルが1の値を持つ確率をとしよう。確率演算子を用いて、次のように書くこともできる。すなわち、普通の期待値演算子を使うとである。同じく、第2ピクセルが1になる確率をとしよう。

今もし、二つのピクセルの値の確率変数が、相互に独立だとしたら、先の同時確率は簡単に求めることができる。すなわち、状態の場合、二つのピクセルはともにゼロであるから、同時確率は、

となる。同様に、状態2の場合は、第1ピクセルが1で第2ピクセルは0であるから

となる。一般に、

と書けば、すべての状態を表していることになる。独立の場合は、これを使えばベイズ分類が可能になり、それは、先の記事で私がやったNaive Bayes Classifierに他ならない。

では、独立性がなく、二つのピクセルの値の間に相関がある場合はどうなるだろうか。それを以下で考えてみよう。

今、第1ピクセルが1を持つ状態は、であるから、それらの状態の同時確率と、この状態の確率(周辺確率と呼んでもいい)との関係は次のようになる。(括弧つきの数字は、式の番号)

(1)

は、二番目のピクセル(確率変数)が1になる確率であり、それは次のようになる。

(2)

さらに、二つの確率変数の間の共分散が存在し、それをであらわそう。この時、状態の同時確率と共分散の間には次のような関係が成立する。

(3)

私と同じくらいの理解力の人間では、上の式を見てすぐには理解できないと思うので、少し説明しよう。

共分散の普通の定義は、確率変数の共分散は、その期待値をそれぞれとして、期待値演算子 を用いて、と表される。言葉で表せば、は、状態の値 に、この状態が現れる確率をかけたものである。

今、(3)式をvの項が4つあると見て、一番最初の項をみてください。これは状態は変数が0,0の値を持っているので、分散の定義では、v_{1}$をかけたもの、すなわち、となったものなのである。他の三つの項も、分散の定義に沿っていることは、お分かりになるだろう。それらの合計が、我々の求めるに他ならない。

さて、もう一つ式がある。これが一番簡単だ。すなわち、同時確率も確率だから、それらを合計したものは1でなければならないということである。

(4)

我々の目的は、 各状態の同時確率の4つの値を求めることだが、それに対応する4つの式を得たことになる。当然、およびは与えられているものとする。実際それは、MNISTの場合、教育用データで標本値を求めることができるので、心配はいらない。

解くための学力は、中学あるいは高校の数学力があればいいだろう。が、4元連立方程式を解くというのはとても面倒くさい。十代の頃は、いつもケアレスミスで、間違ってしまった。でも、やらなければ前に進めないので、やりました。結果は次のようなものとなる。




結果を見て、そのわかりやすさに感激した。つまり、共分散のある場合の同時確率は、共分散がない場合、すなわち、相互に独立の場合の同時確率に、分散を足すか引くかしたもの。両方が1ないしは0の時は、共分散を足し、違う場合は引くということだ。つまり、共分散がある状態は同時に現れたり、消えたりする確率が高くなる。

結果を見て、なるほどと思った。

ただ、これだけでは、一般にどうなるかが予想できない。これは2ピクセルの場合である。しかし、MNISTは、784ピクセルもあるのだ。

もちろん、すぐに一般的な場合を考えても良いのかもしれない。計算手順はだいたいわかったので。しかし、やはりそれは無謀だ。まず、3ピクセルの場合をやってみようと思った。

3変数のベルヌーイ分布:共分散がある場合

3変数の場合も、基本同じだろう。と、思った。まず、テーブルを書いておく。

まず、式を作っていく。一番わかりやすいのは同時確率の合計が1という である。さらに、周辺確率と同時確率との関係から、3本の式ができる。また、相関係数と確率の関係から式が3本できる(相関係数は、対象であることに注意)。あれ、7つしかできない。未知数である同時確率は8個あるのに、式が一つ足りないのだ。これには、かなり悩んだ。訳がわからなかったが、いろいろ調べてわかった。

の三つの確率変数の間の同時共分散を考慮しなければならないのである。今この分散をとしよう。すなわち、である。

この8本の式を全部書くと、ノート1ページを埋めるくらいになった(笑)。さらに、このサイトに書く元気は湧いてこない。かといって、何も書かないわけにもいかないので、行列表現で、この8元連立方程式を書いておこう。ただし、要素の中の式が長くならないようにとしておく。

式はできた。解くための必要条件、式の数=未知数の数、も満たしている。が、さすがにこれを記号的に解く気にはなれない。コンピュータに頼るしかない。記号が入っていない、変数以外は実際の数値ならば、EXCElで簡単に解くことができるが、ここでは係数などが記号のままとかなければならないのである。数式を記号のまま処理して、連立方程式を解けるシステム、まず、頭に浮かんだのはMathematicaで、かつてライセンスを取っていた。ほとんど使わなかったが、今更、アップグレードやら何やらの手続きも面倒だ。トライアル期間を利用して解くという手もあるかと思ったが、そもそも使い方がよくわからない。

そこで、Maximaというのが無料で、よく使われているということを知って、これをMacにインストールした(サイト)。すぐに解くことができた。ただし、式が少し間違っていると、答えを出さなかったり、とんでもなく長い答えを出してくるので、正確に式を入れなければならない。

solve([v1+v2+v3+v4+v5+v6+v7+v8=1,
            v5+v6+v7+v8=p1,
            v3+v4+v7+v8=p2,
            v2+v4+v6+v8=p3,
            p1*p2*(v1+v2)-p1*(1-p2)*(v3+v4)-(1-p1)*p2*(v5+v6)+(1-p1)*(1-p2)*(v7+v8)=s12,
            p1*p3*(v1+v3)-p1*(1-p3)*(v2+v4)-(1-p1)*p3*(v5+v7)+(1-p1)*(1-p3)*(v6+v8)=s13,
            p2*p3*(v1+v5)-p2*(1-p3)*(v2+v6)-(1-p2)*p3*(v3+v7)+(1-p2)*(1-p3)*(v4+v8)=s23,
            -p1*p2*p3*v1+p1*p2*(1-p3)*v2+p1*(1-p2)*p3*v3-p1*(1-p2)*(1-p3)*v4
            +(1-p1)*p2*p3*v5-(1-p1)*p2*(1-p3)*v6
            -(1-p1)*(1-p2)*p3*v7+(1-p1)*(1-p2)*(1-p3)*v8=theta]
,[v1,v2,v3,v4,v5,v6,v7,v8]);

(%o19) [[v1 = (- theta) + s23 + p1 ((- s23) + p3 + p2 (1 - p3) - 1) + s13
 + p2 ((- s13) + p3 - 1) + s12 + p3 ((- s12) - 1) + 1, 
v2 = theta + p1 (s23 + p2 p3 - p3) - s23 + p2 (s13 - p3) - s13 + p3 (s12 + 1), 
v3 = theta + p1 (s23 + p2 (p3 - 1)) - s23 + p2 (s13 - p3 + 1) + p3 s12 - s12, 
v4 = (- theta) + s23 + p1 ((- s23) - p2 p3) + p2 (p3 - s13) - p3 s12, 
v5 = theta + p1 (s23 - p3 + p2 (p3 - 1) + 1) + p2 s13 - s13 + p3 s12 - s12, 
v6 = (- theta) + p1 ((- s23) - p2 p3 + p3) - p2 s13 + s13 - p3 s12, 
v7 = (- theta) + p1 (p2 (1 - p3) - s23) - p2 s13 - p3 s12 + s12, 
v8 = theta + p1 (s23 + p2 p3) + p2 s13 + p3 s12]]

上半分に、solveというコマンドで解くための式を定式化している。下半分の、Maximaのコマンドプロンプト(%o19)の後に、 '[[' と ']]'に挟まれてからまでの答えが書き出されている。素晴らしい!!なお、thetaはのことである。

簡略化が、不完全なので、ぱっと見よくわからないが、整理すると次のようになる。ただし、であることに注意する。2次元の場合の自然な拡張になっていることも直ちに読み取ることができる。








この同時確率の正しさを検証するものとしては、Jozef L. Teugels氏が1990年にJournal of Multivariate Analysis, Vol.32,誌に発表した論文"Some Representations of the Multivariate Bernoulli and Binominal Distributions" (pp.256-268)がある(サイト)。そこに、私のものと同じように、2次の場合と3次の場合の式が書かれている。ただし、3次の場合の式の並びと、ここでの式の並びは対応していない。それを適切に対応させれば、両者が一致していることを確認できるだろう。実は、先のが必要になることは、この論文に教えてもらった。

ここまでくると、一般的な多数のベルヌーイ変数がある場合の同時確率が予想できる。ただし、その一般形を書いた式は、私は見つけていない。いや、書き方が大変なことは予測できる。だからこそ、上記のTeugels氏は、クロネッカー積というものを使わざるを得なかったのだろうと思う。もし、どなたか、多変量ベルヌーイ分布の同時確率の一般系の書き方がわかる方がおられたら是非教えていただきたい(ツイッター @wassiisg まで)。

ただし、共分散を考慮したベルヌーイ分布を使って、一般的な(Naiveではない)ベイズ分類器を作る場合、必ずしも一般式は必要ではないというか、一般的な式のまま計算はできないと思う。どのように使うかは、次以降の記事で書きたい。