平均値の差のベイズ推定
２つの独立なサンプルの場合
２つの独立なサンプルの平均値の差の事後分布を求めるスクリプトを書いてみた。
データは
data1 = [193, 204, 163, 189, 182, 222,
178, 142, 198, 197,
201, 192, 224, 220, 155, 189, 188, 215, 195, 169]
と
data2 = [223, 238, 217, 210, 226, 206,
194, 164, 200, 187,
232, 261, 237, 208, 234, 212, 227]
であるとき、PyMC3を用いて事後分布を求めるスクリプトをリスト１のように作成した。スクリプトにおいてμなどのギリシャ文字が使われているが、これはJupyter Notebookの場合、例えば、「\mu」と打ち込んでから「Tabキー」を押すと、「\mu」は「μ」に変換される。
すなわち、
\mu + Tabキー ―＞ μ
\sigma + Tabキー ―＞ σ
である。
リスト１ ２つの独立なサンプルの平均値の事後分布
import pymc3 as pm
import arviz
as az
import matplotlib.pyplot
as plt
import seaborn
as sb
data1 = [193, 204, 163, 189, 182, 222,
178, 142, 198, 197,
201, 192, 224, 220, 155, 189, 188, 215, 195, 169]
data2 = [223, 238, 217, 210, 226, 206,
194, 164, 200, 187,
232, 261, 237, 208, 234, 212, 227]
with pm.Model()
as two_samples_ind:
μ1 = pm.Normal('μ1',
mu = 0.0, sigma = 1000.0)
δ = pm.Normal('δ',
mu = 0.0, sigma = 1000.0)
μ2 = pm.model.Deterministic('μ2',
μ1 + δ)
σ = pm.HalfNormal('σ',
sigma = 1000.0)
y1 = pm.Normal('y1',
mu = μ1, sigma = σ, observed = data1)
y2 = pm.Normal('y2',
mu = μ2, sigma = σ, observed = data2)
trace = pm.sample(3000)
print(az.summary(trace))
az.plot_trace(trace)
plt.savefig('FigTrace.png')
plt.show()
az.plot_posterior(trace)
plt.savefig('FigPost.png')
plt.show()
sb.kdeplot(trace['μ1'], label = 'μ1')
sb.kdeplot(trace['μ2'], label = 'μ2')
plt.yticks([])
plt.legend()
plt.savefig('FigDistr.png')
plt.show()
リスト１を実行する（Windowsの場合は、Jupyter Notebook上で実行するとよい）と、以下の結果を得る。
図１
図２
図３
WindowsでPyMC3を利用するときは、Jupyter Notebook上で実行するとよい。Jupyter Notebookは、Anacondaにインストールして使用できる。Anacondaのインストールは、このウェブサイトで説明している。インストールしたAnacondaの仮想環境にJupyter Notebookがインストールされていないときは、condaコマンドでインストールすればよい。
conda install jupyter
である。PyMC3もArvizとともに
conda install pymc3 arviz
とインストールすればよい。
PyMC3については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。