Up

２つの独立なデータの比較

ｔ分布に基づくベイズ分析（CmdStanPy対応）

CmdStanPyの準備の簡単な説明を、このウェブサイトで行った。

正規分布の観点から見ると外れ値と見做される値も、ｔ分布のように裾の厚い分布においては他の値と同様に扱われる（Martin, 2018）。

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

は、自由度（degrees of freedom, df）と呼ばれることもあるパラメータである。

パラメータの幾つかの値に対するｔ分布と正規分布のグラフを図１に示す。

図１

パラメータの値が小さいとき裾の厚い分布であるが、の値が大きくなると正規分布に近付いている。

母集団分布としてｔ分布を設定したときのベイズ分析のStanスクリプトをリスト１のように用意した。リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２に示すものを用意した。ファイルはba_t_distri_files.zipにまとめた。自由に利用されたい。

リスト１のスクリプトファイル（t_distri.stan）、リスト２のスクリプトファイル（anal_t_distri_means.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされている環境において、次のコマンドを実行する。CmdStanPyの準備の簡単な説明を、このウェブサイトに用意した。

(stan) *****/ba_t_distri_files$ python anal_t_distri_means.py

上のコマンドを実行すると、入力データファイル名と出力ファイル名の設定が求められる。

(stan) *****/ba_t_distri_files$ python anal_t_distri_means.py

Input data file(*.csv) = sampleN100.csv

Output file(*.txt) = Results.txt

入力データファイルは、図２の形式のCSVファイルとして用意する。

・

・

・

図２　入力データファイル例sampleN100.csv

１行目に変数名を書き、２行目からデータ値を書く。１列目はデータの識別用で、プログラムでは用いられない。２列目と３列目に分析データ値を書くが、データ数が異なるときは少ない方のデータのセルは空にしておけばよい。

上の例では、図２のファイル名が設定されている。

出力アイルは、任意のテキストファイル名（*.txt）を設定する。

入力データファイル名と出力ファイル名を設定すると、データが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリング後、トレース図が表示される（図３）。特に問題はなさそうである。

図３

図３のWindowを閉じると、図４のグラフが表示される。

図４

平均値を表すパラメータmu1とmu2の事後分布のグラフと代表値が表示されている。分布の代表値は中央値である（Gelman et al., 2021）。

図４のWindowを閉じると、図５のグラフが表示される。

図５

タイトル部にmu1>mu2の事後確率が1.00000であること、mu1<mu2の事後確率が0.00000であることが表示されている。ほぼ確実にmu1の方がmu2より大きいと考えられる。確率が1.00000とか0.00000いうのは、事後サンプルに基づく比率を表していることによる。

図５のWindowを閉じると、パラメータsgmの事後分布が表示される（図６）

図６

図６Windowを閉じると、パラメータnuの事後分布が表示される（図７）。

図７

KDEによる推定分布では負の領域にも分布が正の値をとっている領域があるが、ヒストグラムではnuの分布は正の領域に限られていることが分かる。

図７のWindowを閉じると、図８のグラフが表示される。

図８

データのヒストグラムと、事後分布の中央値をパラメータとするｔ分布のグラフが描かれている。Data1のヒストグラムにおいて、左の裾に離れた値がある。これは正規分布を基準にすると外れ値の看做されうるが、ｔ分布では特別扱いする必要のない値であろう。

図８のWindowを閉じると、実行終了である。

端末には、以下のようなメッセージが表示されている。

Results.txt was saved.

(stan) *****/ba_t_distri_files$

出力ファイルResults.txtを開くと、以下のようである。

Input file = {inf_nm}

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -1040.090000 0.053335 1.895290 ... 2129.71 5128.12 1.002630

vmu1 0.682837 0.000208 0.012127 ... 2233.07 13653.20 1.001800

vsgm1 0.215601 0.000560 0.026806 ... 1993.94 9609.74 1.001300

vmu2 0.557414 0.000279 0.015715 ... 2782.38 12859.60 0.999727

vsgm2 0.272808 0.001021 0.042583 ... 2759.78 7290.02 1.000500

vnu1 0.089253 0.003547 0.098515 ... 1154.64 8313.62 1.001750

vnu2 0.160886 0.006110 0.215311 ... 1638.22 6208.54 1.001280

mu1 599.191000 0.122397 7.142490 ... 2233.07 13653.30 1.001790

sgm1 63.494500 0.164981 7.894400 ... 1993.94 9609.74 1.001300

mu2 497.561000 0.150923 8.501900 ... 2782.38 12859.60 0.999727

sgm2 73.794700 0.276076 11.518600 ... 2759.78 7290.02 1.000500

nu1 8.925400 0.354720 9.851480 ... 1154.64 8313.62 1.001750

nu2 16.088700 0.611046 21.531100 ... 1638.22 6208.54 1.001280

[13 rows x 11 columns]

Summary:

mean sd hdi_3% hdi_97% ... mcse_sd ess_bulk ess_tail r_hat

mu1 599.191 7.142 585.949 612.203 ... 0.109 3400.0 2233.0 1.0

mu2 497.561 8.502 480.975 512.747 ... 0.140 3202.0 2782.0 1.0

sgm1 63.495 7.894 47.863 77.552 ... 0.139 2393.0 1994.0 1.0

sgm2 73.795 11.519 53.611 95.485 ... 0.150 1815.0 2760.0 1.0

nu1 8.925 9.851 1.838 20.723 ... 1.109 2070.0 1155.0 1.0

nu2 16.089 21.531 1.546 67.382 ... 0.695 1546.0 1638.0 1.0

[6 rows x 9 columns]

ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。

リスト１　ｔ分布のStanスクリプト（t_distri.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;

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

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

}

transformed parameters {

real mu1;

real<lower=0> sgm1;

real mu2;

real<lower=0> sgm2;

real nu1;

real nu2;

/*

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;

nu1 = 0.0001 + vnu1*100;

nu2 = 0.0001 + vnu2*100;

}

model {

vmu1 ~ beta(1, 1);

vsgm1 ~ beta(1, 1);

vmu2 ~ beta(1, 1);

vsgm2 ~ beta(1, 1);

vnu1 ~ beta(1, 1);

vnu2 ~ beta(1, 1);

x ~ student_t(nu1, mu1, sgm1);

y ~ student_t(nu2, mu2, sgm2);

}

リスト２　リスト１を用いてベイズ分析を行うPythonスクリプト（anal_t_distri_means.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import seaborn as sb

import scipy.stats as ss

import arviz as az

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

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

fout_nm = input('Output file(*.txt) = ')

fout = open(fout_nm, 'w') #　出力用ファイルの用意

fout.write('Input file = {inf_nm}\n')

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="t_distri.stan")

# サンプリング

fit = model.sample(data=Data)

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

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

az.plot_trace(i_data, var_names=['mu1','mu2','sgm1','sgm2','nu1','nu2'],

figsize=(12,8))

plt.tight_layout()

plt.show()

print(az.summary(i_data))

fout.write('\nSummary:\n')

fout.write("{0}\n".format(az.summary(i_data,

var_names = ['mu1','mu2','sgm1','sgm2','nu1','nu2'])))

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

"""

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

"""

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

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

nu1 = np.median(pd_fit_frame['nu1'])

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

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

nu2 = np.median(pd_fit_frame['nu2'])

"""

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

"""

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

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

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

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

plt.legend()

plt.show()

"""

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

"""

diff_mu1_mu2 = pd_fit_frame['mu1'] - pd_fit_frame['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}')

plt.tight_layout()

plt.show()

"""

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

"""

plt.title('Posterior distributions'

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

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

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

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

plt.legend()

plt.show()

"""

パラメータnu1とnu2の事後分布

"""

plt.title('Posterior distributions' +

f'\nnu1={nu1:.2f}, nu2={nu2:.2f}')

plt.hist(pd_fit_frame['nu1'], color='g', alpha=0.5,

density=True, label=x_name)

plt.hist(pd_fit_frame['nu2'], color='orange', alpha=0.5,

density=True, label=y_name)

sb.kdeplot(pd_fit_frame['nu1'], color='g', label=x_name)

sb.kdeplot(pd_fit_frame['nu2'], color='orange', label=y_name)

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

plt.legend()

plt.show()

"""

データのヒストグラムとｔ分布のグラフ

"""

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

plt.hist(y, density=True, color='orange', alpha=0.5, label=y_name)

"""

グラフの横軸の範囲

"""

min_xy = min1 if min1 < min2 else min2

max_xy = max1 if max1 > max2 else max2

coord_xy = np.linspace(min_xy, max_xy, 1000)

density1 = ss.t.pdf(coord_xy, nu1, loc=mu1, scale=sgm1) # t分布

density2 = ss.t.pdf(coord_xy, nu2, loc=mu2, scale=sgm2) # t分布

plt.plot(coord_xy, density1, c='g', label=x_name)

plt.plot(coord_xy, density2, c='r', label=y_name)

plt.title(fr'{x_name}: $\mu$={mu1:.2f}, $\sigma$={sgm1:.2f}, $\nu$={nu1:.2f}' +

'\n' +

fr'{y_name}: $\mu$={mu2:.2f}, $\sigma$={sgm2:.2f}, $\nu$={nu2:.2f}',

fontsize=14)

plt.legend()

plt.tight_layout()

plt.show()

fout.close()

print(fout_nm, 'was saved.')

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University.

Martin, O. (2018). Bayesian analysis with Python, 2^nd Ed. Packt

岡本安晴（2014）「心理学データ分析と測定」勁草書房

岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版

Home