Pythonでベイズ識別 - 元理系院生の新入社員がPythonとJavaで色々頑張るブログ

この記事に書かれていること

アヤメデータのベイズ識別方法
Pythonによるベイズ識別の実装

ベイズ識別とは

少し前にベイズの定理と、確率密度関数の記事を書きました。
<a href="http://emoson.hateblo.jp/entry/2015/02/19/225316">ベイズ理論 - 理系大学生がPythonで色々頑張るブログ</a>
ベイズ理論 - 理系大学生がPythonで色々頑張るブログ
<a href="http://emoson.hateblo.jp/entry/2015/02/20/154114">ベイズ理論(2) - 理系大学生がPythonで色々頑張るブログ</a>
ベイズ理論(2) - 理系大学生がPythonで色々頑張るブログ

今回はこのベイズの定理を識別問題に適用した、ベイズ識別器をPythonで実装します。

ベイズの定理のおさらい

ベイズの定理は2つの確率事象の同時確率が条件付き確率から求まる事を利用して、導出することが出来ました。
$P(A,B) = P(B|A)P(A) = P(A|B)P(B)$

$P(B|A)P(A) = P(A|B)P(B)$

$P(B|A) = \frac{P(A|B)P(B)}{P(A)}$

この式を次の様に考えると、識別問題に適用することが出来ます。

$P(B|A)$ ・・・データAが観測された時に、そのデータの分類クラスがBである確率

$P(A)$ ・・・データAが観測される確率

$P(A|B)$ ・・・分類クラスBにおけるデータAの確率分布
(言い回しが面倒なので厳密ではないですが、BにおけるデータAの生成確率と捉えておけば大丈夫だと思います）

$P(B)$ ・・・分類クラスBの生起確率

AとBというシンボルだと少し分かりにくいので、観測データを $\bf x$ 、分類クラスを $C_{i} (i = 1,...,N)$ とし、書き直すと次のようになります。

$P(C_{i}|\bf x) = \frac{P(\bf x|C_{i})P(C_{i})}{P(\bf x)}$

この式から分かる通り、観測データ $\bf x$ が得られた時、そのデータがクラス $C_{i}$ に分類される確率を求める事が出来ます。

では、この式のそれぞれの項目に付いて、具体的に考えてみます。

観測データとクラスに付いて

識別問題における $\bf x$ は、識別対象であるデータの特徴を並べた特徴ベクトルになります。
動物の分類を行う場合には、「体長、重量、生息地、足の数」等の観測値を正規化し特徴ベクトルとします。

$C_{i}(i=1,...,N)$ は、対象に割り当てる分類クラスです。先ほど出した動物の分類の例とすると、「ライオン、鹿、シマウマ、・・・」となります。 $C_{i}(i=1,...,N)$ の添字 $i$ は、クラスを指す変数となります。「ライオン、鹿、シマウマ、・・・・」という分類クラスがあるとしたら、 $C_{i=1}$ はライオンになります。

事前確率について

$P(C_{i})$ は、クラス $C_{i}$ の生起確率となります。生起確率とは観測データ $\bf x$ に関係なく $C_{i}$ が発生する確率です。

病気の診断を例にすると、問診の内容を観測データとした時、問診に関係なく単純に健康な人の比率と病気の人の比率を考えている事になります。
$P(C_{i})$ が非常に小さいと言う事は、病気で言うと１万人、１０万人に１人の様に、滅多に現れないという事を表しています。

また、観測データ $\bf x$ の生起確率は、観測データ全体の中から、 $\bf x$ が観測される確率なのですが、識別問題においては $\bf x$ に対して、分類クラス $C_{i}$ の全ての $i$ について、比較を行う為、 $i$ を変えても $P(\bf x)$ は変化しません。その為、識別問題においては分母の $P(\bf x)$ は省略されることが多いです。

尤度（条件付き確率)

$P(\bf x|C_{i})$ は、クラス $C_{i}$ に分類されたデータ $\bf x$ の確率分布で、離散的な事象の場合と連続な事象の場合で話が少し異なってきます。
まず、離散的な事象を病気の診断を例に考えてみると、 $C_{1}$ を病気とし、 $\bf x$ を $x_{1}=1$ を目が虚ろ、 $x_{1}=0$ を目が虚ろではないとします。
手元にあるカルテを眺めてみた時、病気と診断された人について、「目が虚であった」と観測された人の比率が $P(x_{1}=1 | C_{1})$ となります。

次に、連続的な事象を動物の分類の例に考えてみると、 $C_{1}$ をライオンとし、 $x_{1}$ を体長とします。
先ほどと同様に、ライオンと判断された動物の体長の確率分布を $P(x_{1}|C_{1})$ となります。

識別方法

観測したデータ $\bf x$ について、全てのクラス $C_{i}(i=1,..,N)$ について、事後確率 $P(C_{i}|\bf x)$ を求めます。
そしてその中から最も確率の高いクラスを識別結果とします。式にすると次のようになります。
$arg \max_{i=0,..,N} P(C_{i}|\bf x)$

アヤメデータとは

アヤメデータとは３種のアヤメ（セトーサ、ベルシカラー、バージニカ）をそれぞれ50個ずつ花弁とガクの長さと幅をサンプリングしたデータであり、分類器のベンチマークデータとしてよく登場します。
アヤメデータを花弁とガクの長さ、幅の4次元ベクトルから、2次元取り出してプロッティングすると次のようになります。
f:id:emoson:20150306141311p:plain

このアヤメデータの4次元の特徴ベクトルを観測ベクトルとし、識別器を作成します。

事前確率について

アヤメデータは3種のアヤメがそれぞれ等しい数存在しているため、生起確率 $P(C_{i})$ は等しくなります。

尤度について

アヤメデータの特徴ベクトルは連続値を取るため、確率密度関数によって尤度を表現します。

Pythonによる実装

__author__ = 'emoson'
import matplotlib.pyplot as plt
import numpy as np


def mnd(_x, _mu, _sig):
    """多次元正規分布"""
    x = np.matrix(_x)
    mu = np.matrix(_mu)
    sig = np.matrix(_sig)
    a = np.sqrt(np.linalg.det(sig)*(2*np.pi)**sig.ndim)
    b = np.linalg.det(-0.5*(x-mu)*sig.I*(x-mu).T)
    return np.exp(b)/a


def load(ayame_id, feature_list, file_path="./iris.dat.txt"):
    """アヤメデータの読み込み"""
    with open(file_path) as f:
        if type(feature_list) is list:
            return [[int(line.split(" ")[feature]) for feature in feature_list] for line in f if int(line.split(" ")[-2]) == ayame_id]
        else:
            return [int(line.split(" ")[feature_list]) for line in f if int(line.split(" ")[-2]) == ayame_id]


def mean_list(*vector):
    """平均ベクトルを返す"""
    return [sum(vec)/len(vec) for vec in vector]


def variance(vector1, vector2):
    """分散を求める"""
    mean = mean_list(vector1, vector2)
    return sum([(x-mean[0])*(y-mean[1]) for x, y in zip(vector1, vector2)])/len(vector1)


def variance_covariance_matrix(*vectors):
    """分散共分散行列を返す"""
    return [[variance(vec1, vec2) for vec2 in vectors] for vec1 in vectors]

if __name__ == "__main__":
    #識別に使う特徴 0:ガクの長さ 1:ガクの幅 2:花弁の長さ 3:花弁の幅
    test_feature_list = [0, 1, 2, 3]
    vectors = list()
    test_vectors = [load(ayame_id=ayame_id, feature_list=test_feature_list) for ayame_id in range(1, 4)]

    #特徴ベクトルの読み込み
    for ayame_id in range(3):
        vec_list = list()
        for feature_id in test_feature_list:
            vec_list.append(load(ayame_id=ayame_id+1, feature_list=feature_id))
        vectors.append(vec_list)

    #平均ベクトルの作成
    means = [mean_list(*vector) for vector in vectors]

    #識別率の確認
    true_count = 0
    false_count = 0
    for ayame_id in range(3):
        for vec in test_vectors[ayame_id]:
            #尤度を求める
            prob = [mnd(vec, _mu=mean_list(*vectors[test_id]), _sig=variance_covariance_matrix(*vectors[test_id])) for test_id in range(3)]
            #最も尤度の高いクラスを分類結果とする
            if prob.index(max(prob)) == ayame_id:
                true_count += 1
            else:
                false_count += 1
                
    print("識別率: %f" % (true_count / (true_count + false_count)))