PythonにはPyMC3というベイズ統計モデリングと確率論的機械学習のためのパッケージがある。ベイズ推定の勉強のためにPyMC3の環境を作成し、コイントスで表が出る確率をベイズ推定してみる。
Bash on Ubuntu on Windows
Jupyter Notebook
はじめに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。
以下のようなグラフが作成される。確率は0.5なので表、裏はほぼ同数になる。
コイントスで表が出る確率をPyMC3で推定してみる。ベイズ推定についてはインターネットや書籍にたくさん情報があるのでそちらにまかせるとして（きちんと説明できるほど理解できていない）、とりあえずは、コイントスのモデルを考える必要がある。
コイントス（ベルヌーイ試行）をn回行ったとき、表が出る回数が従う確率分布が二項分布なので、二項分布を採用する。PyMC3ではBinominalクラスで二項分布を扱える。引数のnは試行回数で、今回は500で固定。もうひとつの引数である確率pを不明なものとして推定する。
確率pの値を推定するために、このpにも分布を当てはめる。今回はpについての情報は何もないものとして一様分布を採用する。PyMC3ではUniformクラスを使う。引数には最小値lowerと最大値upperが設定できるが、pは確率なので最小値0、最大値1とする。コードは以下のようになる。マルコフ連鎖モンテカルロ法（MCMC法）のアルゴリズムとしてメトロポリス・ヘイスティング法を使用した。
結果は以下のようになった。表が出る確率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なので問題なさそう。
環境
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。
%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法）のアルゴリズムとしてメトロポリス・ヘイスティング法を使用した。
%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 件のコメント:
コメントを投稿