【ベイズ推論】コイン投げの実例を通してPymcの使い方を掴む

Posted: , Category: ベイズ統計 , 機械学習 , 確率 , 統計学

pymcはベイジアンモデリングをPythonで行うためのPPL(確率的プログラミング言語)のフレームワークです。

ベイジアンモデリングは近年機械学習の有力な手法として利用されていますが、解析的に解けるわずかな例を除いて解析的に解くことができず、基本的にはMCMCか変分推論のアルゴリズムを用いて近似的にモデルの事後分布や予測分布を求めることになります

pymcはこの近似計算をするためのPythonライブラリです。

似ているライブラリとして、Rに付属するライブラリであるstanや、そのPythonラッパーであるpystanなどが利用されていますが、pymcはPython環境だけで、MCMCによるベイズ推論をすることができるので、大変有用かと思います。Pythonで他に利用されているものとして、Pyroなどもありますが、こちらについてはまた別の機会に解説します。

また、今回は説明しませんが、PyMCでは、MCMCの学習アルゴリズムだけでなく、変分推論(Variational Inference)のAPIも整備されており、自動推論変分推論(ADVI)などのアルゴリズムも利用可能です。

今回は、簡単なコイン投げの実例(ベルヌーイ試行)を通して、Pymcの使い方や、確率的プログラミング言語を利用したベイジアンモデリングの実装方法を理解できるようにすることが目的の記事となっています。

本記事の内容
  • コイン投げの思考をベイズ推論する
  • Pymcを用いてコイン投げの試行におけるパラメータを推定する

コイン投げの試行 〜今回行うベイズ推論〜

今回の記事の目的は、簡単な例を通してPymcの動かし方を学ぶことですので、具体例もできるだけ簡単なものにします。

今回は、コイン投げにおける試行の例を通して、ベイズ推論を行なってみます。

コイン投げの試行とパラメータ推定

今、表と裏が出る確率がわからないコインがあります。

このコインを100回振った時の、表と裏がでた回数がある時、このコインの表、裏が出る確率をベイズ手法で求めたい。

この辺りの問題設定や、共役事前分布を用いて解析的にベイズ推論のパラメータを求める方法については、次の記事をご覧ください。

[ベイズ推論] 数式と具体例でベイズ手法を用いた学習と予測を理解する
今回は、機械学習で非常によく用いられているベイズ手法による学習と予測の理論的背景を理解するために、前提知識とベイズ推論にための学習と予測の方法について一から百までまとめます。 今回を記事を執筆するあたっての参考書は、「ベ […]

コイン投げ思考をベイズ的に定式化する

まず、Pymcを使う前に、このコイン投げの試行をベイズ的に定式化しましょう。

今、コインの面が出る確率を$\mu$とすると、1回の試行でどの面が出る確率分布$p(x)$は、ベルヌーイ分布を用いて次のように表現することができます。

\begin{equation}
\begin{split}
p(x | \mu) &= Bern(x | \mu) \\
&= \mu^{x} (1 - \mu)^{1-x}
\end{split}
\end{equation}

ここで、このパラメータ$\mu$がデータから確率的に定まるようにしたいので、パラメータ$\mu$も確率変数として表現します。ベルヌーイ分布のパラメータは、ベータ分布(Beta Distribution)を用いて次のように表現できます。

\begin{equation}
\begin{split}
p(\mu | a, b ) &= Beta(\mu | a, b) \\
&=  C \mu^{a -1} (1 - \mu) ^ {b-1}
\end{split}
\end{equation}

ここで、Cは正規化定数です。ベータ分布の詳しい解説は、こちらの記事をご覧ください。

ベータ分布の期待値や分散、統計的性質について全部まとめる
ベータ分布は、値が0~1の間を取りうる確率分布で、機械学習や統計学の応用分野で非常に多く登場する確率分布です。 ベータ分布とよく混同しやすいガンマ分布は、取りうる値が正の実数の範囲内なので、ベータ分布とガンマ分布と混同し […]

今回、得られているデータを$\bm{X}$とし、その個数を$N$とします。そして、$i$番目のデータを$X_i$とすると、今回考えているベイズモデルのグラフィカルモデルはこのようになります。

また、今回求めたいパラメータ$\mu$の分布は、ベイズの定理を用いると次のようになります。

今回求めたいパラメータ$\mu$の分布
\begin{equation}
\begin{split}
p(\mu | \bm{X}) = \frac{p(\bm{X} | \mu )p(\mu)}{p(\bm{X}) } \\
= \frac{p(\bm{X} | \mu )p(\mu)}{ \int p(\bm{X}, \mu) d\mu }
\end{split}
\end{equation}

(3)の分布は、実はベルヌーイ分布とその共役事前分布であるベータ分布を利用すると解析的に解けることが知られていますが、今回は、あえてpymcを使って近似的にこれらを解いてみます。

Pymcを用いて、パラメータを推論

それでは、Pymcを用いてベイズモデルにおけるパラメータを推論していきましょう。

Pymcでモデルを組む場合は、確率変数における依存関係をグラフィカルモデルで表現すると、実装が簡単になるのでおすすめです。今回考えるベイズモデルのグラフィカルモデルはこのようになっています。

$X_i$がN回のコイン投げのうちの$i$番目の結果を示しており、パラメータ$\mu$のベルヌーイ分布によってこれらの結果が得られたとする生成モデルを考えています

ベイジアンモデリングにおいては、データの生成過程をこのようにグラフィカルモデルに落とし込むことができたらほぼ勝ち確定です。あとは、PyMCなどの言語でモデリングして推論すれば大丈夫です。

それでは実装に入っていきます。

PyMCでベルヌーイ試行のベイジアンモデリングを実装

まず最初に、必要なライブラリをまとめてインストールします。

import numpy as np
from scipy import stats
import pymc3 as pm
import arviz

arvizが見慣れないパッケージだと思いますが、これはMCMCの結果を可視化するためのライブラリで、MCMCなどをやっていくと頻繁に使うライブラリです。

続いて、事後分布の学習に用いるテストデータを生成します。

今回は、コインの表面が出る確率が0.7であるコインを用いた設定でテストデータを生成します。

実際にPymcで学習した結果、0.7に近い値が出力できれば、うまく学習できていると確かめることができます。

今回は1000個のデータを用意します。データの生成には、statsライブラリを利用します。

N = 1000
mu_real = 0.7
data = stats.bernoulli.rvs(mu_real, size=N)

続いてモデルの定義と学習を行います。ここが結構PyMC独自の記法が登場します。

まず、コードを先に掲載します。

with pm.Model() as model:
    # 事前分布の設定, ベルヌーイ分布のパラメータをp1(0~1)を定義
    p1 = pm.Uniform('p1', lower=0.0, upper=1.0)
    
    # ベルヌーイ分布を定義。observerdにデータを定義
    y = pm.Bernoulli("obs_A", p1, observed=data)

上のコードでは、pm.Model オブジェクトを作成しています。このオブジェクトが、pymcで実際に定義する統計モデルの本体です。

with pm.Model でブロックを作っており、そのブロック内で、統計モデルの中に登場する確率分布やそのパラメータを記述します。

今回はベルヌーイ分布の定義と、そのパラメータを、0~1の範囲における一様事前分布をpm.Uniformで与えています。今回、ベルヌーイ分布のパラメータ$\mu$を求める必要がありますが、その事前分布を何かしらの形で置いてあげる必要があるため、今回は一様分布で与えています。

ここまでで、モデルの定義ができました。

きちんと意図通りに、ベイジアンモデリングができているか確認するために、pm.Modelの中身を確認しましょう。

model.basic_RVs
#=> [p1_interval__ ~ TransformedDistribution, obs_A ~ Bernoulli]

このようになりました。obs_A は観測値で、観測値がベルヌーイ分布に従い、そのパラメータp1_interval__ が、TransformedDistributionという形でモデリングされていることがわかりました。

このTransformedDistributionの意味はひとまず放置して大丈夫です。

ここまででベイジアンモデルリングのモデル定義は終了です。続いて、MCMCのアルゴリズムで、パラメータ$\theta$の事後分布からサンプリングを行なっていきましょう。コードは次にようになります。

with model:
    trace = pm.sample(
        draws=10000,
        tune=1000,
        step=pm.NUTS(),
        chains=3,
        random_seed=128,
        return_inferencedata=True
    )

モデリングの時と同じように、withを用いてpm.Modelオブジェクトのコンテキスト内でサンプリングを行います。(この辺りは、Pymc特有の記法なのですが、何度か使っているうちに勝手に覚えていきます)

sample関数の引数ですが、drawsがサンプリング回数です。tuneは、最初の何個のデータを棄却するかを指定します。MCMCのサンプリングの場合、最初は精度が出ないので、精度がでないデータを何個か捨てるようなことを通常よく行います。

step引数は、MCMCのアルゴリズムで、今回はNUTSを指定しています。return_inferencedata引数は、データの返り値をarviz形式にするかどうかを指定できます。今回はこちらをTrueにしています。

これで、パラメータ$\theta$からのサンプリングは完了です。コードを実行すると時間がかかりますが、しばらくすると実行が終わります。

実行が終わったあとは、次のようにarvizを用いて、事後分布の挙動を確認します。

arviz.plot_trace(trace)

このような結果が出ました。

今回は10000点のデータを事後分布からサンプリングしましたが、p1の分布の形状として、0.7あたりが最大値となっており、うまくパラメータを推論できていることがわかります。

また次のコードで、平均や信頼区間などの値も出力することができます。

arviz.summary(trace)

また事後分布だけを描写する場合は、次のarvizコードを利用します。

arviz.plot_posterior(trace)

【広告】
統計学的にあなたの悩みを解決します。
仕事やプライベートでお悩みの方は、ベテラン占い師 蓮若菜にご相談ください。

機械学習と情報技術