読者です 読者をやめる 読者になる 読者になる

人工知能に関する断創録

人工知能、認知科学、心理学、ロボティクス、生物学などに興味を持っています。このブログでは人工知能のさまざまな分野について調査したことをまとめています。最近は、機械学習・Deep Learningに関する記事が多いです。



フィッシャーの線形判別

機械学習 PRML

今回は、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()

結果は、

f:id:aidiary:20100425100242p:plain

おおー、データの平均だけではなく、共分散も考慮するとちゃんと分類できました。フィッシャーの線形判別は2クラスだけではなく、多クラスにも拡張できるようです(4.1.6)。今回は、省略します。