Up

匂いの順応曲線

ベイズ回帰分析例

2025.11：CmdStanPy対応に改訂

　回帰分析（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（PyStan2）のPythonにおける使用例は、岡本（2019）で説明している。

リスト１のStanスクリプトを利用して図１のデータをベイズ分析するPythonスクリプトをリスト２のように用意した。

リスト１のファイル、リスト２のファイルと図１のデータファイルは、圧縮ファイルfiles.zipとしてまとめた。ダウンロード・展開すれば自由に利用できる。

リスト１のStanスクリプトファイル（reg_exp.stan）、リスト２のPythonスクリプトファイル（AdapCurve.py）、およびデータファイル（Time_course.xls）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされている環境において、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。

(stan) *****/files$ python AdapCurve.py

上のコマンドを実行すると、ファイルが読み込まれ、stanスクリプトのコンパイルが始まる。コンパイルは、多少時間が掛かる。コンパイル後、まずイソ吉草酸のデータに対するMCMCサンプリングが始まる。サンプリングが終了すると、トレース図が表示される（図２）。

図２　イソ吉草酸の場合

図２のWindowを閉じると、次に、β―フェニールエチルアルコールのデータに対するMCMCサンプリングが始まる。サンプリングが終了するとトレース図が表示される。

図３　β―フェニールエチルアルコールの場合

図３のWindowを閉じると、図４のグラフが表示される。

図４

イソ吉草酸とβ―フェニールエチルアルコールに対する回帰曲線（式（１））のグラフである。パラメータaの値は、事後分布の中央値である。分布の代表値として、Gelmanら（2021）は、平均値より安定した中央値を勧めている。

図４のWindowを閉じると、パラメータaの事後分布のグラフと中央値が表示される（図５）。

図５

時刻

イソ吉草酸

Βーフェニールエチルアルコール

tにおける匂いの強さに対する時間後の強さの比は次式で表される。

すなわち、時間後の強さの減少比率は時刻ｔに関わらず一定である。

イソ吉草酸の場合、

であるので、

となる。１分経過ごとに86％に減少する。

β―フェニールエチルアルコールの場合、

であるので、

となる。１分経過ごとに97％に減少する。

引用文献

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University

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;

array[N] real x;

array[N] real y;

}

parameters {

real<upper = 0.0> a;

real<lower = 0.0> sa;

real<lower = 0.0> sb;

}

transformed parameters {

array[N] real<lower = 0.0> sgm;

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

from cmdstanpy import CmdStanModel

import seaborn as sb

import arviz as az

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.array(df_data, dtype=float)

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)

model = CmdStanModel(stan_file='reg_exp.stan')

stan_data = {'N':len(data_time), 'x': data_time, 'y':data_Iso}

fit_Iso = model.sample(data=stan_data, adapt_delta=0.9999,

iter_warmup=10_000)

print(fit_Iso.diagnose())

print(fit_Iso.summary())

inf_Iso_data = az.from_cmdstanpy(fit_Iso) # Arviz DataFrame型に変換

az.plot_trace(inf_Iso_data)

plt.tight_layout()

plt.show()

fit_Iso = fit_Iso.draws_pd() # Pandas DataFrame型に変換

stan_data = {'N':len(data_time), 'x': data_time, 'y':data_B_ph}

fit_B_ph = model.sample(data=stan_data, adapt_delta=0.999)

print(fit_B_ph.diagnose())

print(fit_B_ph.summary())

inf_B_ph_data = az.from_cmdstanpy(fit_B_ph) # Arviz DataFrame型に変換

ax = az.plot_trace(inf_B_ph_data)

plt.tight_layout()

plt.show()

fit_B_ph = fit_B_ph.draws_pd() # Pandas DataFrame型に変換

a_smpls_Iso = fit_Iso['a']

a_hat_Iso = np.median(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')

a_smpls_B_ph = fit_B_ph['a']

a_hat_B_ph = np.median(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(Med.={0:.3f})'.format(a_hat_Iso))

sb.kdeplot(a_smpls_B_ph, c = 'g',

label = r'$\beta$-Phenyl(Med.={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)

Up