平均値の差のベイズ推定
２つの関連あるサンプルの場合
条件１と条件２の２つの条件での値の対としてデータが得られている場合の、条件１と条件２の平均値の差をベイズ推定するスクリプトを書いてみた。ベイズ推定はPyMC3を使って行うが、PyMC3については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
いま、条件１でのデータX1と条件２でのデータX2が、対［X1, X2］の形で得られているとする。このとき、X1 – X2 の平均値の推定を行う。
データが
data = [[198, 198], [225, 193], [213,
191], [145, 178], [165, 161],
[181, 179], [183, 214], [217, 205], [226, 207], [187, 186],
[178, 165], [208, 204], [200, 199], [224, 208], [256, 204],
[224, 188], [136, 145], [188, 181], [233, 187], [217, 203]]
と与えられているとき、条件１と条件２の差diffを
diff = [v[0]
- v[1] for v in data]
と計算して、次のPyMC3のスクリプトでMCMCサンプリングを行う。
with pm.Model()
as model_rel:
μ = pm.Normal('μ',
mu = 0.0, sigma = 1000.0)
σ = pm.HalfNormal('σ',
sigma = 1000.0)
y = pm.Normal('y',
mu = μ, sigma = σ, observed = diff)
trace = pm.sample(3000)
プログラムの完全なスクリプトは、このサイトの後ろにリスト１として示した。
リスト１を実行する（Windowsで実行するときは、Jupyter Notebook上で実行するとよい）と、分析結果が以下のように図示される。
図１
図１は、サンプリングのトレースである。サンプリングに特に問題は無いようである。このグラフは、スクリプト
az.plot_trace(trace)
によって描かれている。
図２
図２は、差（条件１―条件２）の平均値μの事後分布である。この図はスクリプト
az.plot_posterior(trace['μ'], ref_val
= 0.0)
plt.title('μ', fontsize
= 20)
によって描かれている。
94% HPDに０が含まれず、94% HPDは正の側に位置している。また、ref_valとして設定した値０の左側、すなわち負の値の領域の確率は2.3％となっている。差（条件１―条件２）の平均値μは正である、すなわち、条件１の平均値の方が条件２より大きいと考えてよさそうである。
図３
図３は、差（条件１―条件２）の分布のデータと事後予測のグラフである。データの方は分布のKDE（カーネル密度推定）によるプロット（赤の曲線）である。スクリプトは
smpl_indx = np.random.choice(range(len(trace['μ'])), size = 100, replace = False)
x = np.linspace(-100,
100, 1000)
plt.title('Posterior Predictions', fontsize = 20)
for i in smpl_indx:
mu = trace['μ'][i]
sgm
= trace['σ'][i]
y = ss.norm(loc = mu, scale = sgm).pdf(x)
plt.plot(x,
y, 'b-', alpha = 0.1)
sb.kdeplot(diff, c = 'r', label = 'data')
である。事後予測は、MCMCサンプルから100個のμとσのペアをランダムに選んで確率密度関数を描いている（青の曲線群）。事後予測による分布関数とデータの分布関数の関係に問題は無さそうである。
なお、WindowsでPyMC3を利用するときは、Jupyter Notebook上で実行するとよい。Jupyter Notebookは、Anacondaにインストールして使用できる。Anacondaのインストールは、このウェブサイトで説明している。インストールしたAnacondaの仮想環境にJupyter Notebookがインストールされていないときは、condaコマンドでインストールすればよい。
conda install jupyter
である。PyMC3もArvizとともに
conda install pymc3 arviz
とインストールすればよい。
import pymc3 as pm
import arviz
as az
import numpy
as np
import matplotlib.pyplot
as plt
import scipy.stats
as ss
import seaborn
as sb
data = [[198, 198], [225, 193], [213,
191], [145, 178], [165, 161],
[181, 179], [183, 214], [217, 205], [226, 207], [187, 186],
[178, 165], [208, 204], [200, 199], [224, 208], [256, 204],
[224, 188], [136, 145], [188, 181], [233, 187], [217, 203]]
diff = [v[0] - v[1] for v in data]
with pm.Model()
as model_rel:
μ = pm.Normal('μ',
mu = 0.0, sigma = 1000.0)
σ = pm.HalfNormal('σ',
sigma = 1000.0)
y = pm.Normal('y',
mu = μ, sigma = σ, observed = diff)
trace = pm.sample(3000)
print(az.summary(trace))
az.plot_trace(trace)
plt.savefig('Fig_trace.png')
plt.show()
az.plot_posterior(trace['μ'], ref_val
= 0.0)
plt.title('μ', fontsize
= 20)
plt.savefig('Fig_mu.png')
plt.show()
smpl_indx = np.random.choice(range(len(trace['μ'])), size = 100, replace = False)
x = np.linspace(-100,
100, 1000)
plt.title('Posterior Predictions', fontsize = 20)
for i in smpl_indx:
mu = trace['μ'][i]
sgm
= trace['σ'][i]
y = ss.norm(loc = mu, scale = sgm).pdf(x)
plt.plot(x,
y, 'b-', alpha = 0.1)
sb.kdeplot(diff, c = 'r', label = 'data')
plt.legend()
plt.savefig('Fig_Pred.png')
plt.show()