EMアルゴリズムまとめ - 気まぐれエンジニア日記

なんか気分なので, また記事を書こうと思います. 筆者は大学生なので, 定期試験と言うものを受けなければいけません. 今回の定期試験の範囲にEMアルゴリズムっていうのがあるのでそれについてまとめてみたいと思います. パターン認識を専門にしようと考えている以上理解しなければいけないので, 自らの頭の整理もかねて文章に起こそうと思います.

1. 確率モデルと最尤推定

EMアルゴリズムの説明に入る前に, 確率モデルあるいは生成モデルと最尤推定について説明しなければいけません.

確率モデルとは, 確率的に様々な値をとる事象があった時に, それぞれの値がどのくらい発生しやすいかを数式でモデル化したものです. 正規分布や指数分布, 一様分布, 二項分布といったような確率分布の名前を聞いたことがある人も多いと思います. 確率モデルといった場合には確率密度関数のことを指すことが多いですが, 累積分布関数の格好で書いても, どんな形であっても確率を数式で表したものは広く確率モデルと呼ばれます.

さて, 確率モデルにはいくつかパラメタが存在する場合があります. 例えば, 単に正規分布といっても様々なものがあって, 一般にはその確率密度関数は平均 $\mu$ , 分散 $\sigma ^2$ をパラメタとしてもった次の式で書かれます :
$\displaystyle N(x; \mu, \sigma ^2) := \frac{1}{\sqrt{2 \pi \sigma ^2}}e^ {-\frac{(x - \mu) ^2}{2\sigma ^2}}$
パターン認識では, あるモデルに従って発生した値の列がデータとして与えられた時に, そのパラメータを推定したい, なんてことがまま起こります. 集めたデータは正規分布に従って発生しているはずだということはわかるのだけど, その平均と分散を推定したい, などです. その推定の1つの方法として最尤推定があります.

パラメタとして $\theta$ を持つ確率モデル $M(\theta)$ から値 $x_1, x_2, \cdots, x_n$ が観測されたとしましょう. このとき, モデル $M(\theta)$ における, 観測データ $x_1, x_2, \cdots, x_n$ に対する尤度 $L(M(\theta))$ を次で定義します :
$\displaystyle L(M(\theta)) := \prod_{i = 1}^n P(x_i|M(\theta))$
ここに, $P(x_i|M(\theta))$ は, モデル $M(\theta)$ から値 $x_i$ が生起する確率です.

最尤推定とは, 尤度が最大となるようにモデルのパラメタ $\theta$ を推定することを言います. すなわち最尤推定によって求められた最適なモデルのパラメタ $\theta ^*$ は次の式でかけます :
$\theta ^* = \mathrm{argmax} L(M(\theta))$
今, 尤度を定義しましたが, この値はモデル $M(\theta)$ から値 $x_1, x_2, \cdots, x_n$ がこの順に出てくる確率に他ならないことに注意しましょう. すなわち, 最尤推定とは, モデル $M(\theta)$ から値 $x_1, x_2, \cdots, x_n$ がこの順に出てくるという, 今まさに観測した事象が最も起こりやすいようなモデルを最適としましょうといっているのです. とりあえずモデルを動かしてみたら値 $x_1, x_2, \cdots, x_n$ がこの順に出てきた. じゃあこれが最も起こりやすいのが正解のモデルなんじゃないのというスーパー安直な方法なのです.

「なんでこのモデルが最適だと思ったの？」「だって今目の前で起こったことが1番起こりやすいんだもん」「それはそう」という会話なわけです, ええ.

2. 正規分布と最尤推定

確率モデルが正規分布のように単純なものであれば, 最尤推定は簡単に行うことができるということを今からお見せしましょう. $N(x; \mu, \sigma ^2)$ から, 値 $x_1, x_2, \cdots, x_n$ がこの順に出てきたとき, 平均 $\mu$ , 分散 $\sigma ^2$ を推定します. こういう時には, パラメータでのグラディエントを求めて, それが零ベクトルになるような場所を探せばいいわけです. さて, 対数関数は単調増加ですから, 尤度 $L(M(\theta))$ を最大化することは, その対数をとった
$\displaystyle \log{L(M(\theta))} = \sum_{i = 1}^n \log{P(x_i|M(\theta))}$
を最大化することと同値です. このとき, 正規分布では,
$\displaystyle \log{L(\theta)} = \sum_{i = 1}^n \log{\frac{1}{\sqrt{2 \pi \sigma ^2}}e^ {-\frac{(x_i - \mu) ^2}{2\sigma ^2}}} = \sum_{i = 1}^n \log{\frac{1}{\sqrt{2 \pi}} - \frac{1}{2}\log{\sigma ^2} -\frac{(x_i - \mu) ^2}{2\sigma ^2}}$
なので, $\mu$ で微分して,
$\displaystyle \frac{1}{\sigma ^2}\sum_{i = 1}^n (x_i - \mu) = 0$
$\displaystyle \mu = \frac{1}{n}\sum_{i = 1}^n x_i$
$\sigma ^2$ で微分して,
$\displaystyle \sum_{i = 1}^n - \frac{1}{2 \sigma ^2} + \frac{(x_i - \mu) ^2}{2} \cdot \frac{1}{(\sigma ^2) ^2} = 0$
$\displaystyle \sigma ^2 =\frac{1}{n} \sum_{i = 1}^n (x_i - \mu) ^2$
とよく見慣れた式が出てきました. つまり, 正規分布の平均 $\mu$ , 分散 $\sigma ^2$ の最尤推定の結果は, 観測したデータの平均と分散そのものであるということがわかります.

3. EMアルゴリズム

単純な確率モデルの最尤推定であれば, 先のように観測データから得られる尤度の関数をパラメタで微分して極値を求めることで, 解析的に行うことができました. しかしながら, 確率モデルの中に, 観測できないデータがある場合にはそう簡単にはいきません. 例として確率モデルの重ね合わせを考えます.
$M_1(\theta_1), M_2(\theta_2), \cdots, M_k(\theta_k)$ を確率モデルとし, パラメタを推定したいモデルが
$\displaystyle M := \sum_{i = 1}^k w_iM_i(\theta_i)$
ただし,
$\displaystyle \sum_{i = 1}^k w_i = 1\land\forall i\in \mathbb{N}\cap[1, n], w_i\geq 0$
と与えられているとします. このとき, 値 $x_1, x_2, \cdots, x_n$ を観測した場合の最尤推定をしようと思うのですが, モデルが重み付き和で表されているので, 先のように対数をとってもうまくいきそうにありません. そこでこんなことを考えます.
モデル $M$ は確率 $w_i$ でモデル $M_i(\theta_i)$ を選んで, そこから値を取り出していると見なすことができるんだから, 値 $x_1, x_2, \cdots, x_n$ に加えて, どのモデルから出てきたかの情報が得られれば良いのではないか. もし, どのモデルが出てきたかがわかってしまえば, 各モデル $M_i(\theta_i)$ から出てきた値だけに注目して, 先と同じように微分して極値を求められそうではないかと考えるわけです.

そこで, あるモデル
$\displaystyle M := \sum_{i = 1}^k w_iM_i(\theta_i)$
をとりあえず1つ固定して, 値 $x_1, x_2, \cdots, x_n$ が $M_1(\theta_1), M_2(\theta_2), \cdots, M_k(\theta_k)$ のどれから出てきたと考えるべきかの確率
$P(M_j|x_i, M)$
を求めてみよう, その上でもっといいモデル
$\displaystyle M^* := \sum_{i = 1}^k w_i^* M_i(\theta_i^*)$
に少しずつ更新しようと考えるわけです. 実は, この確率がもとまってしまえば, 新しいモデル $M'$ の推定は簡単にできてしまうよ, というのがEMアルゴリズムのキモなのです. 詳しく見ていきましょう.

今, モデル $M$ を1つ固定します. 最尤推定で最大化したい値はもっといいモデル $M^*$ における, 観測データ $x_1, x_2, \cdots, x_n$ に対する尤度 $L(M^*)$ で,
$\displaystyle L(M^*) := \prod_{i = 1}^n P(x_i|M^*)$
あるいはその対数
$\displaystyle \log{L(M^*)} = \sum_{i = 1}^n \log{P(x_i|M^*)}$
です.
今, 各 $x_i$ について,
$\displaystyle \sum_{j = 1}^k P(M_j|x_i, M) = 1$
ですから,
$\displaystyle \log{L(M^*)} = \sum_{i = 1}^n \log{P(x_i|M^*)} = \sum_{i = 1}^n 1\cdot \log{P(x_i|M^*)} = \sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M) \log{P(x_i|M^*)}$
となります.
次がポイントです. こんな式変形をしようと思ったのはなぜだったか思い出してください. 値 $x_1, x_2, \cdots, x_n$ に加えて, $M_1(\theta_1), M_2(\theta_2), \cdots, M_k(\theta_k)$ のうちのどのモデルから出てきたかの情報が得られればいいのになぁと考えていたわけです. すなわち,
$P(x_i, M_j |M^*)$
たちだけの式にしてしまえば, $\log{L(M^*)}$ を最大化するようなパラメータ $\theta_1, \theta_2, \cdots, \theta_k$ たちは簡単に求められそうだという条件下で話をしていたわけです. この値が出てくるように, もう1回だけ変形をしてみましょう. 条件付き確率の定義によって,
$\displaystyle P(x_i|M^*) = \frac{P(x_i, M_j |M^*)}{P(M_j |x_i, M^*)}$
ですから,
$\displaystyle \log{L(M^*)} = \sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M) \log{ \frac{P(x_i, M_j |M^*)}{P(M_j |x_i, M^*)}} = \sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M) \left(\log{P(x_i, M_j |M^*)} - \log{P(M_j |x_i, M^*)}\right)$
ここで,
$\displaystyle Q(M, M^*) := \sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M)\log{P(x_i, M_j |M^*)}$
$\displaystyle H(M, M^*) := -\sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M)\log{P(M_j |x_i, M^*)}$
とおけば,
$\log{L(M^*)} = Q(M, M^*) + H(M, M^*)$
となります.
はあ, やったやったと安心してはいけません. 最初の目的をもう一度思い出しましょう. モデル $M$ を固定して, そこから得られる情報を使って, ちょっといいモデル $M^*$ を導きたいのでした. つまり,
$\log{L(M^*)} - \log{L(M)} > 0$
となるように $M^*$ を決めましょうということになるのです. しかもその差はなるべく大きい方が良さそうですよね.
$\log{L(M^*)} - \log{L(M)} = Q(M, M^*) - Q(M, M) + H(M, M^*) - H(M, M)$
ですが, $H(M, M^*) - H(M, M)$ はダイバージェンスなので, $M^*$ として, $M$ とは違うものを選べば自動的に $H(M, M^*) - H(M, M) > 0$ となります. この説明は後述します. 一方, $Q(M, M^*) - Q(M, M))$ は, $Q(M, M^*)$ を最大化することができれば, 最大化できます. すなわち,
$\displaystyle Q(M, M^*) := \sum_{i = 1}^n \sum_{j = 1}^k P(M_j|x_i, M)\log{P(x_i, M_j |M^*)}$
を最大化するように $M^*$ を選べば良いということがわかります. 嬉しいことに, $P(M_j|x_i, M)$ はすでに求めることができていますし, $P(x_i, M_j |M^*)$ の値の評価は簡単にできるという条件で話をしていましたから, この問題は無事に解けそうそうだということがわかります.

$M$ から $M^*$ へちょっといいモデルに更新ができました. あとはこれを繰り返していけば, どんどんいいモデルになりそうですね. これがEMアルゴリズムです.

EMアルゴリズムの手順をまとめると次のようになります.

初期モデルを決める.
モデル $M$ とモデル化したい確率的な現象の値の列 $x_1, x_2, \cdots, x_n$ から, $P(M_j|x_i, M)$ を計算する.
$P(M_j|x_i, M)$ をもとに, $Q(M, M^*)$ を最大化するようなちょっといいモデル $M^*$ を推定する.
$M$ を $M^*$ に置き換えて, 2に戻る.

3. ダイバージェンスは正

ちょっとちょっと, 終わった気になっていませんか？？まだ, 1つ残っていますよ.

$H(M, M^*) - H(M, M)$ はダイバージェンスなので, $M^*$ として, $M$ とは違うものを選べば自動的に $H(M, M^*) - H(M, M) > 0$ となります.

ほんまか？？証明しましょう.

$p_{ij} := P(M_j|x_i, M), q_{ij} :=P(M_j |x_i, M^*)$ とおいて, 次を証明すれば良さそうです.

命題ダイバージェンスは正
$p_{ij}, q_{ij}$ が
$\displaystyle \sum_i \sum_j p_{ij} = 1\land\forall i, j, p_{ij}\geq 0$
および
$\displaystyle \sum_i \sum_j q_{ij} = 1\land\forall i, j, q_{ij}\geq 0$
を満たすとする. このとき,
$\displaystyle - \sum_i \sum_j p_{ij}\log q_{ij} + \sum_i \sum_j p_{ij}\log p_{ij} \geq 0$
等号成立は $p_{ij}, q_{ij}$ がそれぞれ等しいときである.

証明:
$f(x) = -\log x + x - 1$ とおく.
$\displaystyle f'(x) =-\frac{1}{x} + 1 = 0$
を解くと, $x = 1$ .
しかも $0 < x < 1\Rightarrow f'(x) < 0$ かつ $x > 1\Rightarrow f'(x) > 0$ より,
$f(x)$ は, $x = 1$ で最小値0を取ることがわかる.
したがって $x > 0\Rightarrow -\log x \geq 1 - x$ であり, 等号成立は $x = 1$ とわかる. このとき,

$\displaystyle - \sum_i \sum_j p_{ij}\log q_{ij} + \sum_i \sum_j p_{ij}\log p_{ij} = - \sum_i \sum_j p_{ij}\log\frac{q_{ij}}{p_{ij}}$
$\displaystyle \geq \sum_i \sum_j p_{ij}(1 - \frac{q_{ij}}{p_{ij}}) = \sum_i \sum_j (p_{ij} - q_{ij}) = \sum_i \sum_j p_{ij} - \sum_i \sum_j q_{ij} = 1 - 1 = 0$
等号成立は, $\displaystyle \frac{q_{ij}}{p_{ij}} = 1$ ,
すなわち $p_{ij} = q_{ij}$ のときである.

おしまい！！お疲れした〜！