独立なサンプルの平均値の差のベイズ分析（CmdStanPy対応）

２つの独立なデータの正規分布によるベイズ分析

CmdStanPy対応

正規分布は次式で与えられる（Gelman et al., 2014）。

独立な２のデータ、XとY、が与えられたとき、それらが正規分布に従うとしてベイズ分析を行うStanスクリプトをリスト１のように用意した、パラメータは、データXに関するものは添字１を付け、データYに関するものは添字２を付けている。ファイルはbanormalcspfiles.zipにまとめた。

リスト１　正規分布に従う２つの独立なデータのStanスクリプト（ba_normal.stan）

data {

int n1;

array[n1] real x;

int n2;

array[n2] real y;

real min1;

real max1;

real min2;

real max2;

}

parameters {

real<lower=0, upper=1> vmu1;

real<lower=0, upper=1> vsgm1;

real<lower=0, upper=1> vmu2;

real<lower=0, upper=1> vsgm2;

}

transformed parameters {

real mu1;

real<lower=0> sgm1;

real mu2;

real<lower=0> sgm2;

Parameters rescaled

mu1 = min1 + (max1 - min1) * vmu1;

sgm1 = 0.0001 + (max1 - min1) * vsgm1 / 2.0;

mu2 = min2 + (max2 - min2) * vmu2;

sgm2 = 0.0001 + (max2 - min2) * vsgm2 / 2.0;

}

model {

vmu1 ~ beta(1, 1);

vsgm1 ~ beta(1, 1);

vmu2 ~ beta(1, 1);

vsgm2 ~ beta(1, 1);

x ~ normal(mu1, sgm1);

y ~ normal(mu2, sgm2);

}

乱数の生成の変域を区間に制約して生成している。すなわち

vmu1 ~ beta(1, 1);

vsgm1 ~ beta(1, 1);

vmu2 ~ beta(1, 1);

vsgm2 ~ beta(1, 1);

である。これらの乱数を、データに合わせた変域の乱数に変換して、モデルをデータに適用している。すなわち、

mu1 = min1 + (max1 - min1) * vmu1;

sgm1 = 0.0001 + (max1 - min1) * vsgm1 / 2.0;

mu2 = min2 + (max2 - min2) * vmu2;

sgm2 = 0.0001 + (max2 - min2) * vsgm2 / 2.0;

と変換したパラメータにより、データとモデルの関係が次のように設定されている。

x ~ normal(mu1, sgm1);

y ~ normal(mu2, sgm2);

リスト１のStanスクリプトをCmdStanPyを利用して用い、ベイズ分析を行うPythonスクリプトをリスト２のように用意した。CmdStanPyの簡単な準備方法の説明を、このウェブサイトで行った。

リスト２　リスト１のStanスクリプトをCmdStanPyを利用して用いるPythonスクリプト（ba_normal.py）

from cmdstanpy import CmdStanModel

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import arviz as az

import pandas as pd

import seaborn as sb

inf_nm = input('Input data file(*.xlsx) = ')

dframe = pd.read_excel(inf_nm) # データの読み込み

df_keys = dframe.keys()

x_name = df_keys[1]

y_name = df_keys[2]

print(x_name, y_name)

x = dframe[x_name].values

x = x[np.isnan(x)==False] # 空セルデータの削除

print('X =\n', x)

y = dframe[y_name].values

y = y[np.isnan(y)==False] # 空セルデータの削除

print('Y =\n', y)

mean1 = np.mean(x)

sd1 = np.std(x)

print('mean1 =', mean1, ' sd1 =', sd1)

mean2 = np.mean(y)

sd2 = np.std(y)

print('mean2 =', mean2, ' sd2 =', sd2)

min1 = np.min(x)

max1 = np.max(x)

min2 = np.min(y)

max2 = np.max(y)

print(min1,' ', max1)

print(min2,' ', max2)

Data = {'n1':len(x), 'x':x, 'min1':min1, 'max1':max1,

'n2':len(y), 'y':y, 'min2':min2, 'max2':max2}

model = CmdStanModel(stan_file="ba_normal.stan") # コンパイル

fit = model.sample(data=Data) # サンプリング

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

# トレースプロット

az.plot_trace(inf_data, var_names=['mu1','sgm1','mu2','sgm2'])

plt.tight_layout()

plt.show()

pdf_fit = fit.draws_pd() # Pandas DataFrameに変換

"""

事後分布の中央値を点推定値とする

"""

mu1 = np.median(pdf_fit['mu1'])

sgm1 = np.median(pdf_fit['sgm1'])

mu2 = np.median(pdf_fit['mu2'])

sgm2 = np.median(pdf_fit['sgm2'])

"""

パラメータmu1とmu2の事後分布

"""

plt.title('Posterior distributions' +f'\nmu1={mu1:.2f}, mu2={mu2:.2f}')

sb.kdeplot(pdf_fit['mu1'], label=x_name)

sb.kdeplot(pdf_fit['mu2'], label=y_name)

plt.xlabel('mu', fontsize=14)

plt.legend()

plt.show()

"""

パラメータmu1とmu2の差の事後分布と事後確率

"""

diff_mu1_mu2 = pdf_fit['mu1'] - pdf_fit['mu2']

diff_mu_med = np.median(diff_mu1_mu2)

pos_idx = diff_mu1_mu2 > 0

neg_idx = diff_mu1_mu2 < 0

prop_pos = np.mean(pos_idx) # P(mu1-mu2>0)

prop_neg = np.mean(neg_idx) # P(mu1-mu2<0)

sb.kdeplot(diff_mu1_mu2)

plt.xlabel('mu1-mu2', fontsize=14)

plt.title(f'P(mu1>mu2)={prop_pos:.5f}, P(mu1<mu2)={prop_neg:.5f}' +

f' med.={diff_mu_med:.2f}')

plt.tight_layout()

plt.show()

"""

パラメータsgm1とsgm2の事後分布

"""

plt.title('Posterior distributions'

f'\nsgm1={sgm1:.1f}, sgm2={sgm2:.2f}')

sb.kdeplot(pdf_fit['sgm1'], label=x_name)

sb.kdeplot(pdf_fit['sgm2'], label=y_name)

plt.xlabel('sgm', fontsize=14)

plt.legend()

plt.show()

"""

パラメータsgm1とsgm2の差の事後分布と事後確率

"""

diff_sgm1_sgm2 = pdf_fit['sgm1'] - pdf_fit['sgm2']

diff_sgm_med = np.median(diff_sgm1_sgm2)

prop_pos = np.mean(diff_sgm1_sgm2 > 0) # P(mu1-mu2>0)

prop_neg = np.mean(diff_sgm1_sgm2 < 0) # P(mu1-mu2<0)

sb.kdeplot(diff_sgm1_sgm2)

plt.xlabel('sgm1-sgm2', fontsize=14)

plt.title(f'P(sgm1>sgm2)={prop_pos:.5f}, P(sgm1<sgm2)={prop_neg:.5f}' +

f' med.={diff_sgm_med:.2f}')

plt.tight_layout()

plt.show()

"""

　　　　　　ヒストグラムと密度関数

"""

plt.hist(x, density=True, color='b', alpha=0.3, label=x_name)

plt.hist(y, density=True, color='r', alpha=0.3, label=y_name)

minxy = min1 if min1 < min2 else min2

maxxy = max1 if max1 > max2 else max2

x_coord = np.linspace(minxy, maxxy, 1000)

y_coord1 = ss.norm.pdf(x_coord, loc=mu1, scale=sgm1)

plt.plot(x_coord, y_coord1, color='b', label=x_name)

y_coord2 = ss.norm.pdf(x_coord, loc=mu2, scale=sgm2)

plt.plot(x_coord, y_coord2, color='r', label=y_name)

plt.title(f'{x_name}: mu1={mu1:.1f}, sgm1={sgm1:.1f}' +

f'\n{y_name}: mu2={mu2:.1f}, sgm2={sgm2:.1f}')

plt.legend()

plt.show()

Stanスクリプトのコンパイルを、次のコードで行っている。

model = CmdStanModel(stan_file="ba_normal.stan")

コンパイルされたmodelによるMCMCサンプリングを次のコードで実行している。

fit = model.sample(data=Data)

サンプルの統計を、次のコードで表示している。

print(fit.summary())

MCMCサンプリングのトレースを描くために、MCMCサンプリングの結果をarviz（az）のInferenceData型に変換している。

inf_data = az.from_cmdstanpy(fit)

変換されたデータinf_dataを用いて、次のスクリプトによりトレース図を描いている。

az.plot_trace(inf_data, var_names=['mu1','sgm1','mu2','sgm2'])

plt.tight_layout()

plt.show()

MCMCサンプルを取り出すために、次のコードでサンプルfitをPandasのDataFrame型に変換している。

pdf_fit = fit.draws_pd()

各パラメータのサンプルは、例えば次のスクリプトのようにPandasのDataFrameから取り出せる。

plt.title('Posterior distributions' +f'\nmu1={mu1:.2f}, mu2={mu2:.2f}')

sb.kdeplot(pdf_fit['mu1'], label=x_name)

sb.kdeplot(pdf_fit['mu2'], label=y_name)

plt.xlabel('mu', fontsize=14)

plt.legend()

plt.show()

上のスクリプトは、パラメータmu1とmu2のサンプルを取り出して事後分布のグラフを描くものである。

リスト２のスクリプトを、MiniConda/Ubuntu24/WSL/Windowsにおいて実行した。CmdStanPyの簡単な準備方法は、このウェブサイトで説明した。

リスト１とリスト2のファイル、および入力データファイルを同じフォルダにおいて、リスト２のスクリプトを実行する。

次のコマンド

(stan) *****/banormalcspfiles$ python ba_normal.py

を実行すると、入力データファイル名を聞いてくる。次の例では、data.xlsxを設定している。

Input data file(*.xlsx) = data.xlsx

入力データファイルは、図１a、図１ｂの形式で、Excelファイルとして用意する。

図１a

図１ｂ

ファイルの形式は、1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用であり、分析では用いられない。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。

入力データファイル名を設定してEnterキーを押すと、データが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、時間が掛かるが、コンパイル後、ＭＣＭＣサンプリングが実行される。

MCMCサンプリングが終了すると、トレース図が表示される（図２）。