チョコレート選好の一貫性
星野ら（2020）のTable 4によれば、「ノーズクリップなし」条件でチョコレートの選択が「ノーズクリップあり」条件でも同じである対象者は表１に示す人数であった。
|
表１ チョコレートの選好データ |
||
|
|
ノーズクリップなし |
ノーズクリップありでも 同じチョコレートを選択 |
|
ガーナ |
55 |
31 |
|
ベネズエラ |
24 |
8 |
|
エクアドル |
16 |
9 |
「ノーズクリップなし」条件で選んだチョコレートが、さらに「ノーズクリップあり」条件でも選ばれる頻度は２項分布に従うと考えることができる。すなわち、「ノーズクリップなし」条件であるチョコレートが
人に選ばれ、その
人の内
人が同じチョコレートを「ノーズクリップあり」条件で選んでいるとき、
である。ここで、
はチョコレート
に対して「ノーズクリップなし」で選ばれたときに「ノーズクリップあり」でも選ばれる確率である。
これらの選択確率は実験で用いられたチョコレートという共通点があり、何らかのまとまりがあると考えられる。いま、選択確率をロジット変換した値が共通の正規分布に従っていると考える（階層モデル；hierarchical model）。ロジット変換するのは、確率の変域が０から１の値であり、これらをロジット変換すると値の変域はー∞から＋∞となり、正規分布が使えるからである。
ここで、
である。
式（１）の右辺の正規分布のパラメータ
と
はhyperparameterと呼ばれている。これらのhyperparameterの事前分布をできるだけ平らな事前分布として次式で与える。
これらhyperparameterの事前分布は、hyperprior分布と呼ばれている。出来るだけ平らであると同時に事後分布の安定性を考えて設定した。
上の設定の確率モデルのStanスクリプトをリスト１に示す。
リスト１ 階層モデルのStanスクリプト（persistence_hi.stan）
data {
int
Ng;
int
kg;
int
Nv;
int
kv;
int
Ne;
int
ke;
}
parameters {
real hm;
real<lower = 0.0> hsgm;
real hg;
real hv;
real he;
}
transformed parameters {
real<lower = 0, upper =
1> pg;
real<lower = 0, upper =
1> pv;
real<lower = 0, upper =
1> pe;
pg
= inv_logit(hg);
pv
= inv_logit(hv);
pe
= inv_logit(he);
}
model {
hm
~ normal(0, 10);
hsgm
~ uniform(0, 20);
hg ~ normal(hm, hsgm);
hv
~ normal(hm, hsgm);
he ~ normal(hm, hsgm);
kg ~ binomial(Ng, pg);
kv
~ binomial(Nv, pv);
ke
~ binomial(Ne, pe);
}
リスト１のStanスクリプトを用いて表１のデータを分析するPythonスクリプトをリスト２に示す。
Stanについては、岡本（2019）で説明している。
これらのスクリプトファイルは圧縮ファイルfileschocoPersist.zipにまとめた。ダウンロード・展開すれば、自由に利用できる。
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（AnaPersist.py）
import numpy
as np
import matplotlib.pyplot
as plt
import seaborn
as sb
import pystan
import platform
sm = pystan.StanModel(file
= 'persistence_hi.stan')
print('Compilig
ended.')
Data = {'Ng':55, 'kg':31, 'Nv':24,
'kv':8, 'Ne':16, 'ke':9}
fit = sm.sampling(data
= Data,
n_jobs = 1 if platform.system()
== 'Windows' else -1)
print(fit)
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.title('pv_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['hg'], label = 'hg:Ghana')
sb.kdeplot(fit['hv'],
label = 'hv:Venezuala')
sb.kdeplot(fit['he'], label = 'he:Equador')
plt.legend()
plt.savefig('FigH.png')
plt.show()
sb.kdeplot(fit['hm'],
label = 'Hyper Mu')
plt.legend()
plt.savefig('FigHM.png')
plt.show()
sb.kdeplot(fit['hsgm'],
label = 'Hyper Sigma')
plt.legend()
plt.savefig('FigHS.png')
plt.show()
plt.hist(fit['hsgm'],
label = 'Hyper Sigma')
plt.legend()
plt.savefig('FigHSgmHist.png')
plt.show()
リスト２のスクリプトを実行すると以下の結果を得る。
図１のグラフは、３つのチョコレート、ガーナ、ベネズエラ、エクアドル、における２項分布の確率pg、pv、peの事後分布である。ベネズエラの確率分布が他の２つより小さい、すなわち「ノーズクリップなし」条件から「ノーズクリップあり」条件に変わると、他のチョコレートを選びやすいことが判る。
図１
図２は、pg、pv、peのロジット変換（式（２））の値hg、hv、heの事後分布である。
図２
図３は、ハイパーパラメータの事後分布である。
図３
図４は、ハイパーパラメータの事後分布である。
図４
図４のsigmaの分布が、負の領域まで延びているのは、KDEを用いたグラフであることによる。
ヒストグラムは図５のようになる。
図５
分布は０以上の領域になっている。
星野亜由美、石川華衣、岡田綾、飯田文子（2020）．カカオ豆の香がチョコレート嗜好に及ぼす影響．日本官能評価学会誌、24、29-42．