機械学習・Pythonお勉強（ベイズの識別規則を例にして）

はじめに
　これまでRを使って社会科学系の分析のみを行ってきたが、就職するにあたってPythonと機械学習のお勉強をはじめなければならないことになった。そこで、機械学習・Python関係の勉強に関してアウトプットとしてここに残しておく。現在の勉強法はData Campに登録しData Scientist with Pythonのコースを進めているのと、はじめてのパターン認識を読んでいる。ここでは、その両者を組み合わせた勉強のアウトプットを残す。具体的には、Data Campで学んだPythonの知識を活かしつつ、はじめてのパターン認識で紹介されている様々な分析方法を有名なtitanicデータで動かしてみる。Data CampのData Scientist with Pythonのコースについてはまだ序盤だし、はじめてのパターン認識もしっかりと理解しているのか怪しい箇所もあるが、とりあえずアウトプットしないと「分かった感」がないのでやってよう。今回ははじめてのパターン認識の第3章「ベイズの識別規則」で紹介されているベイズの識別規則を用いてtitanicのsurvived or notを予測してみる。

www.datacamp.com

www.morikita.co.jp

ベイズの識別規則

概要
　ここでは、簡単にベイズの識別規則の概要を示す。なお、詳細は上記の本を参照されたい。
　まず、ベイズの定理は以下の通りである。

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

　これは以下の式から簡単に導出できる。

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

　ここで、仮にAをあるサンプルが所属するクラス、Bをそのサンプルの観測データとすると、ベイズの定理の左辺 $P(A|B)$ はあるサンプルの観測データが得られた上でそのサンプルがクラスAに所属する確率を意味すると解釈できる。この確率の事を事後確率という。さらに右辺 $\frac{P(B|A)P(A)}{P(B)}$ の中の $P(B|A)$ はあるクラスに所属するサンプルから観測データが得られる確率（実際には確率とは言えないが、ここではとりあえず無視する）と考えられる。すなわち、サンプルの所属するクラスが不明な状態において、その不明な所属クラスを求めるべきパラメータとした時に観測データが得られる確率なので、これは尤度と考えることができる。また、P(A)はクラスの生起確率で事前確率と呼ぶ。例えば、サンプルにおいて何も情報がない場合、すなわち観測データが何もない場合、あるサンプルのあるクラスへの所属確率は $P(A)$ と考えるのが妥当であろう。例えば、日本人で男女の確率がそれぞれ0.5であるならば、ある日本人を1人連れてきた時、その日本人の特徴が何もわからなければ、その日本人が男である確率は0.5であるというのが妥当であろう。これは日本人の男女の確率が0.5であるという事前確率をそのまま利用している。ただし、仮にその日本人の名前という情報（観測データ）がわかれば、その名前を勘案してその確率を変化させることができる。例えば、その日本人の名前が「たろう」ならば男性である可能性が高くなる。

　このように、上記のベイズの定理は事後確率 = 尤度 × 事前確率　の形になっている。なお後に述べるように、右辺の分母 $P(B)$ はベイズの識別規則においては特に考えなくても良いのでここでは無視しているが、これは周辺確率と呼ばれる。
　仮に、あるサンプルが所属するクラスがAとaの2パターンあった場合を考えよう。あるサンプルに関する観測データBが得られた状態でそのサンプルがクラスAに所属する確率は以下のように表現できる。

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

　一方で、そのサンプルがクラスaに所属する確率は以下のように表現できる。

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

　これらの確率が大きい方が、そのサンプルが所属するクラスだと言える。見方を変えれば、これらの確率が等しいところが境界線となる（識別境界）。すなわち、識別境界は以下のように表現できる。

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

　この式から分かるように、 $P(B)$ はどちらにも共通して現れている。よって、この値が何であってもこの識別境界には影響しない。したがって、一般的にあるサンプルが所属するクラスをCとすると（ここではC = A or a）、以下の式を最大化するクラスがそのサンプルが観測データBのもとで所属するクラスと考える。これがベイズの識別規則である。

$P(B|C)P(C)$

　なお、これまでは観測データをBのみとしてきたが、これは複数でも良い。例えばBに加えてbがあればベイズの識別規則（以下の式を最大化する）は以下の通りとなる。

$P(B, b|C)P(C)$

　ここで、クラスCのもとでBとbが導出される確率が独立であるとすれば（条件付き独立）、 $P(B, b | C) = P(B|C)P(b|C)$ なので、上記の識別規則は以下の通りとなる。

$P(B|C)P(b|C)P(C)$

$\epsilon(B) = min(P(A|B), P(a|B))$

　ここで、クラスAに識別される領域を $R_1$ 、クラスaに識別される領域を $R_2$ とすると条件付きベイズ誤り率の期待値はその定義から以下のように表現できる。なお、この期待値はベイズ誤り率ともいう。

$E(\epsilon(B)) = \int_{R_1 + R_2} \epsilon(B)p(B) dB$

$\epsilon(B) = min(P(A|B), P(a|B))$ なので、上式は以下のように表現できる。

$E(\epsilon(B)) = \int_{R_1 + R_2} min(P(A|B), P(a|B)) p(B) dB$

　さらにベイズの定理より、

$E(\epsilon(B)) = \int_{R_1 + R_2} min(\frac{P(B|A)P(A)}{P(B)}, \frac{P(B|a)P(a)}{P(B)})p(B) dB$

　これを展開すると以下のようになる。

$E(\epsilon(B)) = \int_{R_1 + R_2} min( P(B|A)P(A), P(B|a)P(a) ) dB$

　ここで、 $R_1$ では必ず $P(B|A)P(A) ＞ P(B|a)P(a)$ 、 $R_2$ では必ず $P(B|a)P(a) ＞ P(B|A)P(A)$ なので、上式は以下のようになる。

$E(\epsilon(B)) = \int_{R_2}{P(B|A) P(A)} dB + \int_{R_1}{P(B|a) P(a)} dB$

　以上が、ベイズの識別規則の概要である。

損失
　次に、損失について考える。例えば医療の現場では、健康の人が病気と判断されるよりも病気の人が健康と判断されることの方が問題である。つまり、誤りを犯すことによって発生する危険性はクラス間で対象ではないのである。そこで、このような事を考慮に入れた識別規則を考える必要がある。そのために以下の損失を考える。

$L_{aA}$

　これは真のクラスがAにも関わらず、aと判断することによって被る損失である。クラス数がAとaの2つだとすると、それは $L_{aA}$ 、 $L_{aa}$ 、 $L_{Aa}$ 、 $L_{AA}$ が考えられる。
　ここで、データBが与えられた上で、本当はAにも関わらずaとした時に被る損失は以下のように表現できる。

$r(a|B) = L_{aA}P(A|B) + L_{aa}P(a|B)$

　なお、ここではクラスがAとaしかない事を前提としているが、仮にクラスがK個あった場合、上の損失は以下の通りとなる。

$r(a|B) = \sum_{k=1}^{K} L_{ak}P(C_k | B)$

　損失を考慮した時の識別規則は、この上式を最小化することである。
　再び、Aもしくはaのクラスしかない事を前提とするとその損失は以下のように表現できる。

$r(B) = min(r(A | B), r(a | B))$

　ここでその損失の期待値は期待値の定義より以下の通りである。

$E(r(B)) = \int_{R_1 + R_2} min(r(A | B), r(a | B ))P(B) dB$

　損失の定義から、

$E(r(B)) = \int_{R_1 + R_2} min( L_{aa}P(a|B) + L_{aA}P(A|B), L_{Aa}P(a|B) + L_{AA}P(A|B)) p(B) dB$

　ベイズの定理から展開して、

$E(r(B)) = \int_{R_1 + R_2} min(L_{aa}P(B | a)P(a) + L_{aA}P(B | A)P(A), L_{Aa}P(B | a)P(a) + L_{AA}P(B|A)P(A))dB$

　ここで、 $R_1$ では必ずAと判断され、 $R_2$ では必ずaと判断されるので上式は以下のようになる。

$E(r(B)) = \int_{R_1}( L_{Aa}P(B | a)P(a) + L_{AA}P(B|A)P(A)) dB + \int_{R_2}(L_{aa}P(B | a)P(a) + L_{aA}P(B | A)P(A)) dB$

　識別規則は $L_{Aa}P(B | a)P(a) + L_{AA}P(B|A)P(A)$ と $L_{aa}P(B | a)P(a) + L_{aA}P(B | A)P(A)$ を比較して、前者が小さければ、すなわちAと判断される時の損失が小さければクラスAに、後者が小さければ、すなわちaと判断される時の損失が小さければクラスaと判断する。

リジェクト

　すでに述べたように、以下が条件付きベイズ誤り率である。

$\epsilon(B) = min(P(A|B), P(a|B))$

　ここで、仮に $P(A|B)$ と $P(a|B)$ が等しい場合、その誤り率は0.5となる。すなわち半分の確率で誤った分類をしてしまう。そこで、このように誤り率が高い場合の対策として、ある閾値を決めておきそれ以上になった場合の判断を避けることがある。これをリジェクトと呼ぶ。仮に、その閾値をtとすると以下の領域となった場合にリジェクトする。

$\epsilon(B) ≧ t$

　仮に、あるサンプルがAと判断された場合は $P(B|A)P(A) ＞ P(B|a)P(a)$ が成立している。したがって、クラスAに判断される時のリジェクトの閾値をtとすると、以下のようになる。

$\epsilon(B) = \frac{P(B|a)P(a)}{P(B)} ≧ t$

分析
　以下では、上記のベイズの識別規則を使って実際にPythonでtitanicのsurvived or notを予測してみる。ここではsurvived or not（生存したら1、生存できなかったら0）が予測すべき所属クラスである。話を簡単にするため、ここではtitanicデータからランダムに重複無しでもともとのNの80%を学習用データ、20%をテストデータとしてその予測精度を評価しよう。なお、はじめてのパターン認識の第2章でも解説されているが、一般にこの方法での誤り率は真の誤り率よりも大きくなることには留意されたい。また、ここでも簡単のためにデータは性別のみを使用する。

まずはデータの読み込み・チェックや変数作成などの準備。


##### パッケージの読み込み ####
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
import random

##### データの読み込み ####
data = pd.read_csv('train.csv') # train.csvという名前でデータ保存

#### データチェック ####
data.head()
data.info()
data.describe()

## Survived ##
#度数
data['Survived'].value_counts()

## sex ##
# 度数
data['Sex'].value_counts()

# Survivedとの関係
sns.pointplot(x = data['Sex'], y = data['Survived'], ci = 95)

f:id:abcxyzonetwothree:20190815135126p:plain — Sex and Survived

　この図をみると明らかに、男性よりも女性の方が生存していることがわかる。ここから直感的に、ある人の生存の有無を尋ねられたらその人の性別を確認して男性ならば死亡、女性ならば生存と答えれば正解の可能性が高くなることがわかる。まあここからほとんど明らかだが、一応事後確率を計算してみる。


#### ベイズの識別規則 ####
N = len(data['Sex']) #サンプル数
train_data_N = round(N * 0.8) #学習データのN
test_data_N = N - learning_data_N #テストデータのN

train_data = data.sample(n = train_data_N, replace = False) #ランダムに80%サンプル
test_data = data.drop(train_data.index) #残りをtest_data

pd.crosstab(train_data['Sex'], train_data['Survived']) #table形式で確認

# 生存者と死亡者の数
n_survived = 0
n_nonsurvived = 0
for x in train_data['Survived'] :
    if x == 1 :
        n_survived += 1
    elif x == 0 :
        n_nonsurvived += 1
print('N of Survived:', str(np.sum(n_survived)))
print('N of Non-Survived:', str(np.sum(n_nonsurvived)))

# 男性と女性の数
n_female = 0
n_male = 0
for x in train_data['Sex'] :
    if x == 'female' :
        n_female += 1
    elif x == 'male' :
        n_male += 1
print('N of female:', str(np.sum(n_female)))
print('N of male:', str(np.sum(n_male)))

# 男性で生存者、死亡者、女性で生存者、死亡者の数
n_survived_male = 0
n_nonsurvived_male = 0
n_survived_female = 0
n_nonsurvived_female = 0
for x, y in zip(train_data['Survived'], train_data['Sex']) :
    if x == 1 and y == 'male':
        n_survived_male += 1
    elif x== 1 and y == 'female':
        n_survived_female += 1
    elif x == 0 and y == 'male':
        n_nonsurvived_male += 1
    elif x == 0 and  y == 'female':
        n_nonsurvived_female += 1
print('男性かつSurvived:', str(n_survived_male))
print('男性かつNon-Survived:', str(n_nonsurvived_male))
print('女性かつSurvived:', str(n_survived_female))
print('女性かつNon-Survived:', str(n_nonsurvived_female))

# 事前確率の計算
Pr_s = n_survived / N #事前確率（Surivived）
Pr_n = n_nonsurvived / N #事前確率（Non-Surivived）
print('事前確率（Surivived）:', str(Pr_s))
print('事前確率（Non-Surivived）:', str(Pr_n))

#周辺確率の計算
Pr_m = n_male / N #周辺確率（男性)
Pr_f = n_female / N #周辺確率（女性)

# 尤度の計算
Pr_m_s = n_survived_male / n_survived #尤度（Survivedで男性）
Pr_m_n = n_nonsurvived_male / n_nonsurvived #尤度（Non-Surivivedで男性）
Pr_f_s = n_survived_female / n_survived #尤度（Surivivedで女性）
Pr_f_n = n_nonsurvived_female / n_nonsurvived #尤度（Non-Surivivedで女性）
print('Survivedで男性:', str(Pr_m_s))
print('Non-Surivivedで男性:', str(Pr_m_n))
print('Surivivedで女性:', str(Pr_f_s))
print('Non-Surivivedで女性:', str(Pr_f_n))

# 事後確率の計算
Pr_s_m = (Pr_m_s * Pr_s) / Pr_m #事後確率（男性でSurvived）
Pr_s_f = (Pr_f_s * Pr_s) / Pr_f #事後確率（女性でSurvived）
Pr_n_m = (Pr_m_n * Pr_n) / Pr_m #事後確率（男性でNon-Survived）
Pr_n_f = (Pr_f_n * Pr_n) / Pr_f #事後確率（女性でNon-Survived）
print('男性でSurvived:', str(Pr_s_m))
print('女性でSurvived:', str(Pr_s_f))
print('男性でNon-Survived:', str(Pr_n_m))
print('女性でNon-Survived:', str(Pr_n_f))

これを回すと以下の通りになる。

N of Survived: 272
N of Non-Survived: 441
N of female: 254
N of male: 459
男性かつSurvived: 86
男性かつNon-Survived: 373
女性かつSurvived: 186
女性かつNon-Survived: 68
事前確率（Surivived）: 0.3052749719416386
事前確率（Non-Surivived）: 0.494949494949495
Survivedで男性: 0.3161764705882353
Non-Surivivedで男性: 0.8458049886621315
Surivivedで女性: 0.6838235294117647
Non-Surivivedで女性: 0.15419501133786848
男性でSurvived: 0.18736383442265794
女性でSurvived: 0.7322834645669292
男性でNon-Survived: 0.8126361655773421
女性でNon-Survived: 0.2677165354330709

ここから分かるように、ある人が男性だった場合その人が生存している可能性は約18.7%であり、死亡している可能性は約81.3%である。また、女性だった場合その人が生存している可能性は約73.2%であり、死亡している可能性は約26.8%である。よって、男性ならば死亡、女性なら生存と判断することがここでのベイズの識別規則からは妥当であると言えよう。
　テストデータでその正解確率を計算してみたら約79.7%であった。なお、そのコードは以下の通りである。


# test_dataで結果の評価
res = []
N_sucess = 0
N_fail = 0
for x, y in zip(test_data['Sex'], test_data['Survived']):
    if x == 'male':
        res.append(0) #男性なら死亡判定
    elif x == 'female':
        res.append(1) #女性なら生存判定
        
    if res[-1] == y:
        N_sucess += 1 #正解
    elif res[-1] != y:
        N_fail += 1 #不正解

rate = N_sucess / (N_sucess + N_fail)    
print('正解率は', str(rate), 'です')

メモ

個人的なメモです。他者にわかりやすく書くよりも未来の自分にわかりやすく書いています。なお、記事内容の正確さは保証できません。勉強中の身ですので、間違い等ご指摘頂けたら幸いです。

機械学習・Pythonお勉強（ベイズの識別規則を例にして）