フィッシャーの線形判別 - 人工知能に関する断創録

今回は、4.1.4のフィッシャーの線形判別を試してみました。これは、他の手法と少し毛色が違う感じがします。まず、D次元の入力ベクトルxを（4.20）で1次元ベクトル（スカラー）に射影します。ベクトル同士の内積なので結果はスカラーで、wはxを射影する方向を表します。

$(4.20) \hspace{50pt} y = \bf{w}^T \bf{x}$

フィッシャーの線形判別は、射影後のデータの分離度をもっとも大きくするようなデータの射影方向wを見つけるという手法だそうです。

クラス1のデータ集合C1の平均ベクトルとクラス2のデータ集合C2の平均ベクトル（4.21）をw上へ射影したクラス間平均の分離度（4.22）を最大にするwを選択するのが1つめのポイントのようです。式（4.22）の左辺はスカラーです（フォントの違いがわかりにくい）。

$(4.21) \hspace{50pt} \displaystyle \bf{m_1} = \frac{1}{N_1} \sum_{n \in C1} \bf{x_n}, \hspace{50pt} \bf{m_2} = \frac{1}{N_2} \sum_{n \in C2} \bf{x_n}$
$(4.22) \hspace{50pt} \displaystyle m_2 - m_1 = \bf{w}^T (\bf{m_2} - \bf{m_1})$

wは単位長であるという制約のもとで（4.22）を最大化するようにラグランジュ未定乗数法で解くと、

$\bf{w} \propto (\bf{m_2} - \bf{m_1})$

という解が得られます（演習4.4）。これは、ベクトルwとベクトル(m2-m1)が平行ってことを言っているだけです。だけど、ここで終わると図4.6にあるような決定境界は描画できません。これは、図4.6からの想像なのですが、決定境界の直線は、wに直交し、クラス平均同士を結ぶ直線の中点を通る直線だと仮定してます。そこで、下図のようにして決定境界の直線の方程式を求めました。PRMLの原書サポートページで配布されている図表を加筆して引用します。

f:id:aidiary:20100425093410p:plain

まず、wと直交するベクトルw'を求めます。これは、直交するベクトル同士の内積が0になることを利用して、

$w = (w_1, w_2)^T$
$w' = (- w_2, w_1)^T$

のように求められます。なので、決定境界の直線の方程式y=ax+bの傾きはw'からyの増加量/xの増加量で

$\displaystyle a = - \frac{w_1}{w2}$

とわかります。傾きがaで点 m = (m2-m1)/2（クラスの平均を結んだ線の中点）を通るので代入するとbも求まります。では、ここまでのところをプログラミングしてみます。

#coding:utf-8
# 4.1.4 フィッシャーの線形判別（p.185）

import numpy as np
from pylab import *
import sys

N = 100  # データ数

def f(x, a, b):
    # 決定境界の直線の方程式
    return a * x + b

if __name__ == "__main__":
    # 訓練データを作成
    cls1 = []
    cls2 = []

    # データは正規分布に従って生成
    mean1 = [1, 3]  # クラス1の平均
    mean2 = [3, 1]  # クラス2の平均
    cov = [[2.0,0.0], [0.0, 0.1]]  # 共分散行列（全クラス共通）
    
    # データ作成
    cls1.extend(np.random.multivariate_normal(mean1, cov, N/2))
    cls2.extend(np.random.multivariate_normal(mean2, cov, N/2))
    
    # 各クラスの平均をプロット
    m1 = np.mean(cls1, axis=0)
    m2 = np.mean(cls2, axis=0)
    plot([m1[0]], [m1[1]], 'b+')
    plot([m2[0]], [m2[1]], 'r+')
    print m1, m2
    
    # 訓練データを描画
    x1, x2 = np.array(cls1).transpose()
    plot(x1, x2, 'bo')
    
    x1, x2 = np.array(cls2).transpose()
    plot(x1, x2, 'ro')
    
    # ラグランジュの未定乗数法で解いた結果
    w = m2 - m1
    
    # 識別境界を描画
    # wは識別境界と直交するベクトル
    a = - (w[0] / w[1])  # 識別直線の傾き
    
    # 傾きがaで平均の中点mを通る直線のy切片bを求める
    m = (m1 + m2) / 2
    b = -a * m[0] + m[1]  # 識別直線のy切片
    
    x1 = np.linspace(-2, 6, 1000)
    x2 = [f(x, a, b) for x in x1]
    plot(x1, x2, 'g-')

    xlim(-2, 6)
    ylim(-2, 4)
    show()

結果は、下のようになり、ちゃんと分類できてません。

f:id:aidiary:20100425094631p:plain

これは、クラス分布の非対角な共分散が強いために起きる問題だそうです。フィッシャーの線形判別は、射影されたクラス平均間の分離度を大きくするだけではなく、各クラス内では小さな分散を与える関数を最大化するというもう1つの条件が追加されます。

クラスCkから射影されたデータのクラス内分散は、

$(4.24) \hspace{50pt} \displaystyle s_k^2 = \sum_{n \in C_k} (y_n - m_k)^2$

なので、フィッシャーの判別基準は、

$(4.25) \hspace{50pt} \displaystyle J(\bf{w}) = \frac{(m_2 - m_1)^2}{s_1^2 + s_2^2} = \frac{\bf{w}^T \bf{S_B} \bf{w}}{\bf{w}^T \bf{S_W} \bf{w}}$
$(4.27) \hspace{50pt} \displaystyle \bf{S_B} = (\bf{m_2} - \bf{m_1})(\bf{m_2} - \bf{m_1})^T$
$(4.28) \hspace{50pt} \displaystyle \bf{S_W} = \sum_{n \in C_1} (\bf{x_n} - \bf{m_1})(\bf{x_n} - \bf{m_1})^T + \sum_{n \in C_2} (\bf{x_n} - \bf{m\2})(\bf{x_n} - \bf{m\2})^T$

で定義されます。S_Bはクラス間共分散行列、S_Wはクラス内共分散行列と呼ぶようです。ここでは、(4.25)が最大になるようなwを求めろってことですね。分子は大きいほどよくて、分母は小さいほどよい、なるほど。これを解くと、

$(4.30) \hspace{50pt} \displaystyle \bf{w} \propto \bf{S_W}^{-1} (\bf{m_2} - \bf{m_1})$

となります。おー、さっきと少し違ってクラス内共分散S_Wの逆行列がくっついた結果になります。これがあるとどうよくなるのかさっそく確かめてみます。

#coding:utf-8
# 4.1.4 フィッシャーの線形判別（p.185）

import numpy as np
from pylab import *
import sys

N = 100  # データ数

def f(x, a, b):
    return a * x + b

if __name__ == "__main__":
    # 訓練データを作成
    cls1 = []
    cls2 = []
    
    # データは正規分布に従って生成
    mean1 = [1, 3]  # クラス1の平均
    mean2 = [3, 1]  # クラス2の平均
    cov = [[2.0,0.0], [0.0, 0.1]]  # 共分散行列（全クラス共通）
    
    # データ作成
    cls1.extend(np.random.multivariate_normal(mean1, cov, N/2))
    cls2.extend(np.random.multivariate_normal(mean2, cov, N/2))
    
    # 各クラスの平均をプロット
    m1 = np.mean(cls1, axis=0)
    m2 = np.mean(cls2, axis=0)
    plot([m1[0]], [m1[1]], 'b+')
    plot([m2[0]], [m2[1]], 'r+')
    print m1, m2
    
    # 総クラス内共分散行列を計算
    Sw = zeros((2, 2))
    for n in range(len(cls1)):
        xn = matrix(cls1[n]).reshape(2, 1)
        m1 = matrix(m1).reshape(2, 1)
        Sw += (xn - m1) * transpose(xn - m1)
    for n in range(len(cls2)):
        xn = matrix(cls2[n]).reshape(2, 1)
        m2 = matrix(m2).reshape(2, 1)
        Sw += (xn - m2) * transpose(xn - m2)
    Sw_inv = np.linalg.inv(Sw)
    w = Sw_inv * (m2 - m1)
    
    # 訓練データを描画
    x1, x2 = np.transpose(np.array(cls1))
    plot(x1, x2, 'bo')
    
    x1, x2 = np.transpose(np.array(cls2))
    plot(x1, x2, 'ro')
    
    # 識別境界を描画
    # wは識別境界と直交するベクトル
    a = - (w[0,0] / w[1,0])  # 識別直線の傾き
    
    # 傾きがaでmを通る直線のy切片bを求める
    m = (m1 + m2) / 2
    b = -a * m[0,0] + m[1,0]  # 識別直線のy切片
    
    x1 = np.linspace(-2, 6, 1000)
    x2 = [f(x, a, b) for x in x1]
    plot(x1, x2, 'g-')
    
    xlim(-2, 6)
    ylim(-2, 4)
    show()