チョコレート選好の一貫性のベイズ分析
2025.11：CmdStanPy対応に改訂
StanのPythonインターフェースCmdStanPy対応に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
ノーズクリップしたときにチョコレートの選択がどのように影響を受けるかのデータに対して、２項分布を用いたベイズ分析をハイパーパラメータを設定して行った。
星野ら（2020）のTable 4によれば、「ノーズクリップなし」条件でチョコレートの選択が「ノーズクリップあり」条件でも同じである対象者は表１に示す人数であった。
|
表１ チョコレートの選好データ |
||
|
|
ノーズクリップなし |
ノーズクリップありでも 同じチョコレートを選択 |
|
ガーナ |
55 |
31 |
|
ベネズエラ |
24 |
8 |
|
エクアドル |
16 |
9 |
「ノーズクリップなし」条件で選んだチョコレートが、さらに「ノーズクリップあり」条件でも選ばれる頻度は２項分布に従うと考えることができる。すなわち、「ノーズクリップなし」条件であるチョコレートが
人に選ばれ、その
人の内
人が同じチョコレートを「ノーズクリップあり」条件で選んでいるとき、
である。ここで、
はチョコレート
に対して「ノーズクリップなし」で選ばれたときに「ノーズクリップあり」でも選ばれる確率である。
これらの選択確率は実験で用いられたチョコレートという共通点があり、何らかのまとまりがあると考えられる。いま、選択確率が共通の母集団からサンプリングされた値であると考え、共通の母集団をベータ分布と置く（Gelmanら、2014）。すなわち、
とする。パラメータ
と
は、ハイパーパラメータである。
「ノーズクリップなし」条件で人に選ばれたチョコレートが、さらに「ノーズクリップあり」条件で人に選ばれたとき、
が成り立つとする。
上記のモデルによるベイズ分析を行うStanスクリプトをリスト１のように用意した。
リスト１ 階層モデルのStanスクリプト（persistence_hi.stan）
data {
int Ng;
int kg;
int Nv;
int kv;
int Ne;
int ke;
}
parameters {
real<lower=0.1>
a;
// hyperparameter
real<lower=0.1>
b;
// hyperparameter
real<lower = 0, upper =
1> pg;
real<lower = 0, upper =
1> pv;
real<lower = 0, upper =
1> pe;
}
model {
a ~ exponential(0.1);
b ~ exponential(0.1);
pg
~ beta(a, b);
pv
~ beta(a, b);
pe ~ beta(a, b);
kg ~ binomial(Ng, pg);
kv
~ binomial(Nv, pv);
ke
~ binomial(Ne, pe);
}
リスト１のStanスクリプトを用いて表１のデータを分析するPythonスクリプトをリスト２に示す。
Stan（PyStan2）については、岡本（2019）で説明している。
これらのスクリプトファイルは圧縮ファイルfileschocoPersist.zipにまとめた。ダウンロード・展開すれば、自由に利用できる。
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（AnaPersist.py）
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
import platform
from cmdstanpy
import CmdStanModel
model = CmdStanModel(stan_file='persistence_hi.stan')
Data = {'Ng':55, 'kg':31, 'Nv':24,
'kv':8, 'Ne':16, 'ke':9}
fit = model.sample(data=Data,
adapt_delta=0.99)
print(fit.diagnose())
print(fit.summary())
inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換
az.plot_trace(inf_data, figsize=(12,8))
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型に変換
pg_hat = np.median(fit['pg'])
pv_hat = np.median(fit['pv'])
pe_hat = np.median(fit['pe'])
sb.kdeplot(fit['pg'],
label = 'Ghana')
sb.kdeplot(fit['pv'],
label = 'Venezuela')
sb.kdeplot(fit['pe'], label = 'Equador')
plt.xlabel('')
plt.title('pg_hat={0:.3f}, pv_hat={1:.3f}, pe_hat={2:.3f}'.format(pg_hat, pv_hat, pe_hat))
plt.legend()
plt.savefig('Fig_persist_hi.png')
plt.show()
sb.kdeplot(fit['a'], label='a')
sb.kdeplot(fit['b'], label='b')
plt.title('Posterior distributions of a and b')
plt.legend()
plt.show()
a_hat = np.median(fit['a'])
b_hat = np.median(fit['b'])
rng = np.random.default_rng()
n = len(fit['a'])
xcoord = np.linspace(0,
1, 1000)
idx_smpls = rng.choice(n,
size=100, replace=False)
for i in idx_smpls:
ycoord
= ss.beta.pdf(xcoord, fit['a'][i],
fit['b'][i])
plt.plot(xcoord, ycoord, c='b', alpha=0.3)
ycoord = ss.beta.pdf(xcoord,
a_hat, b_hat)
plt.plot(xcoord, ycoord, c='r', label=f'a_hat={a_hat:.3f}\nb_hat={b_hat:.3f}')
plt.legend()
plt.title('Beta distributions', fontsize=16)
plt.tight_layout()
plt.show()
リスト１のStanスクリプトファイル（persistence_hi.stan）とリスト２のPythonスクリプトファイル（AnalPersit.py）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、次のコマンドをCmdStanPyのインストールされている環境で実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/fileschocoPersist$ python AnalPersist.py
上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了後、トレース図が表示される（図１）。
図１
左側に事後分布のグラフ、右側にサンプリングのトレース図が表示されている。
図１のWindowを閉じると、図２のグラフが表示される。
図２
３つのチョコレート、ガーナ、ベネズエラ、エクアドル、における２項分布の確率pg、pv、peの事後分布のグラフである。ベネズエラの確率分布が他の２つより小さい、すなわち「ノーズクリップなし」条件から「ノーズクリップあり」条件に変わると、他のチョコレートを選びやすいことが判る。
図２のWindowを閉じると、図３のグラフが表示される。
図３
ベータ分布のパラメータaとbの事後分布のグラフである。
図3のWindowを閉じると、図４のグラフが表示される。
図４
赤の太い実線は、パラメータaとbの事後分布の中央値a_hatとb_hatの値をパタメータ値とするベータ分布のグラフである。Gelmanら（2021）は、事後分布の点推定値として、平均値より安定した中央値を勧めている。
青の細い実線は、MCMCサンプルからランダムに100対を選び、各対(a, b)に対してそれをパラメータとするベータ分布のグラフである。
図４のWindowを閉じると、スクリプトの実行終了である。
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
星野亜由美、石川華衣、岡田綾、飯田文子（2020）．カカオ豆の香がチョコレート嗜好に及ぼす影響．日本官能評価学会誌、24、29-42．