順序カテゴリ変数の平均値の比較
ベイズ分析
質問紙や調査紙では、データはカテゴリ値で得られることが多い。カテゴリ値の確率モデルとして多項分布を設定したときの２つの独立なデータの平均値のベイズ分析を試みた。カテゴリ値データは、「好きである」、「どちらでもない」、「嫌いである」というように非数値であることが普通である。これらのカテゴリ値に例えば数値１、０、―１が当てはめられているときは、それらの平均値を考えることができる。
いま、一般的に以下のように考える（Gelman et al., 2014）。
カテゴリ数はkであるとする。カテゴリ値を
各カテゴリの確率を
各カテゴリに割り当てられた数値を
とおく。このとき、平均値は
で与えられる。
各カテゴリの観測頻度を
とおくと、確率は多項分布
で与えられる。
なお、カテゴリの度数を数える簡単なスクリプト例を用意している。
データ１として頻度
データ２として頻度
が与えられたときのベイズ分析のスクリプトをリスト１とリスト２のように作成してみた。ファイルはまとめて圧縮ファイルba_ordcat_files.zipとした。
リスト１のスクリプトは、Windows上では直接実行できない。Windowsで実行するときは、WindowsにWSL（Windows Subsystem for Linux）をインストールして、WSL上にインストールされたUbuntu上でリスト１のスクリプトを実行する。Ubuntu上での実行は、miniconda上に用意したPython3.12仮想環境で行った。すなわち、リスト１のスクリプトは、Python3.12/miniconda/Ubuntu/WSL/Windows 10で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。
リスト１のスクリプトでMCMCサンプリングを行った後、リスト２のスクリプトでMCMCサンプルの分析を行った。リスト２のスクリプトは直接Windows上で実行した。WSL/Windowsでは描画ソフトが使えないようである（2025.02.14）。
スクリプトの実行例
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312) *****/ba_ordcat_files$ python
ba_ordcat.py
Input data file(*.csv) =
入力データファイルは、CSVファイルの形式で用意する（図１）。
図１ 入力データファイル例sample1.csv
１行目は、２列目からカテゴリ値を設定する。リスト１のスクリプトは、入力ファイルの読み込みにおいて日本語対応にされていない（encodingパラメータが設定されていない）ので注意。
２行目に、カテゴリ値に対応する数値を設定する。
３行目にデータ１（Data1），４行目にデータ２（Data2）の頻度データを設定する。
各行とも、１列目に変数名などを設定して、データ値は２列目以降に設定されていることに注意。
入力データファイルとリスト１のファイルおよびリスト２のファイルは、同じフォルダ内に置いておく。
入力データファイル名を設定してEnterキーを押すとデータが読み込まれる。
(py312) *****/ba_ordcat_files$ python
ba_ordcat.py
Input data file(*.csv) = sample1.csv
上の例では、図１のファイル名sample1.csvが設定されている。
データが読み込まれると、Stanスクリプトのビルディング、MCMCサンプリングが行われ、サンプリングの結果がファイルに出力されて、リスト１の実行は終了する。
Data.pkl was saved.
fit.pkl was saved.
(py312) *****/ba_ordcat_files$
リスト１の実行後、リスト２を実行する。リスト１はminiconda/Ubuntu24/WSL/Windows10のPython 3.12仮想環境で実行したが、Ubuntu/WSLではGUIが使えない（2025.02.14）ので、リスト２はWindows上で直接実行した。
リスト２のスクリプトを実行すると、リスト１のスクリプトの実行で生成されたMCMCサンプルが読み込まれ、まず平均値パラメータmean1とmean2の事後分布のグラフが表示される（図２）。
図２
データ１の平均値の分布がデータ２の平均値の分布の右側にあり、データ１の平均値の方が大きいことが示されている。事後分布の代表値として中央値（Median）が示されている（Gelman et al., 2021）。
図２のWindowを閉じると、図３のグラフが表示される。
図３
平均値パラメータの差mean1―mean2の事後分布である。平均値の差が正、すなわちmean1>mean2である事後確率は1.000であるので、ほぼ確実にmean1の方がmean2より大きいことが示されているといえる。差mean1―mean2の中央値は0.621である。
図３のWindowを閉じると、図４のグラフが表示される。
図４
データ１(Data1)およびデータ２（Data2）における各カテゴリiの選択確率
および
の事後分布の2.5％、50％、97.5％の各点を求め、折れ線で繋いだものである。したがって、帯の幅は、95％CIを表している。カテゴリC5とカテゴリC7において、95％CIに重なりが見られない。C5においてはデータ２の方が確率が高く、C7においてはデータ１の方が確率が高いと言えそうである。
なお、ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。
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.
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版
import pandas as pd
import stan
import pickle
Stan_Code = """
data {
int
k;
// カテゴリ数
array[k]
real cat_v;
// カテゴリ値
array[k] int f1;
// データ１：各カテゴリの頻度
array[k] int f2;
// データ２：各カテゴリの頻度
}
transformed data {
vector[k] a;
for
(i in 1:k)
a[i] = 1.0;
// Dirichlet分布の一様分布パラメータ
}
parameters {
simplex[k] p1;
simplex[k] p2;
}
model {
p1 ~
dirichlet(a);
// p1 depends on data a
p2 ~
dirichlet(a);
// p2 depends on data a
f1 ~
multinomial(p1);
f2 ~ multinomial(p2);
}
generated quantities {
real
mean1;
// データ１の母集団の平均値
real
mean2;
// データ２の母集団の平均値
real
diff_means;
// 母集団の平均値の差
mean1 = 0.0;
mean2 = 0.0;
for
(i in 1:k) {
mean1 += cat_v[i] * p1[i];
mean2 += cat_v[i] * p2[i];
}
diff_means = mean1 - mean2;
}
"""
csvin_nm = input('Input data
file(*.csv) = ')
pd_dframe = pd.read_csv(csvin_nm)
print(pd_dframe)
cat_names = pd_dframe.keys().values[1:]
print(cat_names)
data_names = pd_dframe.values.T[0][1:]
print(data_names)
cat_vals = pd_dframe.values[0][1:]
print(cat_vals)
data1 = pd_dframe.values[1][1:]
print('data1 =', data1)
data2 = pd_dframe.values[2][1:]
print('data2 =', data2)
"""
StanスクリプトによるMCMCサンプリング
"""
sm = stan.build(Stan_Code, data =
{'k':len(cat_names),
'cat_v':cat_vals,
'f1':data1,
'f2':data2})
fit = sm.sample(num_samples=10000)
with open('Data.pkl', 'wb') as f:
pickle.dump(pd_dframe,
f)
# データの保存
print('Data.pkl was saved.')
fit_pd_frame = fit.to_frame()
with open('fit.pkl', 'wb') as f:
pickle.dump(fit_pd_frame,
f) # MCMCサンプルの保存
print('fit.pkl was saved.')
リスト２ MCMCサンプルの分析（draw_ordcat.py）
import pickle
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
with open('Data.pkl', 'rb') as f:
pd_dframe =
pickle.load(f)
# データの読み込み
with open('fit.pkl', 'rb') as f:
fit_pd_frame =
pickle.load(f)
# MCMCサンプルの読み込み
print('pd_dframe =\n', pd_dframe)
print('fit_pd_frame =\n',
fit_pd_frame)
print(fit_pd_frame.keys())
cat_names =
pd_dframe.keys().values[1:] # カテゴリ名
print(cat_names)
data_names = pd_dframe.values.T[0][1:]
print(data_names)
cat_vals =
pd_dframe.values[0][1:]
# カテゴリ値
print(cat_vals)
means1 = fit_pd_frame['mean1']
# データ１の平均値パラメータ値サンプル
means2 = fit_pd_frame['mean2']
# データ２の平均値パラメータ値サンプル
diffs =
fit_pd_frame['diff_means']
# 平均値パラメータ値の差のサンプル
"""
パラメータ値の事後分布
"""
sb.kdeplot(means1, c='g',
label=f'mean1[{data_names[0]}]' +
f', Med.={np.median(means1):.2f}')
sb.kdeplot(means2, c='orange',
label=f'mean2[{data_names[1]}]' +
f', Med.={np.median(means2):.2f}')
plt.xlabel('')
plt.title('Posterior distributions',
fontsize=16)
plt.legend(fontsize=14)
plt.tight_layout()
plt.show()
"""
パラメータ値の差の事後分布
"""
ck_diff = (diffs > 0)
p_pos = np.mean(ck_diff)
diff_med = np.median(diffs)
plt.title(f'P(mean1>mean2)={p_pos:.3f},
Med.of mean1 - mean2={diff_med:.3f}')
sb.kdeplot(diffs)
plt.xlabel(f'$mean1[{data_names[0]}]-mean2[{data_names[1]}]$',
fontsize=14)
plt.tight_layout()
plt.show()
"""
多項分布のカテゴリ確率
"""
ncat = len(cat_names)
cat_vals = np.array(cat_vals,
dtype=float)
"""
カテゴリ確率の2.5%,
50%, 97.5%ポイント
"""
PL1 = np.zeros(ncat)
PM1 = np.zeros(ncat)
PU1 = np.zeros(ncat)
PL2 = np.zeros(ncat)
PM2 = np.zeros(ncat)
PU2 = np.zeros(ncat)
for i in range(ncat):
ps1 =
fit_pd_frame[f'p1.{i+1}']
PL1[i], PM1[i], PU1[i] =
np.percentile(ps1, [2.5, 50, 97.5])
ps2 =
fit_pd_frame[f'p2.{i+1}']
PL2[i], PM2[i], PU2[i] =
np.percentile(ps2, [2.5, 50, 97.5])
plt.plot(cat_vals, PL1, c='g',
ls='--')
plt.plot(cat_vals, PM1, c='g', ls='-',
lw=3, label=data_names[0])
plt.plot(cat_vals, PU1, c='g',
ls='--')
plt.fill_between(cat_vals, PL1, PU1,
color='g', alpha=0.5)
plt.plot(cat_vals, PL2, c='orange',
ls='--')
plt.plot(cat_vals, PM2, c='orange', ls='-',
lw=3, label=data_names[1])
plt.plot(cat_vals, PU2, c='orange',
ls='--')
plt.fill_between(cat_vals, PL2, PU2,
color='orange', alpha=0.5)
plt.xticks(cat_vals, cat_names,
fontsize=12)
plt.xlabel('Categories', fontsize=14)
plt.ylabel('$p_i$', fontsize=14)
plt.plot([], c='k', ls='-', lw=3,
label='Median')
plt.plot([], c='k', ls='--',
label='2.5% or 97.5% points')
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()
print('\n\nPosterior distributions of
pi')
print()
print(data_names[0])
print(f'{'Cat.':>10s}{'2.5%':>15s}{'50%':>15s}{'97.5%':>15s}')
for i in range(ncat):
print(f'{cat_vals[i]:>10}{PL1[i]:>15.3f}{PM1[i]:>15.3f}{PU1[i]:>15.3f}')
print()
print(data_names[1])
print(f'{'Cat.':>10s}{'2.5%':>15s}{'50%':>15s}{'97.5%':>15s}')
for i in range(ncat):
print(f'{cat_vals[i]:>10}{PL2[i]:>15.3f}{PM2[i]:>15.3f}{PU2[i]:>15.3f}')