AIエンジニア/データサイエンティストが使うベイズ統計モデル
Python3 PyMC3 によるMCMC（Markov chain Monte Carlo）
モデリング対象テーマ「メッセージ数に変化はあるか？」
このサイトの記事は，AIエンジニア/データサイエンティストが使う本格的なベイズ統計モデルのベストプラクティスを紹介します。
事後分布の算出にMCMC（Markov chain Monte Carlo）を使います。ベイズ統計モデルにMCMCを使うことは不可分ではないですが，ほぼデフォルトで使います。
モデリング対象テーマは，下で紹介する本にある「メッセージ数に変化はあるか？」です。
日々受信するテキストメッセージ（eメールのこと？）の数は，0件から数十件とバラバラですが，ある日を境に増えているようですが，本当にそうなのかベイズ統計モデルを使って解き明かします。
下の図が毎日の受信数とMCMCで解いた受信数の期待値の図です。明らかに期待値の特異点があります。毎日の受信数のグラフから決して読み取れないと思われます。また，古典的な統計学（頻度主義統計学ともいう非ベイズ統計学）では解けない課題かと思われます。
from IPython.display import Image
Image("./image_common/book_Davidson-Pilon_fig1.jpg")
このサイトの記事は，次の本を参考にしています。
私は，Pythonや統計学の本を50冊ほど持っていますが，その多くがプロになるための実践本ハンズオンとしては物足りたないと感じています。次の本はマチガイなくプロになるための本です。
from IPython.display import Image
Image("./image_common/book_Davidson-Pilon.jpg")
『Pythonで体験するベイズ推論』 ーPyMCによるMCMC入門ー キャメロン デビッドソン=ピロン (著), 玉木 徹 (翻訳)
from IPython.display import Image
Image("./image_common/book_Martin.jpg")
『Pythonによるベイズ統計モデリング』PyMCでのデータ分析実践ガイド Osvaldo Martin (原著), オズワルド マーティン (著), 金子 武久 (翻訳)
本題に入る前にベイズ統計学の一般的な例を解説します。
ベイズ統計モデルは古典的統計学では解けない課題を解く
次の図は，北海道大学 久保拓弥先生の「線形モデルの発展」です。
from IPython.display import Image
Image("./image_common/book_Davidson-Pilon_fig2.jpg")
線形モデルは単回帰でも重回帰でも最小二乗法で正規分布しか扱いません。
一般化線形モデル（GLM ; generalized linear model）は，確率分布とリンク関数を用いてモデリングします。一番有名なものにロジスティック回帰がありますが，二値分類に使います。
一般化線形混合モデル（GLMM ; generalized linear mixed model，混合効果モデルともいう）は，個体差，場所差を扱えるようにしたものです。
GLMMをさらに発展させたものが階層ベイズモデルです。ベイズ統計モデルはMCMCを使うことが不可分ではありませんが，階層ベイズモデルくらい複雑になるともはやMCMCを使うしかありません。
以上の4つは，AIエンジニア/データサイエンティストの必須スキルなのです。
この記事で紹介するモデルは階層ベイズモデルではなくMCMCを使った普通のベイズモデルです。階層ベイズモデルは，本サイトの別の記事で紹介します。
この記事では，MCMCを使わず，解析的に解くベイズモデルも紹介します。最初は，ベイズ統計モデルではなく，ベイズの定理（公式）を使った問題から紹介します。
ベイズ統計学にも段階がある
ベイズ定理の公式で病気の診断の問題を解く
病気の検査・診断法の問題をベイズ定理を使って解きましょう。
1. 実際に病気の人に対して正しく病気であるという診断になるのは98％です。
2. 実際に病気でない人に対して誤って病気であるという診断になるのは5％です。
3. 検査を受ける人の全体では，実際に病気である人と病気でない人の割合は3％，97％です。
4. ある人の検査結果が病気であるという診断のとき本当に病気である確率は？
ベイズ定理は次式です
$$ P\left(A|B\right) = \dfrac {P\left(B|A\right)P\left(A\right)}{P\left(B\right)} $$ $$ where P\left(B\right) = P\left(B|A\right)P\left(A\right) + P\left(B|\overline{A}\right)P\left(\overline{A}\right) $$
$P\left(A\right)$；事象Aが起こる確率
$P\left(B\right)$；事象Bが起こる確率
$P\left(A|B\right)$；事象Bの条件下で事象Aが起こる確率
$P\left(B|A\right)$；事象Aの条件下で事象Bが起こる確率
上の問題をベイズ定理に当てはめると，
事象A；実際に病気である
事象B；診断は病気である（Positive）
問題1.；$P\left(B|A\right)$
問題2.；$P\left(B|\overline{A}\right)$
問題3.；$P\left(A\right)$，$P\left(\overline{A}\right)$
問題4.；$P\left(A|B\right)$
from IPython.display import Image
Image("./image_common/book_Davidson-Pilon_fig3.jpg")
与えられた数字を Confusion Matrix ライクに書きます。
Confusion Matrix のマスの中の記号は下で説明しています。
横軸が予測（検査・診断の結果），縦軸が実際の病気の罹患になります。
赤字が上で与えられた数値であり，わきにベイズ定理の因子を書いてあります。
上のwhereの式の2つの項は次です。
TP = $P\left(B|A\right)P\left(A\right)$
FP = $P\left(B|\overline{A}\right)P\left(\overline{A}\right)$
ベイズ定理をConfusion Matrixのマスの中の記号を使って，古典的確率で解いてみましょう。直感的な式がベイズ定理の公式と同じになります。
解は，P(A|B) = TP / ( TP + FP ) = 0.377（Precision，陽性的中率とか言います）
（Jupyter Notebook の出力の）Confusion Matrixの座標の取り方を説明します。
Predicted 0 1 Actual 0［［TN FP］ 1 ［FN TP］］
昔のコンピュータスクリーンのように左上が原点です。横軸が Pridicted Label であり，右向きに 0, 1 となり，縦軸が Actual Label であり，下向きに 0, 1 となります。
「0」は Negative，「1」は Positive です。医療関係では伝統的に，菌や病原が検出されることを Positive 陽性と表現します。 Pridicted に対して使い，Actual に対しては使いません。
従って，左半分が Negative ，右半分が Positive となり，Pridicted と Actual が一致する左上と右下が True ，一致しない右上と左下が False となります。
これを覚えておくと，TN, FP, FN, TP の位置をマチガイません。
ベイズモデルの尤度と事前分布で事後分布を解析的に計算する
まだ，MCMCの確率的プログラミングを使いません。
ここでもう一度，ベイズ定理を見ます。
$$ P\left(A|B\right) = \dfrac {P\left(B|A\right)P\left(A\right)}{P\left(B\right)}\propto {P\left(B|A\right)P\left(A\right)} $$最後の項を新しく解釈します。
$P\left(B|A\right)$；尤度
$P\left(A\right)$；事前分布
$P\left(A|B\right)$；事後分布
例題：コイントス問題
コイントスの尤度の候補は，2項分布です。事前分布の候補は，ベータ分布です。尤度の2項分布に対してベータ分布は共役事前分布であり，事後分布もベータ分布となります。それで事後分布が解析的に計算できるのです。
正規分布はそれ自体の共役事前分布です。
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import seaborn as sns
plt.style.use('seaborn-darkgrid')
ベイズ定理の公式の記号を書き替えておきます。
$$ p\left(\theta|y\right) \propto {p\left(y|\theta\right)p\left(\theta\right)} $$尤度は2項分布
$$ p\left(y|\theta\right) = \dfrac {N!}{y!\left( N-y\right) !}\theta ^{y}\left( 1-\theta \right) ^{N-y} $$さまざまな2項分布を出力します。
n_params = [1, 2, 4]
p_params = [0.25, 0.5, 0.75]
x = np.arange(0, max(n_params)+1)
f, ax = plt.subplots(len(n_params), len(p_params), sharex=True, sharey=True)
for i in range(3):
for j in range(3):
n = n_params[i]
p = p_params[j]
y = stats.binom(n=n, p=p).pmf(x)
ax[i,j].vlines(x, 0, y, colors='b', lw=5)
ax[i,j].set_ylim(0,1)
ax[i,j].plot(0, 0, label="n = {:3.2f}\np = {:3.2f}".format(n, p), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$y$', fontsize=14)
ax[1,0].set_ylabel('$p(y|\\theta)$', fontsize=14)
ax[0,0].set_xticks(x)
plt.savefig('img103.png', dpi=300, figsize=(5.5, 5.5))
plt.figure()
事前分布はベータ分布
$$ p\left(\theta\right) = \dfrac {\Gamma \left( \alpha +\beta \right) }{\Gamma \left( \alpha \right) \Gamma \left( \beta \right) } \theta ^{\alpha -1}\left( 1-\theta \right) ^{\beta -1} =Beta\left( \alpha ,\beta \right) $$さまざまなベータ分布を出力します。
パラメータが [1, 1]のときが一様分布であり，無情報事前分布となります。
params = [0.5, 1, 2, 3]
x = np.linspace(0, 1, 100)
f, ax = plt.subplots(len(params), len(params), sharex=True, sharey=True)
for i in range(4):
for j in range(4):
a = params[i]
b = params[j]
y = stats.beta(a, b).pdf(x)
ax[i,j].plot(x, y)
ax[i,j].plot(0, 0, label="$\\alpha$ = {:3.2f}\n$\\beta$ = {:3.2f}".format(a, b), alpha=0)
ax[i,j].legend(fontsize=12)
ax[3,0].set_xlabel('$\\theta$', fontsize=14)
ax[0,0].set_ylabel('$p(\\theta)$', fontsize=14)
plt.savefig('img104.png', dpi=300, figsize=(5.5, 5.5))
plt.figure()
事後分布は，上の式に示したとおり，尤度と事前分布の掛け算になります。$\theta$に依存しないすべての項を省略します。
$$ p\left(\theta|y\right) \propto \theta ^{y}\left( 1-\theta \right) ^{N-y} \theta ^{\alpha -1}\left( 1-\theta \right) ^{\beta -1} $$ $$ p\left(\theta|y\right) \propto \theta ^{\alpha -1+y}\left( 1-\theta \right) ^{\beta -1+N-y} $$上の式をよく見ると，ベータ分布になります。このように，事前分布と事後分布が同じになることについて，尤度の2項分布に対してベータ分布は共役事前分布であると言います。
次のように整理します。
$$ \alpha _{pos}=\alpha _{pri}+y , \beta _{pos}=\beta _{pri}+N-y $$ $$ p\left(\theta|y\right) = Beta\left( \alpha _{pri}+y,\beta_{pri}+N-y\right) $$事後分布を計算し，グラフ化する
trialsはコインを投げた回数$N$です。dataは表が出た回数$y$です。beta_paramsは事前分布のベータ分布のパラメータを3種類用意してあります。theta_realはコインの偏りですが，図にバーを書き込むだけに使っていて，もちろん，事後分布の計算に無関係です。
$N$と$y$を順に大きくしてありますが，収束していく様子を図に示したかっただけで，本来は，1組だけ，例えば，$N=100$，$y=35$とするだけで事後分布は得られます。
theta_real = 0.35
trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
beta_params = [(1, 1), (0.5, 0.5), (20, 20)]
dist = stats.beta
x = np.linspace(0, 1, 100)
for idx, N in enumerate(trials):
if idx == 0:
plt.subplot(4, 3, 2)
else:
plt.subplot(4, 3, idx+3)
y = data[idx]
for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')):
p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)
plt.plot(x, p_theta_given_y, c)
plt.fill_between(x, 0, p_theta_given_y, color=c, alpha=0.6)
plt.axvline(theta_real, ymax=0.3, color='k')
plt.plot(0, 0, label="{:d} experiments\n{:d} heads".format(N, y), alpha=0)
plt.xlim(0,1)
plt.ylim(0,12)
plt.xlabel(r'$\theta$')
plt.legend()
plt.gca().axes.get_yaxis().set_visible(False)
plt.savefig('img105.png', dpi=300, figsize=(5.5, 5.5))
plt.figure()
一段目の図は事前分布を表しています。青い線は一様分布です。赤い線は一様分布に似ています。緑線は0.5の周辺に中心があり，そこに集中してます。黒いバーはコインの偏りの位置にあります。これは本来は未知であるものです。
2段目以降の図は事後分布の実験の進行を表しています。色によって収束のスピードが違っています。
このように，ベイズ統計の特定の問題に対しては，MCMCを使わなくて解が得られます。ただし，以上の2つ，つまり，ベイズ定理の公式問題，ベイズモデルを解析的に解く問題は，AIエンジニア/データサイエンティストの現場ではあまり使われないと思います。が，入社試験に出たりします。ベイズ定理の公式の意味は把握しておいた方が良いと思います。
ベイズ統計モデルをMCMC（Markov chain Monte Carlo）で解く
確率分布関数の積分を解析的に解けない場合に乱数シミュレーションを行う
上と同じで，コイントスの問題を続けます。
import pymc3 as pm
import numpy as np
import pandas as pd
from theano import shared
import scipy.stats as stats
import matplotlib.pyplot as plt
2項分布で1回だけ投げたときにはベルヌーイ分布になります。
np.random.seed(123)
n_experiments = 4
theta_real = 0.35 # unkwon value in a real experiment
data = stats.bernoulli.rvs(p=theta_real, size=n_experiments)
data
次にMCMCのモデリングをします。MCMCの典型的な書き方です。
withでモデル構築をします。as our_first_modeでモデルに名前を付けて，with our_first_model:でモデルの構築を続けます。
pm.sampleで事前分布を出力します。withのインデントでpm.sampleの前のコードはすべて事前分布です。少し前の式での尤度は，この場合は，入力データになります。上の式のベルヌーイ分布から出力されるデータが尤度となります。
stepはMCMCの種類です。
traceが事後分布のサンプルです。
with pm.Model() as our_first_model:
# a priori
theta = pm.Beta('theta', alpha=1, beta=1)
# likelihood
y = pm.Bernoulli('y', p=theta, observed=data)
#y = pm.Binomial('theta',n=n_experimentos, p=theta, observed=sum(datos))
start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(1000, step=step, start=start)
burnin = 0 # no burnin
chain = trace[burnin:]
# pm.plots.traceplot(chain, lines={'theta':theta_real});
pm.plots.traceplot(chain, lines={'theta':theta_real});
plt.savefig('B04958_02_04.png', dpi=300, figsize=(5.5, 5.5))
with our_first_model:
step = pm.Metropolis()
multi_trace = pm.sample(1000, step=step)
# multi_trace = pm.sample(1000, step=step, njos=4)
burnin = 0 # no burnin
multi_chain = multi_trace[burnin:]
pm.traceplot(multi_chain, lines={'theta':theta_real});
plt.savefig('B04958_02_06.png', dpi=300, figsize=(5.5, 5.5))
pm.gelman_rubin(multi_chain)
pm.forestplot(multi_chain, varnames=['theta']);
plt.savefig('B04958_02_07.png', dpi=300, figsize=(5.5, 5.5))
pm.summary(multi_chain)
モデリング対象テーマ「メッセージ数に変化はあるか？」
日々受信するテキストメッセージ（eメールのこと？）の数は，0件から数十件とバラバラですが，ある日を境に増えているようですが，本当にそうなのかベイズ統計モデルを使って解き明かします。
最後の図が毎日の受信数とMCMCで解いた受信数の期待値の図です。明らかに期待値の特異点があります。毎日の受信数のグラフから決して読み取れないと思われます。また，古典的な統計学（頻度主義統計学ともいう非ベイズ統計学）では解けない課題かと思われます。
毎日のメッセージ数のグラフ
%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
# count_data = np.loadtxt("data/txtdata.csv")
# print(count_data)
count_data = np.array(
[13, 24, 8, 24, 7, 35, 14, 11, 15, 11, 22, 22, 11, 57, 11, 19, 29, 6,
19, 12, 22, 12, 18, 72, 32, 9, 7, 13, 19, 23, 27, 20, 6, 17, 13, 10,
14, 6, 16, 15, 7, 2, 15, 15, 19, 70, 49, 7, 53, 22, 21, 31, 19, 11,
18, 20, 12, 35, 17, 23, 17, 4, 2, 31, 30, 13, 27, 0, 39, 37, 5, 14,
13, 22,])
figsize(12.5, 3.5)
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data);
$i$日目のメッセージ数を$C_i$とすると，次のように書ける。
$$ C_i \sim \text{Poisson}(\lambda) $$
ポアソン分布のパラメータ$\lambda$は期待値と分散に一致している。そこで上のグラフのメッセージ数がある日から増えていると仮定して，期待値を次のように表す。
$$ \lambda = \begin{cases} \lambda_1 & \text{if } t \lt \tau \cr \lambda_2 & \text{if } t \ge \tau \end{cases} $$
$\lambda _{1}$, $\lambda _{2}$の事前分布を次のように表す。
$$\lambda _{1}\sim Exp\left( \alpha \right)$$ $$\lambda _{2}\sim Exp\left( \alpha \right)$$
$\alpha$はハイパーパラメータと呼ばれる。指数分布の期待値は次の式になるので，これから$\alpha$を決めればよい。
$$\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha}$$
$\tau$は一様分布にする。
$$\tau \sim DiscreteUniform\left( 1,70\right)$$ $$\Rightarrow P\left( \tau = k\right) =\dfrac {1}{70}$$
モデル構築
次のコードは，with pm.Model()によりモデルを構築している。as modelによって，モデルの名称をmodelとしてして，つぎからは，with modelによりモデルの内容を追加している。
$\tau$, $\lambda_1$, $\lambda_2$は確率変数である。
import pymc3 as pm
import theano.tensor as tt
with pm.Model() as model:
alpha = 1.0/count_data.mean() # Recall count_data is the
# variable that holds our txt counts
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
with model:
idx = np.arange(n_count_data) # Index
lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
新しい関数を定義しているが，これが本当の$\lambda$となる確率変数である。
with model:
observation = pm.Poisson("obs", lambda_, observed=count_data)
上の式が尤度に相当するものです。ここで事前分布とデータが渡されます。
それより前の式はすべて事前分布に関するものです。
import scipy.stats as stats
with model:
start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(10000, tune=5000,step=step)
startは，初期値を事後分布のピークにすることを示しています。
stepは，MCMCの種類です。
traceは事後分布のサンプルです。
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
$\tau$, $\lambda_1$, $\lambda_2$の事後分布を表示する
当然のことながら事前分布にまったく似ていない
figsize(12.5, 10)
#histogram of the samples:
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_1$", color="#A60628", normed=True)
plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables
$\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_2$", color="#7A68A6", normed=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")
plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_count_data, alpha=1,
label=r"posterior of $\tau$",
color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data)-20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");
$\lambda_1$はおよそ18, $\lambda_2$はおよそ23である。
最後にメッセージ数の期待値を表示する。あくまでも確率分布であることを忘れてはならない。
figsize(12.5, 5)
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
# ix is a bool index of all tau samples corresponding to
# the switchpoint occurring prior to value of 'day'
ix = day < tau_samples
# Each posterior sample corresponds to a value for tau.
# for each day, that value of tau indicates whether we're "before"
# (in the lambda1 "regime") or
# "after" (in the lambda2 "regime") the switchpoint.
# by taking the posterior sample of lambda1/2 accordingly, we can average
# over all samples to get an expected value for lambda on that day.
# As explained, the "message count" random variable is Poisson distributed,
# and therefore lambda (the poisson parameter) is the expected value of
# "message count".
expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
+ lambda_2_samples[~ix].sum()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
label="observed texts per day")
plt.legend(loc="upper left");
print(expected_texts_per_day)