ベイズの識別規則をRで実装してみた

今学期は授業数が少ない代わり，latex.ltxリーディング，構造生物学論文輪読，そして機械学習の3つの自主ゼミを掛け持ちしているのですが，数学がまったくできないのに機械学習なんぞに手を出したのでそこそこ苦労しています．

はじめてのパターン認識

機械学習のゼミは『はじめてのパターン認識』なる本をテキストにしているのですが，第3章あたりで既に数式アレルギーを発症しました．そんなとき，ある友人の言葉を思い出しました．

「よくわからなかったら，実装しろ」

ということで，『はじめてのパターン認識』第3章で登場するベイズの識別規則を実装してみました．因みに，これが私にとっては"はじめてのR"（単なる電卓としての利用以外）でした．

ベイズの識別規則とは

これをきちんと説明しようとするとこの記事が異様に長くなるところだったのですが，幸いにして既にまとめてくれてる方がいらっしゃいました．

はじめてのパターン認識 第3章 ベイズの識別規則 解説

詳しい説明は上記リンク先にあるので，ここではごくごく簡単な説明に留めます．

まず，観測データが与えられた下で，それがクラスに属する条件付き確率を事後確率といい，次式で求められます．

この式はベイズの定理と呼ばれています．ベイズの識別規則は，この事後確率が最も大きくなるクラスに観測データを分類する識別規則です．

また，識別を誤ったときに生じる損失（危険性）はクラス間で対称であるとは限らないので，ベイズの識別規則に損失（本当はに属するものをと識別することによる損失）の概念を導入することがあります．このとき，損失の期待値が最小となるようにした識別規則が最小損失基準に基づくベイズの識別規則です．

Rでの実装

さきほどのリンク先にも掲載されていますが，『はじめてのパターン認識』で登場したサンプルデータをそのまま借用し，次のデータに基いて識別規則を作成することにします．

	サンプル数	喫煙する（）	飲酒する（）
健康な人（）	800人	320人	640人
健康でない（）	200人	160人	40人

それでは，ベイズの識別規則と最小損失基準に基づくベイズの識別規則を返す関数を定義します．説明されている通りに計算させておけば問題ありません．

# ベイズの識別規則
byzeIR <- function(U,C1,C2) {
  # 事前確率
  P_C1 <- C1[1]/U
  P_C2 <- C2[1]/U
  
  # 各パラメタに対するクラス付き条件確率：(X,C)=(1,1),(0,1),(1,2),(0,2)
  P_S <- c(C1[2]/C1[1],1-(C1[2]/C1[1]),C2[2]/C2[1],1-(C2[2]/C2[1]))
  P_T <- c(C1[3]/C1[1],1-(C1[3]/C1[1]),C2[3]/C2[1],1-(C2[3]/C2[1]))
  
  # クラス付き条件確率：(S,T)=(1,1),(0,1),(1,0),(0,0)
  P_ST.C1 <- c(P_S[1]*P_T[1],P_S[2]*P_T[1],P_S[1]*P_T[2],P_S[2]*P_T[2])
  P_ST.C2 <- c(P_S[3]*P_T[3],P_S[4]*P_T[3],P_S[3]*P_T[4],P_S[4]*P_T[4])
  
  # 同時確率：(S,T)=(1,1),(0,1),(1,0),(0,0)
  P_STC1 <- P_ST.C1*P_C1
  P_STC2 <- P_ST.C2*P_C2
  
  # 周辺確率：(S,T)=(1,1),(0,1),(1,0),(0,0)
  P_ST <- P_STC1+P_STC2
  
  # 事後確率：(S,T)=(1,1),(0,1),(1,0),(0,0)
  P_C1.ST <- P_ST.C1*P_C1/P_ST
  P_C2.ST <- P_ST.C2*P_C2/P_ST
  
  # 識別規則
  IR <- c(
    if (P_C1.ST[1] > P_C2.ST[1]) {"C1"} else {"C2"},
    if (P_C1.ST[2] > P_C2.ST[2]) {"C1"} else {"C2"},
    if (P_C1.ST[3] > P_C2.ST[3]) {"C1"} else {"C2"},
    if (P_C1.ST[4] > P_C2.ST[4]) {"C1"} else {"C2"}
    )
  
  return(IR)
}

# 最小損失基準に基づくベイズの識別規則
byzeSLIR <- function(U,C1,C2,L) {
  SLIR <- c(
    if (P_ST.C1[1]/P_ST.C2[1] > (L[2]-L[3])*P_C2/((L[4]-L[1])*P_C1)) {"C1"} else {"C2"},
    if (P_ST.C1[2]/P_ST.C2[2] > (L[2]-L[3])*P_C2/((L[4]-L[1])*P_C1)) {"C1"} else {"C2"},
    if (P_ST.C1[3]/P_ST.C2[3] > (L[2]-L[3])*P_C2/((L[4]-L[1])*P_C1)) {"C1"} else {"C2"},
    if (P_ST.C1[4]/P_ST.C2[4] > (L[2]-L[3])*P_C2/((L[4]-L[1])*P_C1)) {"C1"} else {"C2"}
  )
  
  return(SLIR)
}

# 全体集合
U <-1000

# クラスデータ：Cx=c(ALL,S,T)
C1 <- c(800,320,640)
C2 <- c(200,160,40)

# 損失基準：L=c((C1,C1),(C1,C2),(C2,C2),(C2,C1))
L <- c(0,20,0,1)

# 識別規則を表示：(S,T)=(1,1),(0,1),(1,0),(0,0)
byzeIR(U,C1,C2)     #=> "C1" "C1" "C2" "C1"
byzeSLIR(U,C1,C2,L) #=> "C2" "C1" "C2" "C2"

感想

2クラスにしか対応していないからかもしれないですが，正直これが何の役に立つのかよくわかりません．何かもう1つ実例を出そうと思ったのですが，何も思いつきませんでした（何か面白い例を思いついたら追記するかもしれません）．

ところで，最後にも求まった最小損失基準に基づくベイズの識別規則はなかなかですね．

被験者「喫煙も飲酒もしていません」
医者「あなたは不健康ですね」
被験者「！？」