Up

平均値の差のベイズ推定

2つの関連あるサンプルの場合

 

条件1と条件2の2つの条件での値の対としてデータが得られている場合の、条件1と条件2の平均値の差をベイズ推定するスクリプトを書いてみた。ベイズ推定はPyMC3を使って行うが、PyMC3については、<岡本安晴いまさら聞けないPythonデータ分析丸善出版>で説明している。

いま、条件1でのデータX1と条件2でのデータ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]]

 

と与えられているとき、条件1と条件2の差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上で実行するとよい)と、分析結果が以下のように図示される。

 

図1

 

図1は、サンプリングのトレースである。サンプリングに特に問題は無いようである。このグラフは、スクリプト

az.plot_trace(trace)

によって描かれている。

 

図2

 

図2は、差(条件1―条件2)の平均値μの事後分布である。この図はスクリプト

az.plot_posterior(trace['μ'], ref_val = 0.0)

plt.title('μ', fontsize = 20)

によって描かれている。

94% HPDに0が含まれず、94% HPDは正の側に位置している。また、ref_valとして設定した値0の左側、すなわち負の値の領域の確率は2.3%となっている。差(条件1―条件2)の平均値μは正である、すなわち、条件1の平均値の方が条件2より大きいと考えてよさそうである。

 

図3

 

図3は、差(条件1―条件2)の分布のデータと事後予測のグラフである。データの方は分布の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個のμとσのペアをランダムに選んで確率密度関数を描いている(青の曲線群)。事後予測による分布関数とデータの分布関数の関係に問題は無さそうである。

 

なお、WindowsPyMC3を利用するときは、Jupyter Notebook上で実行するとよい。Jupyter Notebookは、Anacondaにインストールして使用できる。Anacondaのインストールは、このウェブサイトで説明している。インストールしたAnacondaの仮想環境にJupyter Notebookがインストールされていないときは、condaコマンドでインストールすればよい。

conda install jupyter

である。PyMC3Arvizとともに

conda install pymc3 arviz

とインストールすればよい。

 

 

リスト1 条件の差の平気インチのベイズ推定

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()

 

 

Up