甘味度と温度の関係
糖類間の比較のダミー変数による分析
2025.11：CmdStanPy対応に改訂
ベイズ分析のソフトStanのPythonインターフェースCmdStanPy対応に改訂した。CmdStaPyの簡単な説明を、このウェブサイトに用意した。
２種類の糖、フラクトースとグルコースの甘味度と温度の関係をそれぞれ回帰モデルで表し、２つの回帰モデルを、ダミー変数を用いて１つの回帰モデルにまとめた。
糖の甘味度と温度の関係が前橋（2011）に報告されている。前橋（2011）の第２図からフラクトースとグルコースの温度と甘味度の関係を読み取ったものを表１に示す。
|
表１ 前橋（2011）の第2図からの読み取り |
|||
|
フラクトース |
グルコース |
||
|
温度℃ |
甘味度 |
温度℃ |
甘味度 |
|
5 |
1.43 |
3 |
0.59 |
|
18 |
1.27 |
19 |
0.57 |
|
40 |
1 |
40 |
0.56 |
|
60 |
0.79 |
60 |
0.56 |
甘味度を従属変数、温度を独立変数とする回帰モデルをフラクトースおよびグルコースに対して以下のように設定することができる。
いま、ダミー変数
をフラクトースに対して１、グルコースに対して０をとるものを用意する。
である。
このとき、回帰モデルを次のようにおく。
この回帰式は、フラクトースのときは
となり、グルコースのときは
となる。
このダミー変数を用いた回帰モデルのStanスクリプトをリスト１のように用意した。
スクリプトファイルおよびデータファイルは圧縮ファイルfilesSweetTemp.zipにまとめた。ダウンロード・展開すれば自由に利用できる。
Stan（旧版）の説明は、岡本（2019）で行っている。
リスト１ ダミー変数を用いた回帰モデルの分析Stanスクリプト（ModelDummy.stan）
data {
int N;
array[N] real Y;
array[N] real X;
array[N] real G;
}
parameters {
real a1;
real a2;
real b1;
real b2;
real<lower = 0.0> sgm;
}
transformed parameters {
array[N] real a;
array[N] real b;
for (i
in 1:N) {
a[i] = a1 + a2 * G[i];
b[i] = b1 + b2 * G[i];
}
}
model {
a1 ~ normal(0, 10);
a2 ~ normal(0, 10);
b1 ~ normal(0, 10);
b2 ~ normal(0, 10);
sgm
~ exponential(0.1);
for (i
in 1:N) {
Y[i] ~ normal(a[i] + b[i] * X[i], sgm);
}
}
表１のデータを分析するために、2つのExcelファイルを図１のように用意した。
図１
図１のファイルを読み込んで、リスト１のStanスクリプトによりベイズ分析を行うPythonスクリプトをリスト２のように用意した。温度Xを中心化した値Zを独立変数としていることに注意。
リスト２ ダミー変数を用いた回帰モデルの分析Pythonスクリプト（SweetTempDummy.py）
from cmdstanpy
import CmdStanModel
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
xlfile = pd.ExcelFile('fructose.xlsx')
print(xlfile)
df = pd.read_excel(xlfile)
print(df)
col_names = df.keys()
print(col_names)
Xf = df['Temperature'].values
Yf = df['Sweetness'].values
Gf = np.ones(len(Xf))
# ダミー変数 G=1
print('Xf
=\n', Xf)
print('Yf
=\n', Yf)
print('Gf =\n', Gf)
xlfile = pd.ExcelFile('glucose.xlsx')
print(xlfile)
df = pd.read_excel(xlfile)
print(df)
col_names = df.keys()
print(col_names)
Xg = df['Temperature'].values
Yg = df['Sweetness'].values
Gg = np.zeros(len(Xf))
# ダミー変数 G = 0
print('Xg
=\n', Xg)
print('Yg
=\n', Yg)
print('Gg =\n', Gg)
X = np.concatenate((Xf,Xg))
Y = np.concatenate((Yf, Yg))
G = np.concatenate((Gf,
Gg))
print('X =\n', X)
print('Y =\n', Y)
print('G =\n', G)
model = CmdStanModel(stan_file='ModelDummy.stan')
Z = X - np.mean(X)
# 中心化
Data = {'N':len(Y),
'Y':Y, 'X':Z, 'G':G}
fit = model.sample(data=Data,
iter_warmup=5000, adapt_delta=0.95)
print(fit.summary())
inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換
az.plot_trace(inf_data, var_names=['a1','a2','b1','b2','sgm'],
figsize=(12,8))
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型に変換
sgm_hat = np.median(fit['sgm'])
sb.kdeplot(fit['sgm'])
plt.title('Sigma: median = {0:.2}'.format(sgm_hat))
plt.savefig('FigSgm.png')
plt.show()
Y_hat = np.empty(len(Y))
a1_hat = np.median(fit['a1'])
a2_hat = np.median(fit['a2'])
b1_hat = np.median(fit['b1'])
b2_hat = np.median(fit['b2'])
sb.kdeplot(fit['a1'], label = 'a1')
sb.kdeplot(fit['a2'], label = 'a2')
plt.title('a1: median={0:.2f}, a2: median={1:.2f}'.format(
a1_hat, a2_hat))
plt.legend()
plt.savefig('Figa1a2.png')
plt.show()
sb.kdeplot(fit['b1'], label = 'b1')
sb.kdeplot(fit['b2'], label = 'b2')
plt.title('b1: median={0:.4f}, b2: median={1:.4f}'.format(
b1_hat, b2_hat) + '\n' + f'Prob(b1<0)={np.mean(fit["b1"]<0)}')
plt.xticks(rotation = 45)
plt.legend()
plt.tight_layout()
plt.savefig('Figb1b2.png')
plt.show()
plt.plot(Z, Y, 'o')
print(np.min(Z),
np.max(Z))
plt.plot([np.min(Z),
np.max(Z)], [a1_hat + b1_hat*np.min(Z),
a1_hat + b1_hat*np.max(Z)],
label = 'glucose')
plt.plot([np.min(Z),
np.max(Z)],
[a1_hat + a2_hat + (b1_hat + b2_hat)*np.min(Z),
a1_hat + a2_hat + (b1_hat + b2_hat)*np.max(Z)],
label = 'fructose')
x_coord = np.arange(0,
61, 10)
z_coord = x_coord -
np.mean(X)
plt.xticks(z_coord, x_coord)
plt.xlabel('$^{\circ}C$', fontsize
= 14)
plt.ylabel('Sweetness')
plt.plot([0, 0], [np.min(Y),
np.max(Y)], label = 'Mean of X = {0:.1f}'.format(np.mean(X)))
plt.legend()
plt.tight_layout()
plt.savefig('FigRegressDummy.png')
plt.show()
print('Mean of X = {0:.5f}'.format(np.mean(X)))
ck_b1 = [1 if v < 0 else 0 for v in
fit['b1']]
print('Prop(b1 < 0) =
{0:.5f}'.format(np.sum(ck_b1) / len(ck_b1)))
リスト１のStanスクリプトファイル（ModelDummy.stan）、リスト２のPythonスクリプトファイル（SweetTempDummy.py）、および入力データファイル（fructose.xlsxとglucose.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、次のコマンドを、CmdStanPyがインストールされている環境で実行する。CmdStanPyの簡単な説明は、このウェブサイトに用意した。
(stan) *****/filesSweetTemp$ python SweetTempDummy.py
上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。このコンパイルには多少時間が掛かる。コンパイル後、MCMCサンプルのトレース図が表示される（図２）。
図２
左側に事後分布のグラフ、右側にサンプリングのトレース図が表示されている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
sigmaの事後分布のグラフである。事後分布の代表値として中央値が表示されている。Gelmanら（2021）は、代表値として、平均値より安定した中央値を勧めている。
図４
a1とa2の事後分布のグラフである。a1はグルコースの回帰直線のy切片、a1+a2はﾌラクトースの回帰直線のｙ切片である。回帰直線は、中心化された独立変数Ｚに対して求められている。したがって、データの独立変数の平均値におけるグルコースの回帰直線の値が正の値であること、ﾌラクトースの回帰直線の独立変数の平均値における値がグルコースより高いことが示されている。
図５
パラメータb1とb2の事後分布のグラフである。パラメータb1の中央値はー0.0005であり、事後分布が負である確率は0.91であるので、グルコースの甘味度は、温度とともに減少し、減少の度合いは１℃あたり0.0005であると言える。
ﾌラクトースの温度に対する甘味減少度はb1+b2であるが、b2の事後分布の中央値はー0.0112であり、事後分布（赤色のグラフ）の範囲は明瞭に負の領域にある。ﾌラクトースの甘味減少度は、グルコースより速いと言える。
図５のWindowを閉じると、図６のグラフが表示される。
図６
ﾌラクトースとグルコースの甘味度の温度に対する回帰直線である。中央の赤の垂直線は、温度データの平均値の位置を表している。この測定温度の平均値においてﾌラクトースの方が甘味度が高いことは、図４のグラフに対応している。温度による甘味度の減少がﾌラクトース方が急であることは、図５のグラフに対応している。
図６のWindowを閉じると、スクリプトの実行終了である。
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
前橋健二（2011）．甘味の基礎知識．日本醸造協会誌、106、818-825．