Python
はじパタ

はじパタのベイズの識別規則をpythonで実装する

More than 1 year has passed since last update.

はじパタP23 ~ P24付近のやつです。
例題を理解するためにまとめてみました。
なので、コードの拡張性は皆無です。

ベイズの識別規則

数式はこんな感じです。

P(C_i|x) = \frac{p(x|C_i)P(C_i)}{p(x)}
事後確率 = \frac{クラス条件付き確率*事前確率}{周辺確率}
識別クラス = arg \; max \; P(C_i|x)

例題

はじパタから具体例を示していきます。
ある町の1000人をランダムにサンプルした仮想的なデータを下表に示します。

サンプル数 喫煙する人
($S=1$)
飲酒する人
($T=1$)
健康な人($G=1$) 800人 320人 640人
健康でない人($G=0$) 200人 160人 40人

健康な人($G=1$)800人と健康でない人($G=0$)200人が得られ、それぞれのクラスで喫煙と飲酒の有無を調査した結果が示されています。
ここでは、喫煙ありを$S=1$、なしを$S=0$で、飲酒ありを$T=1$、なしを$T=0$で表します。
目的は、このデータをもとにこの町の他の人の喫煙と飲酒の有無から、その人の健康状態を予測するためのベイズの識別規則を導くことです。
実際にパラメーターを当てはめるとこんな感じです。

P(G|S, T) = \frac{P(S, T|G)P(G)}{P(S, T)}

準備

導出にあたり、上表から事前確率$P(G)$とクラス条件付き確率$P(S|G), P(T|G)$を愚直に求めます。
まず、事前確率$P(G)$です。

健康な人、健康でない人

P(G=1)=\frac{800}{1000}=\frac{4}{5},\; P(G=0)=\frac{200}{1000}=\frac{1}{5}

同様にクラス条件付き確率$P(S|G), P(T|G)$を求めます。

喫煙かつ健康な人、喫煙せず健康な人

P(S=1|G=1)=\frac{320}{800}=\frac{2}{5},\; P(S=0|G=1)=\frac{480}{800}=\frac{3}{5}

喫煙かつ健康でない人、喫煙せず健康でない人

P(S=1|G=0)=\frac{160}{200}=\frac{4}{5},\; P(S=0|G=0)=\frac{40}{200}=\frac{1}{5}

飲酒かつ健康な人、飲酒せず健康な人

P(T=1|G=1)=\frac{640}{800}=\frac{4}{5},\; P(T=0|G=1)=\frac{160}{800}=\frac{1}{5}

飲酒かつ健康でない人、飲酒せず健康でない人

P(T=1|G=0)=\frac{40}{200}=\frac{1}{5},\; P(T=0|G=0)=\frac{160}{200}=\frac{4}{5}

これらの事前確率とクラス条件付き確率をそのままコードへぶち込みます。

    # 健康
    GOOD = 1
    BAD = 0

    # 喫煙
    smoking = 1
    nonSmoking = 0

    # 飲酒
    drinking = 1
    nonDrinking = 0

    # 健康な確率
    _prob_helth = []
    _prob_helth.insert(GOOD, 800/1000)
    _prob_helth.insert(BAD, 200/1000)

    # 健康と喫煙の条件付き確率
    _prob_smoke_G = []
    _prob_smoke_G.insert(smoking, 320/800)
    _prob_smoke_G.insert(nonSmoking, 480/800)
    _prob_smoke_B = []
    _prob_smoke_B.insert(smoking, 160/200)
    _prob_smoke_B.insert(nonSmoking, 40/200)

    _prob_smoke = []
    _prob_smoke.insert(GOOD, _prob_smoke_G)
    _prob_smoke.insert(BAD, _prob_smoke_B)

    # 健康と飲酒の条件付き確率
    _prob_drink_G = []
    _prob_drink_G.insert(drinking, 640/800)
    _prob_drink_G.insert(nonDrinking, 160/800)
    _prob_drink_B = []
    _prob_drink_B.insert(drinking, 40/200)
    _prob_drink_B.insert(nonDrinking, 160/200)

    _prob_drink = []
    _prob_drink.insert(GOOD, _prob_drink_G)
    _prob_drink.insert(BAD, _prob_drink_B)

分類する

準備が整ったので、喫煙($S=1$)かつ飲酒($T=1$)する人の健康状態を分類してみます。

    # どんな人??
    human = [smoking, drinking]

ベイズの識別規則に当てはめた状態が下式です。

P(G=1|S=1, T=1) = \frac{P(S=1, T=1|G=1)P(G=1)}{P(S=1, T=1)}

クラス条件付き確率$P(S=1, T=1|G=1)$と周辺確率$P(S=1, T=1)$を分解しながら求めていきます。

まず、クラス条件付き確率$P(S=1, T=1|G=1)$から見ていきます。
クラス条件付き確率$P(S=1, T=1|G=1)$は$S$と$T$の間に条件付き独立が成り立っていると仮定すると、下式のように書くことができます。

P(S=1, T=1|G=1)=P(S=1|G=1)P(T=1|G=1)

上式の右辺$P(S=1|G=1)$と$P(T=1|G=1)$は先ほど求めたクラス条件付き確率ですね。
具体的に数値を当てはめます。

P(S=1, T=1|G=1)=P(S=1|G=1)P(T=1|G=1)=\frac{2}{5} * \frac{4}{5}=\frac{8}{25}
P(S=1, T=1|G=0)=P(S=1|G=0)P(T=1|G=0)=\frac{4}{5} * \frac{1}{5}=\frac{4}{25}

コードにするとこんな感じです。
パラメーターに応じて元々格納していたクラス条件付き確率を掛け算しているだけですね。

    # 健康なクラス条件付き確率
    goodConditionalProb = exeConditionalProb(human, GOOD, _prob_smoke, _prob_drink, _prob_helth)

    # 健康でないクラス条件付き事後確率
    badConditionalProb = exeConditionalProb(human, BAD, _prob_smoke, _prob_drink, _prob_helth)

    # クラス条件付き確率を計算する
    def exeConditionalProb(human, helth, _prob_smoke, _prob_drink, _prob_helth):
        conditionalProb = _prob_smoke[helth][human[0]] * _prob_drink[helth][human[1]]
        return conditionalProb

次に周辺確率を求めます。
周辺確率は同時確率の和で求められます。

P(S=1, T=1)=P(S=1, T=1, G=1)+P(S=1, T=1, G=0)

となると、同時確率、$P(S=1, T=1, G=1)$と$P(S=1, T=1, G=0)$が欲しいわけですが、これはそもそもクラス条件付き確率の定義を見ると、

P(S, T|G)=\frac{P(S, T, G)}{P(G)}

が、成り立つので、変形させると、

P(S, T, G)=P(S, T|G)P(G)

と、なり、先ほど求めたクラス条件付き確率を使うことで求めることができます。
具体的に数値を当てはめてみます。

P(S=1, T=1, G=1)=P(S=1, T=1|G=1)*P(G=1)=\frac{8}{25}*\frac{4}{5}=\frac{32}{125}
P(S=1, T=1, G=0)=P(S=1, T=1|G=0)*P(G=0)=\frac{4}{25}*\frac{1}{5}=\frac{4}{125}
P(S=1, T=1)=\frac{32}{125}+\frac{4}{125}=\frac{36}{125}

コードはこうです。
先ほどのクラス条件付き確率を使って同時確率を求め、それらを足すことで周辺確率を求めます。

    # 健康な同時確率
    goodJointProb = exeJointProb(goodConditionalProb, GOOD, _prob_helth)

    # 健康でない同時確率
    badJointProb = exeJointProb(badConditionalProb, BAD, _prob_helth)

    # 周辺確率
    marginalProb = goodJointProb + badJointProb

    # 同時確率 = クラス条件付き確率 * 事前確率 を計算する
    def exeJointProb(conditionalProb, helth, _prob_helth):
        jointProb = conditionalProb * _prob_helth[helth]
        return jointProb

以上で必要な値が揃ったのでベイズの識別規則に突っ込みます。

P(G=1|S=1, T=1) = \frac{P(S=1, T=1|G=1)P(G=1)}{P(S=1, T=1)}=\frac{\frac{32}{125}}{\frac{36}{125}}=\frac{8}{9}
P(G=0|S=1, T=1) = \frac{P(S=1, T=1|G=0)P(G=0)}{P(S=1, T=1)}=\frac{\frac{4}{125}}{\frac{36}{125}}=\frac{1}{9}

上式の結果から、喫煙($S=1$)かつ飲酒($T=1$)する人は健康($G=1$)と分類されます。

コードは既に計算した同時確率と周辺確率を使って事後確率を求めるだけですね。
その後、大きい方の値を分類結果として返します。

    # 事後確率
    goodPosterior = goodJointProb / marginalProb
    badPosterior = badJointProb / marginalProb

    # 大きい方へ分類する
    result = goodPosterior if goodPosterior >= badPosterior else badPosterior

以上となります。
はじめにhumanで定義したパラメータを変えることでその人が健康か健康でないかを分類することができます。

最後に

結果的にパワーコードで押し切る形になってしまいましたが、参考にさせていただいたまとめでもおっしゃっているように、コードにしてみることで理解できていること、できていないことがハッキリしました。
役に立つかはわかりませんが、全体のコードはgitHubに置いておきます。

参考

実装に当たって下記まとめを参考にさせていただきました。
最も基本的なベイズ定理のpython実装