Up

独立な２つのサンプルの平均値と分散のベイズ分析

 

２つのデータセットとがそれぞれ独立に２つの母集団からサンプリングされたとする。それぞれの母集団の確率分布は、平均と、標準偏差との正規分布であるとする。すなわち、

である。

　このモデルのPyStan 3対応Stanスクリプトをリスト１のように用意した。なお、スクリプトファイルなどはまとめて圧縮ファイルDiffIndMeansFiles.zipとした。ダウンロード展開すれば自由に利用できる。岡本安晴「いまさら聞けないPythonでデータ分析」ではPyStan 2について説明している。PyStan 3とPyStan 2の違いについては、このウェブサイトなどで説明している。

 

 

リスト１　Stanスクリプト（diff_model.stan）

 

data {

    int n1;

    array[n1] real x1;

    int n2;

    array[n2] real x2;

    //

    //     Values for rescaling

    //

    real dmean1;

    real dsgm1;

    real dmean2;

    real dsgm2;

}

parameters {

    real vmu1;

    real vmu2;

    real<lower=0.0> vsgm1;

    real<lower=0.0> vsgm2;

}

transformed parameters {

    real mu1;

    real mu2;

    real sgm1;

    real sgm2;

    //

    //      Rescaling of parameters

    //

    mu1 = vmu1 * dsgm1 + dmean1;

    mu2 = vmu2 * dsgm2 + dmean2;

    sgm1 = vsgm1 * dsgm1 + 0.001;

    sgm2 = vsgm2 * dsgm2 + 0.001;

}

model {

    vmu1 ~ normal(0.0, 10.0);

    vmu2 ~ normal(0.0, 10.0);

    vsgm1 ~ exponential(0.1);

    vsgm2 ~ exponential(0.1);

    x1 ~ normal(mu1, sgm1);

    x2 ~ normal(mu2, sgm2);

}

generated quantities {

    real diff;

    real dsgm;

    diff = mu1 - mu2;

    dsgm = sgm1 - sgm2;

}

 

 

PyStan 3では、事前分布の分散が大き過ぎると警告メッセージが出される。これを避けるために、事前分布を

 

    vmu1 ~ normal(0.0, 10.0);

    vmu2 ~ normal(0.0, 10.0);

    vsgm1 ~ exponential(0.1);

    vsgm2 ~ exponential(0.1);

 

と設定してから、母集団のパラメータの変動範囲をデータに合わせて

 

    mu1 = vmu1 * dsgm1 + dmean1;

    mu2 = vmu2 * dsgm2 + dmean2;

    sgm1 = vsgm1 * dsgm1 + 0.001;

    sgm2 = vsgm2 * dsgm2 + 0.001;

 

とスケーリングしている。

リスト１のStanスクリプトを用いてベイズ分析を行うスクリプトを、リスト２のスクリプトのように用意した。リスト２のスクリプトをWindowsで実行する場合は、WindowsにUbuntuをインストールして、UbuntuにインストールされたAnacondaにPython仮想環境を用意して実行できる。以下の実行例は、Python 3.11仮想環境でのものである（2024.02.04）。WindowsへのUbuntuのインストールなどは、このウェブサイトなどで説明している。

ベイズ分析については、岡本安晴（2019）、岡本安晴（2014）などで説明している。

リスト２のスクリプトを実行すると、まず、入力データファイル名と出力ファイル名を聞いてくる。入力データファイルはExcel形式のファイルを用意し、出力ファイルは適当なテキストファイル名を設定する。

 

(py311) *****/DiffIndMeansFiles$ python AnalIndMeans.py

入力データファイル（*.xlsx）= data.xlsx

出力ファイル（*.txt）= temp.txt

 

入力データファイルの例として図１のものを考える。

図１　Excelデータ（data.xlsx）

 

A市とB市からランダムにそれぞれ１５人と２０人を選び、血圧を測ったものである（仮想データ）。１行目に変数名を書く。1列目はデータの識別値であるが、分析には用いられない。２列目と３列目にデータの値を書く。

出力ファイル名を設定すると、データの読み込みが始まり、ヒストグラムが表示される（図２）。

図２

 

図２のウィンドウを閉じると、Stanスクリプト（リスト１）のビルドが始まり、MCMCサンプリングが行われる。サンプリング後、トレース図が表示される（図３）。

図３

 

図３のウィンドウを閉じると、パラメータとの事後分布のグラフが表示れる（図４）。

図４

 

２つの事後分布は、明瞭に分離されている。タイトルの上段にの事後分布の四分位数、下段にの事後分布の四分位数が表示されている。分布の代表値としてGelmanら（2021）は、平均値は不安定であるとして、中央値を勧めている。

図４のウィンドウを閉じると、の事後分布のグラフが表示される（図５）。

図５

 

の事後分布の95％HDIが［-11.54, -6.04］となっていて、０より負の領域にあることが明瞭である。すなわち、母集団の平均値はよりの方が大きいことが明白であると言える（図４参照）。

図５のウィンドウを閉じると、との事後分布のグラフが表示される（図６）。

図６

 

２つの事後分布のグラフはほとんど重なっている。

図６のウィンドウを閉じると、パラメータの差の事後分布のグラフが表示される（図７）。

図７

 

95％HDIの中ほどに原点0が位置している。すなわち、２つのパラメータの差は０と考えてもよさそうである。

図７のウィンドウを閉じると、プログラムの実行終了である。端末画面には、出力ファイルtemp.txtが保存されたことが表示されている。

 

temp.txt was saved.

 

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

 

入力データファイル: data.xlsx

 

A市

[126, 122, 120, 132, 129, 124, 124, 125, 124, 124, 121, 121, 130, 125, 129]

 

B市

[126, 139, 137, 130, 139, 132, 128, 137, 134, 137, 131, 133, 136, 131, 130, 134, 141, 135, 136, 133]

 

A市

平均値 = 125.07　　標準偏差 = 3.43

 

B市

平均値 = 133.95　　標準偏差 = 3.81

 

Summary...

          mean     sd   hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

vmu1    -0.003  0.302   -0.573    0.560      0.005    0.005    3340.0    2265.0    1.0

vmu2     0.000  0.245   -0.483    0.449      0.004    0.004    3469.0    2492.0    1.0

vsgm1    1.140  0.248    0.745    1.606      0.005    0.004    3188.0    2139.0    1.0

vsgm2    1.098  0.191    0.783    1.464      0.003    0.003    3551.0    2576.0    1.0

mu1    125.058  1.036  123.098  126.989      0.018    0.013    3340.0    2265.0    1.0

mu2    133.951  0.935  132.106  135.664      0.016    0.011    3469.0    2492.0    1.0

sgm1     3.918  0.850    2.561    5.518      0.016    0.012    3188.0    2139.0    1.0

sgm2     4.187  0.728    2.986    5.584      0.013    0.010    3551.0    2576.0    1.0

diff    -8.894  1.389  -11.524   -6.283      0.025    0.018    3251.0    2507.0    1.0

dsgm    -0.270  1.127   -2.362    1.860      0.020    0.017    3224.0    2712.0    1.0

 

Posterior distributions of mu1 and mu2

A市:  Q1 = 124.41, Q2 = 125.06, Q3 = 125.71

B市:  Q1 = 133.33, Q2 = 133.94, Q3 = 134.55

 

Posterior distribution of mu1 - mu2...

mu1-m2: Q1 = -9.77, Q2 = -8.89, Q3 = -7.97

95% HDI = [-11.54, -6.04]

 

Posterior distributions of sgm1 and sgm2...

A市:  Q1 = 3.33, Q2 = 3.77, Q3 = 4.38

B市:  Q1 = 3.68, Q2 = 4.10, Q3 = 4.58

 

Posterior distribution of sgm1 - sgm2...

sgm1-sgm: Q1 = -0.99, Q2 = -0.31, Q3 = 0.43

95% HDI = [-2.35, 2.11]

 

 

プログラムの実行中に表示されるグラフは、グラフのウィンドウの保存ボタンをクリックすれば保存される。

 

 

リスト２　独立なサンプルのベイズ分析（AnalIndMeans.py）

 

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import japanize_matplotlib

import scipy.stats as ss

import stan

import arviz as az

import seaborn as sb

 

finnm = input('入力データファイル（*.xlsx）= ')

txtoutnm = input('出力ファイル（*.txt）= ')

txtout = open(txtoutnm, 'w')

#

#      Excelファイルの読み込み

#

xlsx_fl = pd.ExcelFile(finnm)                       #   ExcelファイルからExcelFile型オブジェクトを生成

data_xlsx = pd.read_excel(xlsx_fl,  header = None)  #   ExcelFile型オブジェクトからDataFrame型オブジェクトを生成

data = data_xlsx.values                             #   Dataframe型オブジェクトからnumpyの配列オブジェクトを生成

print('\n', data)

 

txtout.write('\n入力データファイル: {}\n'.format(finnm))

 

VarNameA = data[0][1]

VarNameB = data[0][2]

 

dataA = []

dataB = []

for i in range(1,len(data)):

    if not pd.isnull(data[i][1]):

        dataA.append(data[i][1])

    if not pd.isnull(data[i][2]):

        dataB.append(data[i][2])

print('dataA =\n', dataA)

print('dataB =\n', dataB)

 

txtout.write('\n{0}\n'.format(VarNameA))

txtout.write('{0}\n'.format(dataA))

txtout.write('\n{0}\n'.format(VarNameB))

txtout.write('{0}\n'.format(dataB))

 

meanA = np.mean(dataA)

meanB = np.mean(dataB)

sdA = np.std(dataA)

sdB = np.std(dataB)

print(VarNameA)

print('平均値 = {0:.2f}　　標準偏差 = {1:.2f}'.format(meanA, sdA))

print(VarNameB)

print('平均値 = {0:.2f}　　標準偏差 = {1:.2f}'.format(meanB, sdB))

txtout.write('\n' + VarNameA + '\n')

txtout.write('平均値 = {0:.2f}　　標準偏差 = {1:.2f}\n'.format(meanA, sdA))

txtout.write('\n' + VarNameB + '\n')

txtout.write('平均値 = {0:.2f}　　標準偏差 = {1:.2f}\n'.format(meanB, sdB))

 

plt.title(('{0}： mean = {1:.2f}  sd = {2:.2f}'.format(VarNameA, meanA, sdA)) + \

          '\n{0}： mean = {1:.2f}  sd = {2:.2f}'.format(VarNameB, meanB, sdB),

          fontsize = 16)

plt.hist(dataA, alpha = 0.5, label = VarNameA)

plt.hist(dataB, alpha = 0.5, label = VarNameB)

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

plt.ylabel('度数', fontsize = 14)

plt.legend()

plt.tight_layout()

plt.show()

 

Data = {'n1':len(dataA), 'x1':dataA, 'n2':len(dataB), 'x2':dataB,

        'dmean1':meanA, 'dsgm1':sdA, 'dmean2':meanB, 'dsgm2':sdB}

with open('diff_model.stan') as f_stan:

    sm = stan.build(f_stan.read(), data=Data)

fit = sm.sample()

 

i_data = az.from_pystan(posterior=fit, posterior_model=sm)

 

smry = az.summary(i_data)      #     Summary of statistics of MCMC sampling

print('Summary...\n', smry)

txtout.write(f'\nSummary...\n{smry}\n')

#

#         Trace graph of MCMC sampling

#

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

              figsize=(10,7))

plt.tight_layout()

plt.show()

#

#         Posterior distributions of parameters mu1 and mu2

#

d_frame = fit.to_frame()

mu1_para = d_frame['mu1']

mu2_para = d_frame['mu2']

q1A, q2A, q3A = np.percentile(mu1_para, [25, 50, 75])

q1B, q2B, q3B = np.percentile(mu2_para, [25, 50, 75])

strA = f'{VarNameA}:  Q1 = {q1A:.2f}, Q2 = {q2A:.2f}, Q3 = {q3A:.2f}\n' + \

       f'{VarNameB}:  Q1 = {q1B:.2f}, Q2 = {q2B:.2f}, Q3 = {q3B:.2f}'

sb.kdeplot(mu1_para, color='g', label=VarNameA)

sb.kdeplot(mu2_para, color='b', label=VarNameB)

plt.xlabel('mu1 or mu2', fontsize=16)

plt.title(strA)

plt.legend()

plt.show()

txtout.write('\nPosterior distributions of mu1 and mu2\n' + strA + '\n')

#

#           Posterior distribution of difference mu1 - mu2

#

diff_para = d_frame['diff']

print(diff_para)

diff_para = np.array(diff_para)

print(diff_para)

q1D, q2D, q3D = np.percentile(diff_para, [25, 50, 75])

LvUv = az.hdi(diff_para, hdi_prob=0.95)

print('LvUv =', LvUv)

strD = f'mu1-m2: Q1 = {q1D:.2f}, Q2 = {q2D:.2f}, Q3 = {q3D:.2f}\n' + \

      f'95% HDI = [{LvUv[0]:.2f}, {LvUv[1]:.2f}]'

txtout.write('\nPosterior distribution of mu1 - mu2...\n' + strD + '\n')

sb.kdeplot(diff_para)

plt.title(strD)

plt.xlabel('mu1-m2', fontsize = 16)

plt.tight_layout()

plt.show()

#

#           Posterior distributions of parameters sgm1 and sgm2

#

sgm1_para = d_frame['sgm1']

sgm2_para = d_frame['sgm2']

q1DA, q2DA, q3DA = np.percentile(sgm1_para, [25, 50, 75])

q1DB, q2DB, q3DB = np.percentile(sgm2_para, [25, 50, 75])

strSgm = f'{VarNameA}:  Q1 = {q1DA:.2f}, Q2 = {q2DA:.2f}, Q3 = {q3DA:.2f}\n' + \

       f'{VarNameB}:  Q1 = {q1DB:.2f}, Q2 = {q2DB:.2f}, Q3 = {q3DB:.2f}'

txtout.write('\nPosterior distributions of sgm1 and sgm2...\n' + strSgm + '\n')

sb.kdeplot(sgm1_para, label=VarNameA)

sb.kdeplot(sgm2_para, label=VarNameB)

plt.xlabel('sgm1 or sgm2', fontsize=16)

plt.title(strSgm)

plt.legend()

plt.tight_layout()

plt.show()

#

#          Posterior distribution of difference sgm1 - sgm2

#

dsgm_para = d_frame['dsgm']

dsgm_para = np.array(dsgm_para)

q1DS, q2DS, q3DS = np.percentile(dsgm_para, [25, 50, 75])

LvUvS = az.hdi(dsgm_para, hdi_prob=0.95)

print('LvUvS =', LvUvS)

strD = f'sgm1-sgm: Q1 = {q1DS:.2f}, Q2 = {q2DS:.2f}, Q3 = {q3DS:.2f}\n' + \

      f'95% HDI = [{LvUvS[0]:.2f}, {LvUvS[1]:.2f}]'

txtout.write('\nPosterior distribution of sgm1 - sgm2...\n' + strD + '\n')

sb.kdeplot(dsgm_para)

plt.title(strD)

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

plt.legend()

plt.tight_layout()

plt.show()

 

txtout.close()

print(f'\n{txtoutnm} was saved.\n')

 

 

 

Home