PyMCで学ぶベイズ統計モデリング｜イカサマコインの判別にチャレンジ

PyMCは、Pythonでベイズ統計モデリングを行うためのライブラリです。この記事では、イカサマコインの判別を、PyMCを利用してやってみます。これを通じて、回帰問題などではなく、よりベイズモデリングらしい例を通じてPyMCの使い方を解説します。

あわせて読みたい

PyMC（ベイズ統計モデリング）で回帰問題を解く｜Python

PyMCとは

PyMCは、Pythonでベイズ統計モデリングを行うためのライブラリです。PyMCは、マルコフ連鎖モンテカルロ（MCMC）サンプリングや変分推論などのアルゴリズムを使用して、ベイズ統計モデリングのための様々な手法を提供しています。

ベイズ統計モデリングでは、未知のパラメータの確率分布をモデルを使って表現します。そして、観測したデータを使って、モデルの未知のパラーメータの事後分布を推定します。これにより、分布の予測を行うのがベイズ統計モデリングになります。

PyMCでは、モデルを柔軟かつ効率的に記述する機能と、モデルを使ってパラメータを推定する機能が用意されています。このライブラリを使うと、ベイズ統計モデリングを使った推定を比較的手軽に行うことが可能です。

今回は、このライブラリを使って「イカサマコインかどうか」を判定してみます。

conda create -c conda-forge -n pymc_env "pymc>=5" "numpy=1.25.2"

PyMCによるイカサマチェック

※記事の末尾に、コード全体を置いています。コピー＆ペーストして実行してみる場合は、そちらを利用してください。

ライブラリのインポート

必要なライブラリをインポートします。今回は、「推論時に与えるコイントスの結果」をランダムに生成するために、”scipy”を利用するので、scipyもインストールしておきます。

import numpy as np
from scipy import stats
import pymc as pm
import matplotlib.pyplot as plt

コイントスの結果（仮）を生成する

今回は、「すでに100回コイントスをしていて、その結果が分かっている」とします。イカサマコインを想定しますので、ここでは、表が出やすいコインをシミュレートします。

ここでは、結果の作成にベルヌーイ分布の乱数を用いました。引数のmuが0.7なので、表と裏が7:3で出現するイカサマコインをトスした結果をシミュレーションしたことになります。

N = 100
mu = 0.7
data = stats.bernoulli.rvs(mu, size=N)

乱数を使っているの変化しますが、data = 1とdata = 0の割合は、例えば以下のようなコードで確認できます。

print(len(data), sum(data == 0), sum(data == 1))
# 100 32 68

上記の例では、おおよそ裏（０）が3, 表（１）が7になっています。

モデリング

今回のモデルは以下のようになります。コインの表・裏の出る確率pは一様分布（0~1）として、コイントスの結果はベルヌーイ分布に基づくと定義しています。

ベルヌーイ分布とは、「成功か失敗か」「表か裏か」のように２種類のみの結果しか得られないような試行の結果を0と1で表した分布です。ちょうど今回のコイントスのような状態が当てはまります。

確率pは、とりあえずどの確率も一様に発生するとしてモデル化しています。「均等に出るように作ったけど、製造時のばらつきがある」といった状態を想定する場合は、また違う分布でモデリングするのも良いかもしれません。

with pm.Model() as model:
    p = pm.Uniform('p', lower=0.0, upper=1.0)
    y = pm.Bernoulli("y", p, observed=data)

今回はモデルが単純なので、たった２行ですが、複雑なモデルの場合はモデルが大きくなります。

サンプリング

モデルが定義できたので、作成したモデルでサンプリングを行います。ここではサンプリング回数（draws）を3000回、最初の1000回を捨てる（tune=1000）ようにしています。

最初の試行は暴れやすいのでtune=xxxで捨てます

また、乱数のシードを42に初期化しています。

with model:
    trace = pm.sample(
        draws = 3000, 
        tune = 1000,
        random_seed = 42,
    )

上記コードにより、サンプリングが実行されます。

結果表示

サンプリングが完了したら結果を表示させます。

pm.plot_trace(trace)
plt.show()

pm.plot_posterior(trace)
plt.show()

print(pm.summary(trace))

plot_traceで経過を確認する

plot_trace()で、経過を確認できます。右側のグラフが試行回数毎のpの値です。これが一定のばらつきで動いている時は、モデルによる推定がうまく動いている時です。

今回のplot_traceの結果を見ると、グラフは均一に乱れているので、pはうまく推定できてそうです。

plot_posteriorで事後分布を確認する

plot_posteriorでグラフィカルに事後分布を表示することが可能です。

このグラフを見ると、おおよそ0.7付近にピークがあることがわかります。

summaryで結果を確認する

summaryを表示すると「結果のまとめ」を表示できます。

    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat
p  0.677  0.046   0.587    0.761      0.001      0.0    5773.0    7940.0    1.0

今回は「イカサマコインかどうか」と確認したいので、「信頼区間が0.5を含むかどうか」と確認します。もし、0.5を含んでいればイカサマコインの可能性は低く、0.5を含まなければイカサマコインの可能性が高いという結論になります。

96%の信頼区間を調べるには、結果のhdi_3%とhdi_97%の値を確認します。これを見ると、0.587~0.761です。この範囲が信頼区間となりますが、今回の結果では0.5は含まれていません。

従って、「このコインはイカサマコインの可能性が高い」という結果になりました。

事前の試行をイカサマコインで行ったので、この結果は正解になります。また、今回のイカサマコインは７割が表になるコインでしたが、0.7が信頼区間に含まれているので、予測した確率の範囲も正確であることがわかります。

2回目の実行結果

サンプリングをやり直した例です。この例でも、コインはイカサマコインの可能性が高い結果となりました。

イカサマがない場合は？

イカサマコインではない場合も試してみます。mu=0.5として、公平なコインをシミュレーションしてみます。

N = 100
mu = 0.5
data = stats.bernoulli.rvs(mu, size=N)

100回しかコイントスをシミュレーションしていないので、事前のコイントスの結果に偏りが発生していました。

しかしながら、結果を見ると二回とも0.50を含む結果となっており、「コインはイカサマコインとは言えない」という正しい結果を導き出せませした。

全コード

以下、全コードです。

import numpy as np
from scipy import stats
import pymc as pm
import matplotlib.pyplot as plt

N = 100
mu = 0.7
data = stats.bernoulli.rvs(mu, size=N)

with pm.Model() as model:
    p = pm.Uniform('p', lower=0.0, upper=1.0)
    y = pm.Bernoulli("y", p, observed=data)

with model:
    trace = pm.sample(
        draws = 3000, 
        tune = 1000,
        random_seed = 42,
    )


pm.plot_trace(trace)
plt.show()

pm.plot_posterior(trace)
plt.show()

print(pm.summary(trace))

まとめ

PyMCを使って、コインがイカサマコインかどうかの判定をやってみました。事前情報が少ないので偶然かもしれませんが、ある程度の予測を行うことがきるのがベイズの面白いところです。

ある/Aru

IT＆機械学習エンジニア/ファイナンシャルプランナー（CFP®)

XContact

専門分野は並列処理・画像処理・機械学習・ディープラーニング。プログラミング言語はC, C++, Go, Pythonを中心として色々利用。現在は、Kaggle, 競プロなどをしながら悠々自適に活動中