ベイズ分析の計算入門
岡本安晴
2019.7、2025.3
ベイズ分析では、データとモデル（パラメータ）との組み合わせに対して確率（同時確率）を考え、その確率を次式で
で与える。確率
は、データを固定してモデルの関数と見たときは尤度関数と呼ばれるものであることに注意。また、モデルは、あらゆるモデルすべてを考えるのではなく、ある枠内のものが対象とされる。例えば、N回の硬貨投げでは二項分布というモデルの範囲で考え、モデルは表の出る確率
で表される。すなわち、モデルは想定されている確率分布のパラメータで表される。このとき、式（１）は
と書くことができる。
統計モデルとして簡単である二項分布を例として、ベイズ分析の入門解説を試みる。説明においてプログラム（スクリプト）が示されることがあるが、スクリプトファイルは、まとめて圧縮ファイルbysn2files.zipとした。
事前分布と事後分布
二項分布は、硬貨投げの結果が表であるか裏であるかというベルヌーイ試行を
試行繰り返したときの表の出る回数の確率を与えるものである。各試行における表の出る確率をとおくと、表の出る回数が
である確率は、次式
で与えられる。二項分布を用いた場合、モデルはパラメータによって決まると考え、式（１）は
と書くことができる。このとき、表の出た回数がであるというデータが与えられたという条件の下でのの確率分布（事後分布と呼ばれている）は
である。表の出た回数がであるというデータが与えられたという条件の下では回数は固定されていて、事後確率分布
はの関数として考えることができる。このとき、
は定数となるので、式（３）は
と書くことができる。
確率分布
は事前確率と呼ばれているが、特に制約を置かないときは無情報事前分布として一様分布を設定することができる。すなわち、
とおく。このとき、式（４）は
と置くことができる。
いま、
回の試行において
回表が出たとすると、式（５）は
となる。事後分布（６）のグラフは、次のスクリプトで描くことができる（リスト１）。
def N10k7(p):
return (p ** 7) * ((1 - p)
** 3)
p_values = np.arange(0.0,
1.0, 0.01)
p_post = N10k7(p_values)
plt.plot(p_values, p_post)
plt.yticks([])
plt.show()
上のスクリプトを実行すると図１のグラフを得る。
事後予測
ベイズ分析では、事後分布に基づいて分析結果の検討が行われる。
いま、得られた事後分布
が、データである表の出た回数
とどの程度の整合性があるか調べる。事後分布からをランダムに選び、そのの下で二項分布
（式（２））に従ってランダムにのサンプルをとる。このとのサンプリングを繰り返して得られるの分布をデータの値と比較する。このサンプリングのシミュレーションを行うスクリプトを、事後分布がベータ分布
であることに注意して、以下のように用意した（リスト２）。
a = 7 + 1
b = 10 - 7 + 1
N = 10
k_values = []
p_samples = []
for t in range(10000):
p = scipy.stats.beta.rvs(a,
b)
k = scipy.stats.binom.rvs(N,
p)
p_samples.append(p)
k_values.append(k)
リストp_samplesにのサンプル値が、リストk_valuesにのサンプル値が格納されている。サンプルされたのヒストグラムとKDEプロットを、比較のためのベータ分布のグラフとともに描くために以下のスクリプトを用意した。
p_values = np.arange(0.0,
1.0, 0.01)
p_beta = scipy.stats.beta.pdf(p_values, a, b)
plt.plot(p_values, p_beta, linestyle = '--', color =
'k',
linewidth = 3, label = 'Beta(7+1, 3+1)')
plt.hist(p_samples,
bins = 20, density = True, color = 'b', alpha = 0.5)
seaborn.kdeplot(p_samples,
linewidth = 3, color = 'r',label = 'KDE plot')
plt.yticks([])
plt.xlabel('p', fontsize
= 14)
plt.legend(loc =
'upper left', fontsize = 14)
plt.title('Sample Distribution of p', fontsize = 16)
plt.show()
上のスクリプトを実行すると図２のグラフを得る。
サンプル値のヒストグラムとKDEプロット（赤の曲線）が、事後分布であるベータ分布の破線のグラフとよく一致していることが分かる。サンプル値の分布は問題ないようである。
サンプル値の分布の描画は、以下のスクリプトで行った。
plt.hist(k_values,
bins = np.arange(-0.25, 10.75, 0.5), density = True)
plt.vlines(7, 0, 0.4, color = 'r', linewidth =
3, label = 'k = 7')
plt.xlabel('k', fontsize
= 16)
plt.legend(loc =
'upper left', fontsize = 16)
plt.title('Posterior predictive distribution of
k', fontsize = 16)
plt.show()
実行すると図３のグラフを得る。データ値の位置が、赤の縦線で示されている。の事後予測分布は、このデータ値の位置を中心にして分布しており、問題ないと思われる。
追加データの分析
上の試行の二項分布における表の出た試行数が試行であるというデータを得た後で、さらにまた同じ硬貨を用いて追加の試行を行ったところ表の出た試行数が
試行であったとする。このデータを先の表の出た試行数が試行であったデータの分析結果に基づいて分析するため、先の分析で得られた事後分布（６）を今回の分析における事前分布としてベイズ分析を行う。
ベイズ分析の式
において、事前分布を先の分析の事後分布（６）
とおき、今回のデータに対する確率（尤度関数）を
とおく。このとき、先の分析結果に基づく、すなわち先の分析結果での事後分布を今回の分析の事前分布とするベイズ分析の事後分布
は
となる。
定数（積分して１になるように調整するためのもの）を無視すると、分布
はベータ分布
であることに注意して、事前分布と事後分布のグラフを描くスクリプトは以下のようになる（リスト３）。
p_values = np.arange(0.0,
1.0, 0.01)
p0_values = scipy.stats.beta.pdf(p_values, 7+1, 3+1)
p1_values = scipy.stats.beta.pdf(p_values, 10+1, 10+1)
plt.plot(p_values,
p0_values, label = 'Prior')
plt.plot(p_values,
p1_values, label = 'Posterior')
plt.yticks([])
plt.xlabel('p', fontsize
= 16)
plt.legend()
plt.show()
実行すると図４を得る。
先のデータの分析結果の事後分布が、新しく追加されたデータの分析において事前分布として用いられ、その事前分布（青い曲線）が新しく追加されたデータの分析結果である事後分布（橙色の曲線）として更新されている様子が表されている。
ベイズ分析では、先の分析結果を事前分布に反映させて、新しく得られたデータの分析を容易に行うことができる。
ベイズ分析は、事後分布がMonte Carlo Markov chain
(MCMC) というシミュレーションで推測する方法が開発され、その実用性が飛躍的に高まり、従来の最尤法では分析が難しいモデルもMCMCの利用によってベイズ分析が行われるようなった。MCMCのソフトは無料のものもいろいろ公開されているが、PyStan（Python用のStan）あるいはPyMCについては拙著「いまさら聞けないPythonでデータ分析」の第３部「ベイズ分析」（pp. 147-236）で説明した。
引用文献
岡本安晴（2019）．いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――．丸善出版．
import matplotlib.pyplot
as plt
import numpy
as np
def N10k7(p):
return (p ** 7) * ((1 - p)
** 3)
p_values = np.arange(0.0,
1.0, 0.01)
p_post = N10k7(p_values)
plt.plot(p_values, p_post)
plt.yticks([])
plt.show()
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.stats
import seaborn
a = 7 + 1
b = 10 - 7 + 1
N = 10
k_values = []
p_samples = []
for t in range(10000):
p = scipy.stats.beta.rvs(a,
b)
k = scipy.stats.binom.rvs(N,
p)
p_samples.append(p)
k_values.append(k)
p_values = np.arange(0.0,
1.0, 0.01)
p_beta = scipy.stats.beta.pdf(p_values, a, b)
plt.plot(p_values, p_beta, linestyle = '--', color =
'k',
linewidth = 3, label = 'Beta(7+1, 3+1)')
plt.hist(p_samples,
bins = 20, density = True, color = 'b', alpha = 0.5)
seaborn.kdeplot(p_samples,
linewidth = 3, color = 'r',label = 'KDE plot')
plt.yticks([])
plt.xlabel('p', fontsize
= 14)
plt.legend(loc =
'upper left', fontsize = 14)
plt.title('Sample Distribution of p', fontsize = 16)
plt.show()
plt.hist(k_values,
bins = np.arange(-0.25, 10.75, 0.5), density = True)
plt.vlines(7, 0, 0.4, color = 'r', linewidth =
3, label = 'k = 7')
plt.xlabel('k', fontsize
= 16)
plt.legend(loc =
'upper left', fontsize = 16)
plt.title('Posterior predictive distribution of
k', fontsize = 16)
plt.show()
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.stats
p_values = np.arange(0.0,
1.0, 0.01)
p0_values = scipy.stats.beta.pdf(p_values, 7+1, 3+1)
p1_values = scipy.stats.beta.pdf(p_values, 10+1, 10+1)
plt.plot(p_values,
p0_values, label = 'Prior')
plt.plot(p_values,
p1_values, label = 'Posterior')
plt.yticks([])
plt.xlabel('p', fontsize
= 16)
plt.legend()
plt.show()