はじめに
今回の記事は下記の記事２つの続きに位置しています。
この記事からいきなり入れる人は、既に下記の知見を得ているものだと思われます。そうでない場合は上記の記事を参考にしてください。
- 同時分布の設計とグラフィカルモデルが対応付けられる
- 設計した同時分布を乗法定理で変形し、事後分布を導出できる
- 事前分布を潜在変数・パラメータの定義域に値を持つ確率分布で設定できる
混合モデル
単一分布おさらい
前回の記事で、正規分布しているように見える$N$ 個のデータ $X ^ N$ に対するモデリングを
$$ p(X ^ N, \mu, \sigma \mid \mathbf \alpha) = p(\mu \mid a, b)p(\sigma \mid b,c)p(X ^ N \mid \mu, \sigma) $$
と行いました。ここで $\mathbf \alpha = (a, b, c, d)$ はハイパーパラメータを格納した確率変数ベクトルです。ハイパーパラメータを何らかの方法で推定しないのであれば、モデルの中に明示的に書いておく必要は無いので、
$$ p(X ^ N, \mu, \sigma) = p(\mu)p(\sigma)p(X ^ N\mid \mu, \sigma) $$
とより簡略された式を使ってみていきましょう。
混合モデル
混合モデルとは、何らかの単純な分布が幾つも集まってできているような分布を指します。他の単純な分布とやらを $p$ として一旦グラフィカルモデルで視覚的に混合モデルの例を見てみましょう。まず簡単のため確率変数 $X$ が１個だけ観測される様子をグラフィカルモデルで描いてみます。
さてこのように書かれたグラフィカルモデルは下記のような式にひとまずは書き下すことができます。
$$ p(X, Z, \mu _ 1, \sigma _ 1, \mu _ 2, \sigma _ 2, \mu _ 3, \sigma _ 3) = p(\mu _ 1)p(\sigma _ 1)p(\mu _ 2)p(\sigma _ 2)p(\mu _ 3)p(\sigma _ 3)p(Z)p(X \mid Z, \mu _ 1, \sigma _ 1, \mu _ 2, \sigma _ 2, \mu _ 3, \sigma _ 3) $$
このように同時分布の設計を行った場合、観測変数以外の$Z$ や $\mu _ k, \sigma _ k$ らはすべて互いに独立に事前分布が設定されることになります。これはもちろん任意性があります。一方で $X$ を生成する分布である $p(X \mid Z, \mu _ 1, \sigma _ 1, \mu _ 2, \sigma _ 2, \mu _ 3, \sigma _ 3)$ は具体的にはどのような形で与えると良いでしょうか。もちろんこれも任意性がある話です。
今は、何らかの単純な分布が幾つも集まってできているような分布を考えています。そして、観測される $X$ は幾つかの単純な分布の中から1つが選ばれ、そして、その分布から確率的に値が生成されるという状況になります。 そのようなことをどうやって表現すれば良いのかが問題です。
さて、当然天下り的に確率変数を準備したのですから、これらの変数を使えば上記の生成の様子を表現する分布を作ることができます。キーポイントは $Z$ が、幾つかの単純な分布から1つを選択するような変数として扱うことです。
今、3つの単純な分布 $p(X \mid \mu _ 1, \sigma _ 1), p(X \mid \mu _ 2, \sigma _ 2), p(X \mid \mu _ 3, \sigma _ 3)$ を考えます。仮に下記のような掛け算を実施したとしましょう。
$$ p(X \mid \mu _ 1, \sigma _ 1)p(X \mid \mu _ 2, \sigma _ 2)p(X \mid \mu _ 3, \sigma _ 3) $$
これは我々が作りたいと思っている分布でしょうか？違いますね。我々が作りたいのは、$Z$ によって1つの分布が選ばれるようなモデルです。これを手続き的な動作で書くならば
$$ \begin{align} p(X \mid \mu _ 1, \sigma _ 1) \ \ (if \ Z = 1)\\ p(X \mid \mu _ 2, \sigma _ 2) \ \ (if \ Z = 2)\\ p(X \mid \mu _ 3, \sigma _ 3) \ \ (if \ Z = 3) \end{align} $$
となっていれば嬉しいのです。これを数式で書ければ嬉しいのです。これを成り立たせるために上の式を
$$ \begin{align} p(X \mid \mu _ 1, \sigma _ 1) ^ 1 p(X \mid \mu _ 2, \sigma _ 2) ^ 0 p(X \mid \mu _ 3, \sigma _ 3) ^ 0 \ \ (if \ Z = 1)\\ p(X \mid \mu _ 1, \sigma _ 1) ^ 0 p(X \mid \mu _ 2, \sigma _ 2) ^ 1 p(X \mid \mu _ 3, \sigma _ 3) ^ 0 \ \ (if \ Z = 2)\\ p(X \mid \mu _ 1, \sigma _ 1) ^ 0 p(X \mid \mu _ 2, \sigma _ 2) ^ 0 p(X \mid \mu _ 3, \sigma _ 3) ^ 1 \ \ (if \ Z = 3) \end{align} $$
と書いてみます。これは、$1$ 乗されている因子だけが残って、$0$ 乗されている因子は $1$ になるので、先ほど書いた if 文の式と全く等価なものになっています。そして、この$1$ 乗したり $0$ 乗したりしている部分を $Z$ が担ってくれれば、if 文での表現をほどくことができます。そのような表現を得るためには $Z = (Z1, Z2, Z3)$ という確率変数ベクトルであると考えなおして、成分のいずれか1つだけが $1$ の値となり、他の値は $0$ になるような確率変数だと考えればよいです（これは、one-hot表現と呼ばれる方法です）。このように確率変数ベクトル $Z$ を準備しておけば
$$ p(X \mid \mu _ 1, \sigma _ 1) ^ {Z _ 1} p(X \mid \mu _ 2, \sigma _ 2) ^ {Z _ 2} p(X \mid \mu _ 3, \sigma _ 3) ^ {Z _ 3} = \prod _ {k = 1} ^ 3 p(X \mid \mu _ k, \sigma _ k) ^ {Z _ k} $$
と書くことで、先程までの if 文表現を1つの数式で表現できるようになりました。またこのようにいずれかのインデックスを選択するような振る舞いをする $Z$ はカテゴリカル分布に従う事前分布を設定することで準備ができます。こうすると $Z$ は3つの成分を持つベクトルになっているわけですし、 $\mu, \sigma$ もそれぞれ3つずつ準備されているのですから、$\mu = (\mu _ 1, \mu _ 2, \mu _ 3)$ と $\sigma = (\sigma _ 1, \sigma _2, \sigma _ 3)$ という形で格納しておくと記述がすっきりします（ベクトルと言いつつ、要するに別々の変数として一々書くのが面倒なので、配列に一緒に入れておきます、くらいの雰囲気である）。
こうして、$X$ を生成する分布は
$$ p(X\mid Z, \mu , \sigma) = \prod _ {k = 1} ^ 3 p(X \mid \mu _ k, \sigma _ k) ^ {Z _ k} $$
と数式として具体的にモデル化できたことになり、同時分布全体は最初の単純作業で書き下した式に代入して
$$ p(X, Z, \mu , \sigma) = p(\mu _ 1)p(\sigma _ 1)p(\mu _ 2)p(\sigma _ 2)p(\mu _ 3)p(\sigma _ 3)p(Z) \prod _ {k = 1} ^ 3 p(X \mid \mu _ k, \sigma _ k) ^ {Z _ k} $$
と表現ができました。これを踏まえて、グラフィカルモデルの方も簡略表記にしておきましょう。
ガウス混合モデル
上記の式はガウス混合モデルを念頭において式を作ってきました。より一般にK個の正規分布が集まってできている分布を想定すれば、
$$ p(X, Z, \mu , \sigma) = p(\mu)p(\sigma)p(Z) \prod _ {k = 1} ^ K \mathcal N(X \mid \mu _ k, \sigma _ k) ^ {Z _ k} $$
とガウス混合モデルを書くことができます。ここで $p(\mu) = p(\mu _ 1)p(\mu _ 2)...p(\mu _ K)$ を略記しているだけであり、$\sigma$ の方も同様です（同時分布の設計を最初からやり直すことを考えれば $p(\mu)$ という複数個の平均パラメータを格納したベクトルの事前分布を考えるとき、必ずしもすべての分布の平均が互いに独立していると仮定しなくても良い！なぜなら、混合モデルの幾つかの分布は互いに似通っており、平均パラメータには何らかの相関があることも考えられるからである。それは分析者がデータを見てどのようにモデリングをするのかを考えるべきところである。今回の記事では、最初にグラフィカルモデルを書いた時点で、各平均パラメータや分散パラメータは互いに線でつながっておらず、独立であると"仮定"しただけの話である）。
ポアソン混合モデル
単純な分布の集まりが混合モデルなのであったから、単純な分布とやらがガウス分布である道理はありません。 データ $X$ が非負整数値になっている、すなわちカウントデータであるのならばポアソン分布を単純な分布として選択しても良いのです。その場合、ポアソン分布は平均パラメータ $\lambda$ のみによって分布形状が定まり、その分布で $X$ を生成するのですから
$$ p(X, Z, \lambda) = p(\lambda)p(Z) \prod _ {k = 1} ^ K {\rm Poi}(X \mid \lambda _ k) ^ {Z _ k} $$
という表記の仕方になります。当然 $\lambda = (\lambda _ 1, ..., \lambda _ K)$ と$K$ 種類のパラメータを格納しているという設定になります。 これは単にデータの生成の部分を書き換えているだけですし、事前分布も必要なパラメータに敷いているだけの話であり、混合モデルの考え方は共通です。 混合モデルは非常に強力で、単純な分布を集めることで峰が複数あるような分布を表現できるので、必ず抑えておきたい基本的な手法となります。
混合モデルの事後分布
さて、単純な分布を今、具体的な分布で与えずに混合モデルを考えていきましょう。単純な分布、すなわち $X$ を生成する分布のパラメータを $\theta$ と仮においてしまうことにします。すなわち同時分布は
$$ p(X, Z, \theta) = p(\theta)p(Z) \prod _ {k = 1} ^ K p(X \mid \theta _ k) ^ {Z _ k} $$
となります。ここまでは単一の $X$ が生成される様子を記述してきましたが、推論を開始する際には複数個、$N$個のデータ $X ^ N$ が集まっている状況を想定します。その場合には同時分布はどのようになるでしょうか。 ここは注意が必要です。まずある一つの $X$ に対して、どの分布が選ばれるのかを決める確率変数 $Z$ は個々に対応づいているはずです。もしもデータ $X$ に対して同じ数だけこの確率変数 $Z$ を準備しなかった場合、データ $X$ は単純な分布の集まりの中から、1回だけ1つの分布を選んで、ずっとその分布から生成されることになってしまいます。毎回、どの分布から生成するかを選びなおしているのですから、データの個数と $Z = (Z _ 1, ..., Z _ K)$ の個数は一緒でなければならないのです。
今、$Z$ は確率変数ベクトルであり、下添字はその成分であるとしているので、表記が厄介になりますが、データの個数に対するインデックスは上添字で $X ^ {(i)}$ とつけることにします。それを踏まえて同時分布は
$$ p(X ^ N, Z, \theta) = p(\theta)p(Z) \prod _ {i = 1} ^ N \left [ \prod _ {k = 1} ^ K p(X ^ {(i)} \mid \theta _ k) ^ {Z _ k ^ {(i)}} \right] $$
と書かれることになります。怖気づかないでください。今まで単一のデータを生成するような同時分布を書いてきた中で、データを生成している分布 $ \prod _ {k = 1} ^ K p(X \mid \theta _ k) ^ {Z _ k}$ をデータの数だけ繰り返しているということにすぎません。
実はこれ、私は混合モデルの最終的な式だけを始めてみた時、$\prod _ {i = 1} ^ N \prod _ {k = 1} ^ K p(X _ i \mid \theta _ k) ^ {Z _ {ik}} $ みたいな感じの表記であったため、非常に混乱したのです。あれ、$Z$ って行列だったんだっけ？？？などとしょうもない疑問を抱いたり、苦労した記憶があるので、回りくどく単独のデータ生成から入り、このように添字やカッコを付けました。また、通常の混合モデルの説明では $Z$ をいきなり導入することはせず、単に単純な分布の線形和 $\sigma _ k \pi _ k p (X \mid \theta _ k)$ から入ることが多いです。これは潜在変数 $Z$ というものを考えるのが特殊なテクニックであると考えている場合での導入方法で、EMアルゴリズムのお膳立てとして出現してくるという流れになります。これは、複雑な階層モデリングをする上で潜在変数を導入することが当たり前な古今では遠回りな導入だと感じています。また、EMアルゴリズムは勉強という観点では非常に重要なアルゴリズムなのですが、事実上、EMアルゴリズムは、複雑なモデルではそれほど良い推測にならない最尤推定の逐次推定アルゴリズムであり、よほど性質の良い問題でない限り（それこそ混合ガウスモデル）、積極的には使われません。複雑なモデルでは事前分布を導入した同時分布の設計を行い、階層ベイズモデリングを実施する意義が十分にあるため、その場合は変分ベイズ法を利用することで殆ど一切EMアルゴリズムは顔を出しません。
さて、事後分布はいつもどおり、乗法定理（ベイズの定理）を活用すれば、同時分布を観測変数 $X ^ N$ の分布で割れば良いだけであり、上記までで仮定した同時分布を代入して
$$ p(\theta, Z\mid X ^ N) = \frac{p(\theta)p(Z) \prod _ {i = 1} ^ N \left [ \prod _ {k = 1} ^ K p(X ^ {(i)} \mid \theta _ k) ^ {Z _ k ^ {(i)}} \right]}{p(X ^ N)} $$
となります。一応、複数個の観測を伴ったグラフィカルモデルも最後に乗せておきましょう。
$\alpha, \gamma$ はそれぞれ事前分布に対する超パラメータです。どんな事前分布を、どんなハイパーパラメータで設定するのかは分析者に任せます。
更に階層を増やす
さて、ここまでを見てこのモデルに疑問はあるでしょうか。例えば、ある資料ではガウス混合モデルの個々のガウス分布の平均パラメータは個々に事前分布を設定するが、分散パラメータは共通の事前分布を設定するというような場合もあるようです（それはデータに応じて柔軟に設計すれば良い）。そういった意味では、今回のモデルのほうが、個々の分布に多様性を許していることになります（その分推論は大変だ）。
もう1つは $Z$ に対するハイパーパラメータの設定です。既に述べたとおり、$Z$ はどの分布を選択するかを決める変数であり、これはカテゴリ分布から生成されることになります。今、このカテゴリ分布を事前分布として設定しているわけですが、これには1つ問題があるようです。なぜなら、仮にカテゴリ分布の超パラメータとして、すべてのインデックスが当確率で選ばれるだとか、$2:3:5$ で選ばれるだとかを設定したならば、推論の中でその値は固定値として用いられます。つまり、例えば当確率で超パラメータが設定されているのであれば、それを前提に、推論が行われ、$X _ i$ に対して $Z _ i$ がどうであるのかが割り当てられることになります。もしも $X _ i$ の多くがどれかの分布に偏っているのであれば、そのように超パラメータを設定すべきですし、できることならばどの分布が選ばれやすいのか、その割り振られる割合（つまりいま設定した超パラメータ自体）も推論の中に入れてやりたくなるわけです。
そういった場合には、超パラメータも勝手に設定する固定値ではなく、事前分布を設定してやり、超パラメータではなくしてあげるということができます（これを超事前分布を設定すると表現するが、そんなこと言い出したら、超事前分布の超超パラメータが表れ、更に超超事前分布を設定したくなり…となってしまいます。というわけで、どこまでを推論の枠組みに入れたいかを考え、推論の枠組みに入れたいものはすべて確率変数として扱い、それを決め終えたあとに、尚、固定値として与えなければならないものを超パラメータと呼べばいいと思います）。
このような方法を階層モデリングなどと言うわけですが、とにかく、$Z$ の超パラメータに相当する図の $\alpha$ を確率変数とみなしてやることにしてしまえば良いのです。グラフィカルモデルでは
と表記されることになります。この時新たに現れた $\eta$ は $p(\alpha)$ の超パラメータであります。$Z$ はカテゴリ分布から生成されるので、$\alpha$ はカテゴリ分布のパラメータになっているべきで、すなわち区間 $[0, 1)$の中の値を取るスカラーが$K$個格納されているベクトル $\alpha = (\alpha _ 1,..., \alpha _ K)$ となっているべきでありかつ、$\sum _ k \alpha _ k= 1$ という要請が入ります。それを満たすようなベクトルを生成する分布としてはディリクレ分布が知られているので、これを $\alpha$ の事前分布として構えてあげて、ディリクレ分布のパラメータを超パラメータ $\eta$ として準備するという形になるでしょう。
今、同時分布は
$$ p(X, Z, \theta, \alpha) = p(\alpha)p(\theta)p(Z \mid \alpha) \prod _ {k = 1} ^ K p(X \mid \theta _ k) ^ {Z _ k} $$
と修正されたことになります。ここにハイパーパラメータの存在を明記すると、
$$ p(X, Z, \theta, \alpha \mid \gamma, \eta) = p(\alpha\mid \eta)p(\theta\mid \gamma)p(Z \mid \alpha) \prod _ {k = 1} ^ K p(X \mid \theta _ k) ^ {Z _ k} $$
となります。
最後に
本当は時系列モデルを書こうと思ったのですが、隠れマルコフモデルなどに繋げる上では混合モデルの知識が必須となってくるため、先に混合モデルを説明しました。 次回こそ時系列＋状態空間モデルの説明に入るかもしれません…。
