Up

チョコレートの選好比率

 

 

３種類のチョコレートの嗜好を、ノーズクリップの脱着によりレトロネーザルアロマの影響をコントロールして調べたデータが星野ら（2020）に報告されている。３種のチョコレート、ガーナ、ベネズエラ、エクアドルからの選択を、ノーズクリップなしの条件とノーズクリップあり（レトロネーザルアロマの関与の除外）の２条件で行われた結果が星野ら（2020）のTable 4に掲載されているが、Table 4から各条件ごとのチョコレートの選択頻度を抜き出して表１にまとめた。各条件におけるチョコレートの選択頻度は、多項分布に従うと見なすことができる。

表１　チョコレートの２条件における選択（星野ら（2020）より）

 

ガーナ

ベネズエラ

エクアドル

ノーズクリップなし

55

24

16

ノーズクリップあり

38

27

30

 

　多項分布は、２項分布の拡張と見做すことができる。３種のチョコレートから１つを選ぶという実験において、それぞれのチョコレートが選ばれた頻度を、、とおき、それぞれのチョコレートが選ばれる確率を、、とおけば、多項分布は次式で表される。

ここで、

である。

式（２）の条件を満たすパラメータは事前分布としてDirichlet分布を用いると得ることができる。Dirichlet分布は３変量の場合、次式で表される。

ここで、

である。

パラメータの値を

と置けば、確率分布（３）は一様分布を表す。

　表1のデータを、「ノーズクリップなし」条件と「ノーズクリップあり」条件の２条件における多項分布として、各条件にモデル（１）を適用するベイズ分析モデルのStanスクリプトをリスト１のように用意した。多項分布のパラメータの事前分布はパラメータ値がであるDirichlet分布を設定している。

Stanについては、岡本（2019）で説明している。

 

リスト１　多項分布による分析のStanスクリプト（choco_choice.stan）

data {

    int freq_1[3];

    int freq_2[3];

}

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）

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import pystan

import platform

 

sm = pystan.StanModel(file = 'choco_choice.stan')

 

freq_noclip = [55, 24, 16]

freq_clip = [38, 27, 30]

 

Data = {'freq_1': freq_noclip, 'freq_2':freq_clip}

 

fit = sm.sampling(data = Data,

                  n_jobs = 1 if platform.system() == 'Windows' else -1)

print(fit)

 

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

 

p1_smpls = fit['p1']

p2_smpls = fit['p2']

 

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)

 

sb.kdeplot(entropies1, label = 'No-clip')

sb.kdeplot(entropies2, label = 'Clip')

plt.title('Entropy: No-clip={0:.2f}, Clip={1:.2f}'.format(ent1_hat, ent2_hat))

plt.legend()

plt.savefig('FigEntropies.png')

plt.show()

 

p1_1_hat = np.median(fit['p1'].T[0])

p1_2_hat = np.median(fit['p1'].T[1])

p1_3_hat = np.median(fit['p1'].T[2])

sb.kdeplot(fit['p1'].T[0], label = 'p1[1]:ガーナ')

sb.kdeplot(fit['p1'].T[1], label = 'p1[2]:ベネズエラ')

sb.kdeplot(fit['p1'].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(fit['p2'].T[0])

p2_2_hat = np.median(fit['p2'].T[1])

p2_3_hat = np.median(fit['p2'].T[2])

sb.kdeplot(fit['p2'].T[0], label = 'p2[1]:ガーナ')

sb.kdeplot(fit['p2'].T[1], label = 'p2[2]:ベネズエラ')

sb.kdeplot(fit['p2'].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()

 

e1_smpls = np.sort(entropies1)

e2_smpls = np.sort(entropies2)

 

def SmallerEqProp(a, v):

    i = 0

    while v > a[i]:

        i += 1

        if i >= len(a):

            break

    prop = i / len(a)

    return prop

 

sum_prop = 0.0

for v in e1_smpls:

    sum_prop += (1 - SmallerEqProp(e2_smpls, v))

prop_Larger = sum_prop / len(e1_smpls)

print('Prop_Larger =', prop_Larger)

 

リスト２のスクリプトを実行すると、以下の結果を得る。

まず、３つの確率、、に対するエントロピー

の事後分布のグラフが表示される（図１）。

図１

 

ノーズクリップあり条件のエントロピーの方が大きいので、ノーズクリップあり条件の３つのチョコレートの選好確率の違いは、ノーズクリップなし条件の場合より顕著ではないと言える。

図１のフォームを閉じると、図２のフォームが表示される。「ノーズクリップなし」条件での３つのチョコレートの選好確率の事後分布が表示されている。

図２

 

ガーナが他の２つより好まれていることが明瞭に示されている。

図２のフォームを閉じると、図３のフォームが表示される。

図３

 

「ノーズクリップあり」条件での３つのチョコレートの選好確率の事後分布が示されている。図２の分布と比べると、３つのチョコレートの選好確率の違いは不明瞭になっている。このことは、図１のエントロピーの事後分布の違いに対応するものである。

 

 

上のスクリプトでは、エントロピーの計算をPythonスクリプト（リスト２）で行っているが、これをStanスクリプトで行うようにしたものが、リスト３のStanスクリプトである。

 

リスト３　エントロピーの計算も行うStanスクリプト（choco_choice.stan）

functions {

    real neg_p_log(real p) {

        if (p > 0.0) {

            return -p * log2(p);

        }

        else {

            return 0.0;

        }

    }

}

data {

    int freq_1[3];

    int freq_2[3];

}

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);

}

generated quantities {

    real entropy1;

    real entropy2;

 

    entropy1 = neg_p_log(p1[1]) + neg_p_log(p1[2]) + neg_p_log(p1[3]);

    entropy2 = neg_p_log(p2[1]) + neg_p_log(p2[2]) + neg_p_log(p2[3]);

}

 

リスト３のStanスクリプトを用いるPythonスクリプトはリスト４のように用意した。

これらのファイルは圧縮ファイルfiles_freq_ubuntu.zipにまとめた。ダウンロード・展開すれば自由に利用できる。

リスト３およびリスト４のスクリプトは、pystan2.19.1.1/Python 3.9.13/Anaconda 3/ubuntu 20/WSL/Windows 10で実行した。

pystan 2.17.1.0/Python 3.7.13/Anaconda 3/Windows 10で実行するとエラーになる。WSLについては、このウェブページで説明している。

 

 

リスト４　リスト３のStanスクリプトを用いるPythonスクリプト（Anal_freqs.py）

import numpy as np

import matplotlib.pyplot as plt

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

 

sm = pystan.StanModel(file = 'choco_choice.stan')

 

freq_noclip = [55, 24, 16]

freq_clip = [38, 27, 30]

 

Data = {'freq_1': freq_noclip, 'freq_2':freq_clip}

 

fit = sm.sampling(data = Data)

print(fit)

 

ent1_hat = np.median(fit['entropy1'])

ent2_hat = np.median(fit['entropy2'])

 

sb.kdeplot(fit['entropy1'], label = 'No-clip')

sb.kdeplot(fit['entropy2'], label = 'Clip')

plt.title('Entropy: No-clip={0:.2f}, Clip={1:.2f}'.format(ent1_hat, ent2_hat))

plt.legend()

plt.savefig('FigEntropies.png')

plt.show()

plt.close()

 

p1_1_hat = np.median(fit['p1'].T[0])

p1_2_hat = np.median(fit['p1'].T[1])

p1_3_hat = np.median(fit['p1'].T[2])

sb.kdeplot(fit['p1'].T[0], label = 'p1[1]:Ghana')

sb.kdeplot(fit['p1'].T[1], label = 'p1[2]:Venezuela')

sb.kdeplot(fit['p1'].T[2], label = 'p1[3]:Equador')

plt.legend()

plt.title('No-clip\nGhana={0:.2f},  Venezuela={1:.2f},  Equador={2:.2f}'.format(

            p1_1_hat, p1_2_hat, p1_3_hat))

plt.savefig('Figp1.png')

plt.show()

plt.close()

 

p2_1_hat = np.median(fit['p2'].T[0])

p2_2_hat = np.median(fit['p2'].T[1])

p2_3_hat = np.median(fit['p2'].T[2])

sb.kdeplot(fit['p2'].T[0], label = 'p2[1]:Ghana')

sb.kdeplot(fit['p2'].T[1], label = 'p2[2]:Venezuela')

sb.kdeplot(fit['p2'].T[2], label = 'p2[3]:Equador')

plt.legend()

plt.title('Clip\nGhana={0:.2f},  Venezuela={1:.2f},  Equador={2:.2f}'.format(

            p2_1_hat, p2_2_hat, p2_3_hat))

plt.savefig('Figp2.png')

plt.show()

plt.close()

 

ck_count = 0

e1_smpls = fit['entropy1']

e2_smpls = fit['entropy2']

e1_smpls = np.sort(e1_smpls)

e2_smpls = np.sort(e2_smpls)

 

def SmallerEqProp(a, v):

    i = 0

    while v > a[i]:

        i += 1

        if i >= len(a):

            break

    prop = i / len(a)

    return prop

 

sum_prop = 0.0

for v in e1_smpls:

    sum_prop += (1 - SmallerEqProp(e2_smpls, v))

prop_Larger = sum_prop / len(e1_smpls)

print('Prop_Larger =', prop_Larger)

 

 

pystan2.19.1.1/Python 3.9.13/Anaconda 3/ubuntu 20/WSL/Windows 10で実行するときは、matplotlibによるグラフ表示はないが、savefigメソッドでグラフを保存することができる。また、複数のグラフを描くときは、各グラフ描画の最後でclose()メソッドを実行する必要がある。メソッドclose()を実行せずに次のグラフ描画に進むと、複数のグラフが重ね書きされる。

 

 

引用文献

星野亜由美、石川華衣、岡田綾、飯田文子（2020）．カカオ豆の香がチョコレート嗜好に及ぼす影響．日本官能評価学会誌、24、29-42．

 

 

Up