「機械学習理論入門」第7章 EMアルゴリズム：最尤推定法による教師なし学習

今回は教師無し学習のクラスタリングアルゴリズムとして最尤推定法を利用したEMアルゴリズムを行う

教師あり学習と教師なし学習の違い

学習	特長
教師あり学習	ｌすでに目的変数がわかっているデータを分析することで未知のデータに対して目的変数の値を推測するルールを導き出す
教師なし学習	与えられたデータを明示的に分類する目的変数がわからない状態でデータ間の類似性を発見していく

具体的には手書き文字の分類問題を取り扱う

今回やること

代表文字の作成
数字の分類
サンプルコードの解説

代表文字の作成

まずがぞうデータの分類は後回しにして、まず特定の数字の手書き画像データについて
それらを平均化した「代表文字」を作成する
例）複数の顔写真を合成して「平均的な顔」の写真を作る

合成方法

「全ての画像を重ね合わせて平均をとる」
手書き文字画像のピクセル一つ一つが黒いか白いかで0,1に振り分けそれをベクトルにする
つまりベクトルのn要素目がnピクセル目が黒いか白いかに対応している

ベクトル化を全ての画像について行いそれら全てを平均する

    μ = \frac{1}{N}\sum_{n=1}^{N}X_n \tag{7.1}

すると全てのベクトルが0-1の値をとる
これをそのピクセルの濃淡とする
そうすると図7.2のようにそれっぽい結果となる

しかし、「なぜこれでうまくいくのか」という事は説明できない
よってそれが説明できるモデルを作成する必要がある

画像生成器による最尤推定法の適用

画像生成器

まず、ランダムに画像を生成する「画像生成器」を用意する

これは図7.2のような濃淡のある画像ファイル
このピクセル1つ1つの濃淡がそのままそのピクセルが黒になる確率と考える

つまり濃淡のある画像ベクトルμの第i成分をμiとして
「確率μiでi番目のピクセルを黒にする」というルールでランダムに画像を生成する
そうすると、図7.3の用に手書き文字を生成できる

最尤推定法の適用

トレーニングセットとなる数字の手書きか画像がN個あるとする
次にその、画像生成器を用意してN個の画像を生成する
そのとき「トレーニングセットと全く同じデータ群が生成される確率」を考える

※その確率は非常に低いものになるが、とにかく何らかの確率は計算される
この確率が最大になるような画像生成器を見つけ出せば、トレーニングセット
を代表する画像になっていると期待できる

先ほどのベクトルにおいて、i番目のピクセルに注目した場合、
そのピクセルが黒い確率pi（xi = 1）は

    P_i = μ_i \tag{7.2}

そのピクセルが白い確率pi（xi = 1）は
math P_i = 1 - μ_i \tag{7.3}

黒い確率と白い確率をまとめて書くと

    P_i = μ_i^{x_i}(1 - μ_i)^{1-μ_i} \tag{7.4}

黒と白のパターンがすべてのピクセルで同じ場合の確率

    p(x) = \prod_{i=1}^{D}p_i = \prod_{i = 1}^{D}μ_i^{x_i}(1-μ_i)^{1-x_i} \tag{7.5}

さらにトレーニングセットに含まれる全てのデータを考えると、その全てに合致する確率は

    p(x) = \prod_{i=1}^{D}p_i(x_n) = \prod_{n = 1}^{N}\prod_{i = 1}^{D}μ_i^{[x_n]_i}(1-μ_i)^{1-[x_n]_i} \tag{7.5}

これがこのモデルにおける尤度関数
これを最大にするμを計算すると7.1が得られる

今回の場合7.4が数学的にはベルヌーイ分布と呼ばれる確率分布になるため
「ベルヌーイ分布」を用いた最尤推定法となる

混合分布を用いた最尤推定法

さっきは特定の数字の手書き画像データから代表文字を作成した
次に、複数の数字の手書き画像データが混在している場合において、文字の種類毎に画像データを分類する

今回の場合、それぞれの数字に対応する画像生成器を用意する必要があるため、K個作成する

使用する画像生成器の選択にも確率を取り入れる
つまり「どれか1つの画像生成器をランダムに選択して、新しい画像を生成する」
この時、k番目の画像生成器を選ぶ確率をπkとする

このπkは以下の条件を満たす

\sum_{k=1}^Kπ_k = 1

このような操作によって特定の画像Xが得られる確率は次式で表される

p(x) = \sum_{k=1}^{K}π_kp_{nk}(x) \tag{7.12}

で全てのデータに対して完璧なパターンを返す確率は

P = \prod_{n=1}^{N}p(x_n) = \prod_{n=1}^{N}\sum_{k=1}^{K}π_kp_{μ_k}(x_n) \tag{7.13}

これがこのモデルの尤度関数
これを最大化する

EMアルゴリズム

こんかいの場合はいつものように対数尤度をとっても最大値が求まらない
EMアルゴリズムはこんな形式の尤度関数を最大にするパラメータを求める手続き

どれか1つの画像生成器を確率に従って選択して、新しい画像を生成して、画像xnが得られる確率は

p(x_n) = \sum_{k=1}^{K}π_kp_{nk}(x_n) \tag{7.14}

これはk個ある画像生成器全てにおいてxnができる確率を足し合わせたもの
そこでk番目の画像生成器から画像xnが得られる確率を計算すると

r_{nk} = \frac{π_kp_{nk}(x_n)}{\sum_{k=1}^{K}π_kp_{nk}(x_n)} \tag{7.15}

k平均法の場合には、もっとも距離の近い代表点に所属するという条件で、どの代表点に所属するかを示す変数r_{nk}の値を設定した。
今回のr_{nk}もそれに該当する。

その画像が、どの画像生成器に所属する割合が決まったら
その割合に基づいて新しく画像生成器を作成する。

μ_k = \frac{\sum_{n=1}^{N}r_{nk}X_n}{\sum_{n=1}^{N}r_{nk}} \tag{7.16}

これは7.1を混合文字の場合に改造したもの
7.1とは異なり確率の平均をとっている

また、その画像生成を選択する確率は

π_k = \frac{\sum_{n=1}^{N}r_{nk}}{N} \tag{7.17}

7.17はそれぞれの画像生成器を使用する割合を「それぞれの画像生成器に所属する画像の量」に比例するように再設定している

この後は計算しなおした値を用いて7.15を計算
そしてさらにまた7.16,17を計算する
というのを繰り返すと最終的に極大値になることが証明されている

図7.6を参照

サンプルコード

実際にEMアルゴリズムの手続きを実行してみる

# -*- coding: utf-8 -*-
#
# 混合ベルヌーイ分布による手書き文字分類
#
# 2015/04/24 ver1.0
#

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame
from numpy.random import randint, rand

#------------#
# Parameters #
#------------#
K = 3   # 分類する文字数
N = 10  # 反復回数


# 分類結果の表示
def show_figure(mu, cls):
    fig = plt.figure()
    for c in range(K):
        subplot = fig.add_subplot(K,7,c*7+1)
        subplot.set_xticks([])
        subplot.set_yticks([])
        subplot.set_title('Master')
        subplot.imshow(mu[c].reshape(28,28), cmap=plt.cm.gray_r)
        i = 1
        for j in range(len(cls)):
            if cls[j] == c:
                subplot = fig.add_subplot(K,7,c*7+i+1)
                subplot.set_xticks([])
                subplot.set_yticks([])
                subplot.imshow(df.ix[j].reshape(28,28), cmap=plt.cm.gray_r)
                i += 1
                if i > 6:
                    break
    fig.show()

# ベルヌーイ分布(7.10)
def bern(x, mu):
    r = 1.0
    for x_i, mu_i in zip(x, mu):
        if x_i == 1:
            r *= mu_i 
        else:
            r *= (1.0 - mu_i)
    return r

# Main
if __name__ == '__main__':
    # トレーニングセットの読み込み
    df = pd.read_csv('sample-images.txt', sep=",", header=None)
    data_num = len(df)

    # 初期パラメータの設定
    mix = [1.0/K] * K #πk その画像生成器が選ばれる確率
    mu = (rand(28*28*K)*0.5+0.25).reshape(K, 28*28)#μ適当に作成した画像
    for k in range(K):
        mu[k] /= mu[k].sum()

    fig = plt.figure()
    for k in range(K):
        subplot = fig.add_subplot(K, N+1, k*(N+1)+1)
        subplot.set_xticks([])
        subplot.set_yticks([])
        subplot.imshow(mu[k].reshape(28,28), cmap=plt.cm.gray_r)
    fig.show()

    # N回のIterationを実施
    for iter_num in range(N):
        print "iter_num %d" % iter_num

        # E phase　(7.15)rnkを計算する
        resp = DataFrame()
        for index, line in df.iterrows():#画像の枚数だけループ
            tmp = []
            for k in range(K):
                a = mix[k] * bern(line, mu[k])#(7.15)rnkの分子
                if a == 0:
                    tmp.append(0.0)
                else:
                    s = 0.0 #(7.15)rnkの分母
                    for kk in range(K):
                        s += mix[kk] * bern(line, mu[kk])
                    tmp.append(a/s)
            resp = resp.append([tmp], ignore_index=True)#(7.15)rnk

        # M phase　(7.16)(7.17)を計算する部分
        mu = np.zeros((K, 28*28))
        for k in range(K):
            nk = resp[k].sum()
            mix[k] = nk/data_num #(7.17)
            for index, line in df.iterrows():
                mu[k] += line * resp[k][index]
            mu[k] /= nk#(7.16)

            #画像を生成する部分
            subplot = fig.add_subplot(K, N+1, k*(N+1)+(iter_num+1)+1)
            subplot.set_xticks([])
            subplot.set_yticks([])
            subplot.imshow(mu[k].reshape(28,28), cmap=plt.cm.gray_r)
        fig.show()

    # トレーニングセットの文字を分類
    cls = []
    for index, line in resp.iterrows():
        cls.append(np.argmax(line[0:]))

    # 分類結果の表示
    show_figure(mu, cls)

図7.7
全ての画像データが正しく分類されているわけではない

図7.10 図7.11
クラスターを増やしてみるとこんな感じ
結果の解釈によっては恣意的になる可能性も…

クラスタリングをすつ際には、クラスター数を変えながら何度か実行してトレーニングセットの特徴を探索していく必要があります。
これは、教師無し学習の特徴である
さらに、クラスタリングで得られた結果をどう判断するかは利用者の主観が入る
このような点が理解できていないと特定の条件下で得られたクラスタリングを絶対的に正しいと思い込む可能性がある
機械学習のアルゴリズムは、それぞれの理論的背景や特性を理解した上で活用する事が大切である。

【読書会】機会学習理論入門_第7章