チョコレートの選好比率
エントロピーによる比較
2025.11：CmdStanPy対応に改訂
StanのPythonインターフェースCmdStanPy対応に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
３種類のチョコレートの嗜好を、ノーズクリップの脱着によりレトロネーザルアロマの影響をコントロールして調べたデータが星野ら（2020）に報告されている。３種のチョコレート、ガーナ、ベネズエラ、エクアドルからの選択を、ノーズクリップなしの条件とノーズクリップあり（レトロネーザルアロマの関与の除外）の２条件で行われた結果が星野ら（2020）のTable 4に掲載されているが、Table 4から各条件ごとのチョコレートの選択頻度を抜き出して表１にまとめた。各条件におけるチョコレートの選択頻度は、多項分布に従うと見なすことができる。
|
表１ チョコレートの２条件における選択（星野ら（2020）より） |
|||
|
|
ガーナ |
ベネズエラ |
エクアドル |
|
ノーズクリップなし |
55 |
24 |
16 |
|
ノーズクリップあり |
38 |
27 |
30 |
多項分布は、２項分布の拡張と考えることができる。３種のチョコレートから１つを選ぶという実験において、それぞれのチョコレートが選ばれた頻度を
、
、
とおく。それぞれのチョコレートが選ばれる確率を
、
、
とおけば、多項分布は次式で表される。
ここで、
である。
式（２）の条件を満たすパラメータ
は、事前分布としてDirichlet分布を用いると得ることができる。Dirichlet分布は３変量の場合、次式で表される。
ここで、
である。
パラメータの値を
と置けば、確率分布（３）は一様分布を表す。
表1のデータにおいて、「ノーズクリップなし」条件と「ノーズクリップあり」条件の２条件における多項分布のそれぞれの３つの確率、、の比較をエントロピーで行う。
エントロピーは、次式
で与えられ、３つの確率の値が等しいとき最大値をとり、どれか１つに確率が集まっている（すなわち、どれかの確率が１である）とき最小値をとる。
上の確率モデルのベイズ分析を行うStanスクリプトをリスト１のように用意した。Stan（PyStan2）については、岡本（2019）で説明している。
リスト１ 多項分布によるベイズ分析のStanスクリプト（choco_choice.stan）
data {
array[3] int freq_1;
array[3] int freq_2;
}
transformed data {
vector[3]
alpha;
alpha[1]
= 1; alpha[2] = 1; alpha[3] = 1;
}
parameters {
simplex[3] p1;
simplex[3] p2;
}
model {
p1 ~ dirichlet(alpha);
p2 ~ dirichlet(alpha);
freq_1 ~ multinomial(p1);
freq_2 ~ multinomial(p2);
}
リスト１のStanスクリプトを用いて表１のデータを分析するPythonスクリプトを、リスト２のように用意した。スクリプトファイルは圧縮ファイルfiles_freq.zipにまとめたので、ダウンロード・展開すれば自由に利用できる。
リスト２ 表１のデータのベイズ分析スクリプト（AnalFreq.py）
from cmdstanpy
import CmdStanModel
import numpy
as np
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
import japanize_matplotlib
model = CmdStanModel(stan_file='choco_choice.stan')
freq_noclip = [55, 24, 16]
freq_clip = [38, 27, 30]
Data = {'freq_1': freq_noclip,
'freq_2':freq_clip}
fit = model.sample(data=Data)
print(fit.summary())
inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換
az.plot_trace(inf_data)
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型に変換
p1_smpls = []
p2_smpls = []
for i in
range(3):
p1_smpls.append(fit[f'p1[{i+1}]'])
p2_smpls.append(fit[f'p2[{i+1}]'])
p1_smpls = np.array(p1_smpls).T
p2_smpls = np.array(p2_smpls).T
def calc_entropy(p):
return -p.T[0]*np.log2(p.T[0]) -p.T[1]*np.log2(p.T[1]) \
-p.T[2]*np.log2(p.T[2])
entropies1 = calc_entropy(p1_smpls)
entropies2 = calc_entropy(p2_smpls)
ent1_hat = np.median(entropies1)
ent2_hat = np.median(entropies2)
prob_e2_e1 = np.mean(entropies2
- entropies1 > 0)
sb.kdeplot(entropies1, label = 'No-clip')
sb.kdeplot(entropies2, label = 'Clip')
plt.title(('Entropy: No-clip={0:.2f},
Clip={1:.2f}\n' +
'Prob(Clip > No-clip)={2:.3f}').
format(ent1_hat, ent2_hat, prob_e2_e1))
plt.legend()
plt.savefig('FigEntropies.png')
plt.show()
p1_1_hat = np.median(p1_smpls.T[0])
p1_2_hat = np.median(p1_smpls.T[1])
p1_3_hat = np.median(p1_smpls.T[2])
sb.kdeplot(p1_smpls.T[0], label = 'p1[1]:ガーナ')
sb.kdeplot(p1_smpls.T[1], label = 'p1[2]:ベネズエラ')
sb.kdeplot(p1_smpls.T[2], label = 'p1[3]:エクアドル')
plt.legend()
plt.title('No-clip\nガーナ={0:.2f}, ベネズエラ={1:.2f}, エクアドル={2:.2f}'.format(
p1_1_hat, p1_2_hat, p1_3_hat))
plt.savefig('Figp1.png')
plt.show()
p2_1_hat = np.median(p2_smpls.T[0])
p2_2_hat = np.median(p2_smpls.T[1])
p2_3_hat = np.median(p2_smpls.T[2])
sb.kdeplot(p2_smpls.T[0], label = 'p2[1]:ガーナ')
sb.kdeplot(p2_smpls.T[1], label = 'p2[2]:ベネズエラ')
sb.kdeplot(p2_smpls.T[2], label = 'p2[3]:エクアドル')
plt.legend()
plt.title('Clip\nガーナ={0:.2f},
ベネズエラ={1:.2f}, エクアドル={2:.2f}'.format(
p2_1_hat, p2_2_hat, p2_3_hat))
plt.savefig('Figp2.png')
plt.show()
リスト２において、グラフにおいて日本語表記するためにjapanize_matplotlibをimportしていることに注意。
リスト２のPythonスクリプトファイル（AnalFreq.py）とリスト１のStanスクリプトファイル（choco_choice.stan）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされている環境で、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/files_freq$ python AnalFreq.py
上のコマンドを実行すると、Stanスクリプトのコンパイルが始まる。このコンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了後、トレース図が表示される（図１）。
図１
左に事後分布のグラフ、右にサンプリングのトレース図が表示されている。
図１のWindowを閉じると、図２のグラフが表示される。
図２
ノーズクリップなし（No-clip）とノーズクリップあり（Clip）のそれぞれの条件におけるエントロピーの事後分布のグラフである。ノーズクリップありのエントロピーの方がノーズクリップなしのエントロピーより大きい事後確率が0.992であることが示されている。ノーズクリップあり条件のエントロピーの方が大きいので、ノーズクリップあり条件での３つのチョコレートの選好確率の違いは、ノーズクリップなし条件の場合より顕著ではないと言える。
図２のWindowを閉じると、図３のグラフが表示される。
図３
「ノーズクリップなし」条件での３つのチョコレートの選好確率の事後分布のグラフが表示されている。ガーナが他の２つより好まれていることが明瞭に示されている。
図３のWindowを閉じると、図４のWindowが表示される。
図４
「ノーズクリップあり」条件での３つのチョコレートの選好確率の事後分布のグラフが示されている。図３の分布と比べると、３つのチョコレートの選好確率の違いは不明瞭になっている。このことは、図２のエントロピーの事後分布の違いに対応するものである。
図４のWindowを閉じると、スクリプトの実行終了である。
星野亜由美、石川華衣、岡田綾、飯田文子（2020）．カカオ豆の香がチョコレート嗜好に及ぼす影響．日本官能評価学会誌、24、29-42．