2018年3月18日日曜日

Python3とPyMC3でベイズ推定

PythonにはPyMC3というベイズ統計モデリングと確率論的機械学習のためのパッケージがある。ベイズ推定の勉強のためにPyMC3の環境を作成し、コイントスで表が出る確率をベイズ推定してみる。


環境


Bash on Ubuntu on Windows


Jupyter Notebook



PyMC3のインストール


はじめにpython3-tkとliblapack-devをインストールしておく。そうしておかないとPyMC3のインストールに失敗する。


続いてpipでPyMC3のインストール。PyMC3といっしょにTheano、NumPy、SciPy、Matplotlibも自動でインストールされる。


バージョンの確認。


PyMC3の機能をフルに使うためにpandasとpatsyもインストールしておくと良いらしい。



コイントスの試行データ生成


コイントスをして表が出るか裏が出るかというような2つのうちどちらかしか起こらない事象をベルヌーイ試行という。ベルヌーイ分布はその試行結果を0と1の2値で表したもので、表が出る確率をpとすると、裏が出る確率は1-pとなる。

まずはコイントス500回の試行データをベルヌーイ分布に従う確率変数から生成し、グラフ化してデータを確認する。普通のコインを想定するので確率は0.5。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%matplotlib inline
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
 
# 試行回数
n = 500
# ベルヌーイ分布に従う500の確率変数をランダムに取得
# 確率は0.5
trials = stats.bernoulli.rvs(0.5, size=n)
# 表（1）の数
h = sum(trials)
 
# コインの表（1）裏（0）の度数をグラフ化
plt.subplots(figsize = ( 11 , 9 ))
plt.bar([0,1], [h, n-h])
plt.tick_params(
    axis='x',
    which='both',
    bottom='off',
    top='off',
    labelbottom='on')
plt.xticks([0, 1], ['表（1）', '裏（0）'])
plt.show()

以下のようなグラフが作成される。確率は0.5なので表、裏はほぼ同数になる。


コイントスで表が出る確率の推定


コイントスで表が出る確率をPyMC3で推定してみる。ベイズ推定についてはインターネットや書籍にたくさん情報があるのでそちらにまかせるとして（きちんと説明できるほど理解できていない）、とりあえずは、コイントスのモデルを考える必要がある。

コイントス（ベルヌーイ試行）をn回行ったとき、表が出る回数が従う確率分布が二項分布なので、二項分布を採用する。PyMC3ではBinominalクラスで二項分布を扱える。引数のnは試行回数で、今回は500で固定。もうひとつの引数である確率pを不明なものとして推定する。

確率pの値を推定するために、このpにも分布を当てはめる。今回はpについての情報は何もないものとして一様分布を採用する。PyMC3ではUniformクラスを使う。引数には最小値lowerと最大値upperが設定できるが、pは確率なので最小値0、最大値1とする。コードは以下のようになる。マルコフ連鎖モンテカルロ法（MCMC法）のアルゴリズムとしてメトロポリス・ヘイスティング法を使用した。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
%matplotlib inline
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc3 as pm
 
# 試行回数
n = 500
# ベルヌーイ分布に従う500の確率変数をランダムに取得
# 確率は0.5
trials = stats.bernoulli.rvs(0.5, size=n)
# 表（1）の数
h = sum(trials)
 
with pm.Model() as model:
    # 表が出る確率の分布（一様分布）
    # 確率なので範囲は0.0～1.0
    p = pm.Uniform('p', lower=0.0, upper=1.0)
     
    # 表が出る回数の分布（二項分布）
    # n: ベルヌーイ試行の回数
    # p: 試行の成功確率
    # observed: 観測値（ここではプログラムで生成した確率変数を使用）
    y = pm.Binomial('y', n=n, p=p, observed=h)
 
    # サンプルの生成
    # 10000: 求めるサンプル数（バーンインは含まない）
    # tune: バーンイン
    # step: MCMC法のアルゴリズム
    # njobs: チェーン数
    trace = pm.sample(10000, tune=2000, step=pm.Metropolis(), njobs=3)
 
# サンプルのプロット
pm.traceplot(trace)
plt.show()
 
# サンプルの要約統計量
print(pm.summary(trace))

結果は以下のようになった。表が出る確率pの事後分布の平均（mean）はほぼ0.5になっており、コイントス試行データ生成時の確率0.5と一致している。95％信用区間（hpd_2.5、hpd_97.5）は0.46～0.54と約0.1の範囲にあるので、有効な結果と言えそう。推定値が収束しているかどうかの指標であるRハットの値（Rhat）は、1.1より小さい値なら問題ないとされている。今回は1.000098なので問題なさそう。






投稿者 時刻:
ラベル:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)