匂いの順応曲線
ベイズ回帰分析例
回帰分析（regression analysis）では、回帰モデル（regression model）と呼ばれる数値データである従属変数（dependent
variable, outcome）と独立変数（independent variable,
predictor/explanatory variable, feature）の関係が求められる（Guttag, 2016, p.373; Raschka &
Mirjalili, 2019, p.4）。回帰分析の例として、匂いの順応曲線を求めてみる。
高島（1997）のFig. 3に２つの匂いの順応データが示されている。高島のFig. 3から１分ごとの値を目で読み取ったものを図1に示す。最初の匂いの強さを100としたときの、イソ吉草酸とβ―フェニールエチルアルコールの匂いの強さを１分ごとに表したものである。
図１
回帰モデルを、指数関数を用いて次のように設定する。
ｔ分後における匂いの強さ
を
とおく。
パラメータの事前分布として一様分布を設定したときのStanスクリプトをリスト.1のように用意した。
StanのPythonにおける使用例は、岡本（2019）で説明している。
リスト１のStanスクリプトを利用して図１のデータをベイズ分析するPythonスクリプトをリスト２のように用意した。
リスト１のファイル、リスト２のファイルと図１のデータファイルは、圧縮ファイルfiles.zipとしてまとめた。ダウンロード・展開すれば自由に利用できる。
リスト２のスクリプトを実行すると、MCMCサンプリングの後、まず図２の回帰曲線が表示される。
図２
図２のフォームを閉じると、回帰曲線のパラメータaの事後分布のグラフ（図３）が表示される。
図３
図２および図３のグラフは、ファイルとして自動的に保存されている。
時刻
イソ吉草酸 Βーフェニールエチルアルコール
tにおける匂いの強さに対する
時間後の強さの比は次式で表される。
すなわち、時間後の強さの減少比率は時刻ｔに関わらず一定である。
イソ吉草酸の場合、
であるので、
となる。１分経過ごとに86％に減少する。
β―フェニールエチルアルコールの場合、
であるので、
となる。１分経過ごとに97％に減少する。
Guttag, J. V. (2016). Introduction
to computation and programming using Python with application to understanding
data, 2nd ed. The MIT Press.
岡本安晴（2019）．いまさら聞けないPythonでデータ分析．丸善出版．
Raschka, S., & Mirjalili, V. (2019). Python machine learning, 3rd ed. Packt.
高島靖弘（1997）．香りと官能評価．日本官能評価学会誌、１，10-17．
リスト１ 順応曲線の回帰分析のStanスクリプト（ファイル名：reg_exp.stan）
data {
int
N;
real x[N];
real y[N];
}
parameters {
real<upper = 0.0> a;
real<lower = 0.0> sa;
real<lower = 0.0> sb;
}
transformed parameters {
real<lower = 0.0> sgm[N];
for (i
in 1:N) {
sgm[i] = sa
+ sb * i;
}
}
model {
a ~ uniform(-2, 0);
sa
~ uniform(0, 10);
sb
~ uniform(0, 10);
for (i
in 1:N) {
y[i] ~ normal(100 * exp(a * x[i]), sgm[i]);
}
}
リスト２ リスト１のStanスクリプトを利用したPythonスクリプト例（ファイル名：AdapCurve.py）
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import pystan
import pickle
import seaborn
as sb
def EstMap(samples,
a = 0.025, n_points = 10000):
import numpy as
np
import scipy.stats
as ss
Lp,
Up = np.percentile(samples, [100 * a/2, 100 * (1 -
a/2)])
coord
= np.linspace(Lp, Up, n_points)
est_pdf
= ss.gaussian_kde(samples).pdf(coord)
map_idx
= np.argmax(est_pdf)
MAP_Est
= coord[map_idx]
return MAP_Est
xlfile = pd.ExcelFile('Time_course.xlsx')
df = pd.read_excel(xlfile)
col_names = df.keys()
print('col_names:\n',
col_names)
col_names = [v for v in col_names]
print('col_names:\n',
col_names)
df_data = df.values
print('df_data\n',
df_data)
data = np.asfarray(df_data)
print(data)
data_time = data.T[0]
data_B_ph = data.T[1]
data_Iso = data.T[2]
for t, b, i
in zip(data_time, data_B_ph,
data_Iso):
print(t, b, i)
sm = pystan.StanModel(file
= 'reg_exp.stan')
with open('sm.pickle',
'wb') as f:
pickle.dump(sm, f)
stan_data = {'N':len(data_time), 'x': data_time, 'y':data_Iso}
fit_Iso = sm.sampling(data
= stan_data, n_jobs = 1)
print(fit_Iso)
a_smpls_Iso = fit_Iso['a']
a_hat_Iso = EstMap(a_smpls_Iso)
pred_Iso = 100*np.exp(a_hat_Iso * data_time)
plt.plot(data_time, pred_Iso, c = 'orange')
plt.plot(data_time, data_Iso, '^', c = 'orange', label = 'Isovaleric')
with open('sm.pickle',
'rb') as f:
sm
= pickle.load(f)
stan_data = {'N':len(data_time), 'x': data_time, 'y':data_B_ph}
fit_B_ph = sm.sampling(data
= stan_data, n_jobs = 1)
print(fit_B_ph)
a_smpls_B_ph = fit_B_ph['a']
a_hat_B_ph = EstMap(a_smpls_B_ph)
pred_B_ph = 100*np.exp(a_hat_B_ph * data_time)
plt.plot(data_time, pred_B_ph, c = 'g')
plt.plot(data_time, data_B_ph, 'o', c = 'g', label = r'$\beta$-Phenyl')
plt.xlabel('Time')
plt.ylabel('Subj. intensity')
plt.legend()
plt.title('Data and Regression Curves')
plt.savefig('FigReg.png')
plt.show()
sb.kdeplot(a_smpls_Iso,
c = 'orange', label = 'Isovaleric(MAP={0:.3f})'.format(a_hat_Iso))
sb.kdeplot(a_smpls_B_ph,
c = 'g', label = r'$\beta$-Phenyl(MAP={0:.3f})'.format(a_hat_B_ph))
plt.legend()
plt.xlabel('a', fontsize
= 14)
plt.title('Posterior Distributions')
plt.savefig('FigDistri.png')
plt.show()
print('\na_hat_Iso
=', a_hat_Iso)
print('a_hat_B_ph
=', a_hat_B_ph)