反応時間の分析
Shifted Gamma, ex-Gaussian, ex-Wald分布によるベイズ分析
反応時間の分布に対して、shifted Gamma分布、ex-Gaussian (exponentially modified Gaussian)分布、ex−Wald分布の当てはまりがよかったとPalmerら（Palmer, Horowitz, Torralba, & Wolfe,2011, J Exp Psychol Hum Percept Perform, 37, 58-71）は報告している。本ウェブサイトでは、以下にshifted Gamma分布によるベイズ分析用Stanスクリプト例、ex-Gaussian分布によるベイズ分析用Stanスクリプト例、ex-Wald分布によるベイズ分析用Stanスクリプト例を用意した。Python用StanであるPyStanについては、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。反応時間の実験プログラム例は別のウェブサイトに用意した。
なお、Okamoto (1981, 1982, 1980) は、単純反応時間において２状態の区別を報告している。反応時間のメカニズムの分析では状態の区別が必要であるが、本ウェブサイトではこの区別を行わず、単純に１つの分布を当てはめている。
本ウェブサイトにおける反応時間データ形式は共通であるので、ここで説明しておく。データファイルはCSV形式で用意する。データの形式は図１のようである。
・
・
・
図１
第１行目に変数名を設定し、データは２行目から書く。第１列目は各データの識別文字（図１では通し番号が使われている）、第２列目に反応時間が設定されている。データ入力後、CSV形式のファイルとして保存するが、これはファイル拡張子として「.csv」を選べばよい。本ウェブサイトの例では、「Data.csv」というファイル名である。
Shifted
Gamma分布による分析
確率分布shifted Gammaの確率密度関数は次式で与えられる。
上式の確率分布関数に対応してStanスクリプトをリストA.1のように用意した。
リストA.1 Shifted Gamma分布用Stanスクリプト（ファイル名shiftedGamma.stan）
data {
int N;
real RT[N];
real max_a;
real max_b;
real max_x;
}
parameters {
real<lower = 0.0, upper =
max_a> param_a;
real<lower = 0.0, upper =
max_b>param_b;
real<lower = 0.0, upper =
max_x> param_x;
}
model {
param_a ~ uniform(0.0,
max_a);
param_b ~ uniform(0.0,
max_b);
param_x ~ uniform(0.0,
max_x);
for(i in 1:N){
RT[i] - param_x ~ gamma(param_a, param_b);
}
}
リストA.1のStanスクリプトを実行するPythonスクリプトをリストA.2のように用意した。
リストA.2 Shifted Gamma分布による分析用Pythonスクリプト（ファイル名Main.py）
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import pystan
import pickle
#
#
データの読み込み
#
flnm = input('Data file name(*.csv) =
')
fin = open(flnm, 'r')
data = fin.readlines()
fin.close()
N = len(data) - 1
DataID = []
RTs = []
for i in range(N):
str = data[i+1].split(',')
DataID.append(str[0])
RTs.append(float(str[1]))
for i in range(N):
print(DataID[i],' ', RTs[i])
RTs = np.array(RTs)
xi_hat = RTs.min()
# シフト値の最大値
print('xi_hat = ', xi_hat)
RTs_xi = RTs - xi_hat
#
# モーメント法に依るパラメータ値の推定
#
beta_hat = RTs_xi.mean() /
RTs_xi.var()
print('beta_hat = ', beta_hat)
alpha_hat = beta_hat * RTs_xi.mean()
print('alpha_hat = ', alpha_hat)
max_a = alpha_hat * 100.0
max_b = beta_hat * 100.0
max_x = xi_hat
Data = {'N': len(RTs), 'RT': RTs,
'max_a': max_a, 'max_b': max_b, 'max_x': max_x}
sm = pystan.StanModel(file =
'shiftedGamma.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
param_a = fit['param_a']
param_b = fit['param_b']
param_x = fit['param_x']
plt.hist(param_a)
plt.xlabel(r'$\alpha$', fontsize = 14)
plt.title(r'Posterior Distribution of
$\alpha$', fontsize = 18)
plt.show()
plt.hist(param_b)
plt.xlabel(r'$\beta$', fontsize = 14)
plt.title(r'Posterior Distribution of
$\beta$', fontsize = 18)
plt.show()
plt.hist(param_x)
plt.xlabel(r'$\xi$', fontsize = 14)
plt.title(r'Posterior Distribution of
$\xi$', fontsize = 18)
plt.show()
#
# 中央値
#
med_a = np.percentile(param_a, 50)
med_b = np.percentile(param_b, 50)
med_xi = np.percentile(param_x, 50)
print('med_a = ', med_a, ' med_b = ', med_b, ' med_xi = ', med_xi)
min_x = RTs.min() - (RTs.max() -
RTs.min()) * 0.1
max_x = RTs.max() + (RTs.max() -
RTs.min()) * 0.1
rt = np.arange(min_x, max_x, (max_x -
min_x) / 1000)
#
# ガンマ分布
#
pdf_values = ss.gamma.pdf(rt, med_a,
med_xi, 1/med_b)
plt.plot(rt, pdf_values, linewidth =
3, label = 'gamma(ξ,α,β)')
#
# データのヒストグラム
#
plt.hist(RTs, density = True, label =
'RT')
plt.xlabel('Time', fontsize = 14)
plt.title('Shifted Gamma
Distribution\n' +
'ξ =
{0:<.3f}, α = {1:<.3f}, β
= {2:<.3f}\n'.format(med_xi, med_a, med_b),
fontsize = 16)
plt.legend(fontsize = 14)
plt.tight_layout()
plt.show()
リストA.1、リストA.2およびデータファイルData.csvは圧縮ファイルRT_shift_gamma.zipとしてまとめたので、ダウンロード解凍すれば利用できる。
リストA.2のスクリプトMain.pyを実行すると入力データファイル名の設定が求められる（図A.1）。
図A.1
入力データファイル名を設定すると、計算が始まる。Stanによるサンプリングが終わると、パラメータαの事後分布が表示される。
図A.2
図A.2のフォームを閉じると、パラメータβの事後分布が表示される（図A.3）。
図A.3
図A.3のフォームを閉じると、パラメータξの事後分布が表示される（図A.4）。
図A.4
図A.4のフォームを閉じると、反応時間のデータのヒストグラムと、事後分布における中央値として得られたパラメータの点推定値によるshifted Gamma分布のグラフが描かれたフォームが提示される（図A.5）。
図A.5
図A.5を閉じるとプログラムの終了である。
Ex-Gaussian分布による分析
確率分布ex-Gaussianの確率密度関数は次式で与えられる（Palmer et al., 2011）。
パラメータの変換
と次式の関係
を用いるとStanでの次式を得る。
さらに、パラメータの変換
と変数変換
を行うとscipy.statsでの次式を得る。
Ex-Gaussian分布用のStanスクリプトをリストB.1のように用意した。
リストB.1 Ex-Gaussian分布用Stanスクリプト（ファイル名Ex_Gaussian.stan）
functions {
real half_normal_lpdf(real
y, real sgm){
if
(y > 0.0){
return normal_lpdf(y | 0.0, sgm);
}
else {
return log(0.0);
}
}
}
data {
int N;
real RT[N];
}
parameters {
real mu;
real<lower = 0.0> sgm;
real<lower = 0.0> lmd;
}
model {
mu ~ normal(0.0, 1000);
sgm ~ half_normal(1000);
lmd ~ half_normal(1000);
for(i in 1:N){
RT[i] ~ exp_mod_normal(mu,
sgm, lmd);
}
}
リストB.1のStanスクリプトを利用するPythonスクリプトをリストB.2のように用意した。リストB.1、リストB.2およびサンプルデータファイルを圧縮ファイルRT_Ex_Gaussian.zipとしてまとめたので、ダウンロード解凍すれば利用することができる。
リストB.2 リストB.1のStanスクリプトを利用するPythonスクリプト（ファイル名Main.py）
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import pystan
import pickle
#
#
データの読み込み
#
flnm = input('Data file name(*.csv) =
')
fin = open(flnm, 'r')
data = fin.readlines()
fin.close()
N = len(data) - 1
DataID = []
RTs = []
for i in range(N):
str = data[i+1].split(',')
DataID.append(str[0])
RTs.append(float(str[1]))
for i in range(N):
print(DataID[i],' ', RTs[i])
RTs = np.array(RTs)
sm = pystan.StanModel(file =
'Ex_Gaussian.stan')
Data = {'N': N, 'RT': RTs}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
mu = fit['mu']
sgm = fit['sgm']
lmd = fit['lmd']
plt.hist(mu)
plt.title('Posterior Distribution of
$\mu$', fontsize = 16)
plt.xlabel('$\mu$', fontsize = 14)
plt.show()
plt.hist(sgm)
plt.title('Posterior Distribution of
$\sigma$', fontsize = 16)
plt.xlabel('$\sigma$', fontsize = 14)
plt.show()
plt.hist(lmd)
plt.title('Posterior Distribution of
$\lambda$', fontsize = 16)
plt.xlabel('$\lambda$', fontsize = 14)
plt.show()
mu_med = np.percentile(mu, 50)
sgm_med = np.percentile(sgm, 50)
lmd_med = np.percentile(lmd, 50)
K = 1.0 / (sgm_med * lmd_med)
x_min = RTs.min() - (RTs.max() -
RTs.min()) * 0.1
x_max = RTs.max() + (RTs.max() -
RTs.min()) * 0.1
x = np.arange(x_min, x_max, (x_max -
x_min) / 100)
RT_crv = ss.exponnorm.pdf(x, K,
mu_med, sgm_med)
plt.plot(x, RT_crv, linewidth = 3,
label = 'Ex-Gaussian')
plt.title('ex-Gaussian($\mu$,$\sigma$,$\lambda$)\n'
+\
'$\mu$ = {0:<.3f}
$\sigma$ = {1:<.5f}
$\lambda$ = {2:<.3f}'.
format(mu_med, sgm_med, lmd_med), fontsize = 16)
plt.hist(RTs, density = True, label =
'RT')
plt.xlabel('Time', fontsize = 14)
plt.legend()
plt.show()
リストB.2のスクリプトを実行すると、入力データファイル名の設定が求められる。
図B.1
図１の形式のデータのファイル名を設定すると計算が始まる。Stanによるサンプリングが終わると、まずパラメータμの事後分布が表示される（図B.2）。
図B.2
図B.2のフォームを閉じると、パラメータσの事後分布のグラフが表示される（図B.3）。
図B.3
図B.3のフォームを閉じると、パラメータλの事後分布のグラフが表示される（図B.4）。
図B.4
図B.4のフォームを閉じると、データのヒストグラムと事後分布の中央値をパラメータの点推定値とするex-Gaussian分布のグラフが表示される（図B.5）。
図B.5
図B.5のフォームを閉じるとプログラムの実行終了となる。
Ex-Wald分布による分析
確率分布ex-Wald分布の確率密度関数は次式で与えられる。
上式において、パラメータ
、
、
を
倍しても値は変わらない。すなわち、
である。このパラメータ値の不定性を除くため、
と選び、
とおいても、一般性は失われない。プログラミングはこのの設定の下で行った。Stanスクリプト例をリストC.1に示す。
リストC.1 Ex-Wald分布用Stanスクリプト例（ファイル名Ex_Wald.stan）
//
// Yasuharu
Okamoto, 2019.10
//
functions {
real half_normal_lpdf(real
y, real sgm){
if
(y > 0.0){
return normal_lpdf(y | 0.0, sgm);
}
else {
return
log(0.0);
}
}
real Ex_Wald_lpdf(real t,
real mu, real a, real g){
real
g_exp;
real
F;
real
k1;
k1 =
sqrt(mu*mu - 2.0 * g);
g_exp = g * exp(-g * t + a * (mu - k1));
F =
normal_cdf((k1 * t - a)/sqrt(t), 0.0, 1.0) +
exp(2.0 * a * k1) *
normal_cdf(-(k1 * t + a)/ sqrt(t), 0.0, 1.0);
return log(g_exp * F);
}
}
data {
int N;
real RT[N];
}
parameters {
real<lower = 0.0> k;
real<lower = 0.0> a;
real<lower = 0.0> g;
}
transformed parameters {
real<lower = 0.0> mu;
mu = sqrt(k*k + 2.0*g);
}
model {
k ~ normal(0.0, 1000);
a ~ half_normal(1000);
g ~ half_normal(1000);
for(i in 1:N){
RT[i] ~ Ex_Wald(mu, a, g);
}
}
リストC.1のStanスクリプトを利用するPythonスクリプトを、リストC.2に示す。リストC1およびリストC.2のファイルとサンプルデータファイルは、圧縮ファイルRT_Ex_Wald.zipとしてまとめた。ダウンロード解凍すれば利用できる。
リストC.2 リストC.1のStanスクリプトを利用するPythonスクリプト（ファイル名Main.py）
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import pystan
import pickle
#
#
データの読み込み
#
flnm = input('Data file name(*.csv) =
')
fin = open(flnm, 'r')
data = fin.readlines()
fin.close()
N = len(data) - 1
DataID = []
RTs = []
for i in range(N):
str = data[i+1].split(',')
DataID.append(str[0])
RTs.append(float(str[1]))
for i in range(N):
print(DataID[i],' ', RTs[i])
RTs = np.array(RTs)
sm = pystan.StanModel(file =
'Ex_Wald.stan')
Data = {'N': N, 'RT': RTs}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
k = fit['k']
a = fit['a']
g = fit['g']
mu = fit['mu']
plt.hist(a)
plt.title(r'Posterior Distribution of
$\alpha$', fontsize = 18)
plt.xlabel(r'$\alpha$', fontsize = 14)
plt.show()
plt.hist(g)
plt.title(r'Posterior Distribution of
$\gamma$', fontsize = 18)
plt.xlabel(r'$\gamma$', fontsize = 14)
plt.show()
plt.hist(mu)
plt.title(r'Posterior Distribution of
$\mu$', fontsize = 18)
plt.xlabel(r'$\mu$', fontsize = 14)
plt.show()
k_med = np.percentile(k, 50)
a_med = np.percentile(a, 50)
g_med = np.percentile(g, 50)
mu_med = np.percentile(mu, 50)
import math
def ex_wald( t, mu, a, g ):
w = []
for v in t:
if v
<= 0.0:
w.append(0.0)
else:
k = (mu**2
- 2*g) ** 0.5
g_exp = g * math.exp(-g * v + (a * (mu - k)))
F = ss.norm.cdf((k * v - a) / (v**0.5)) +\
math.exp(2.0 * a * k) *\
ss.norm.cdf(-(k * v + a) / (v**0.5))
w.append(g_exp * F)
return np.array(w)
L_x = RTs.min() - (RTs.max() -
RTs.min()) * 0.1
U_x = RTs.max() + (RTs.max() -
RTs.min()) * 0.1
x = np.arange(L_x, U_x, 0.001)
y = ex_wald(x, mu_med, a_med, g_med)
plt.plot(x, y, linewidth = 3, label =
'ex-Wald')
plt.hist(RTs, bins = 10, density =
True, label = 'RT')
plt.xlabel('Time', fontsize = 14)
plt.title(r'ex-Wald($\mu$,$\sigma$=1,$\alpha$,$\gamma$)'
+ '\n' +\
(r'$\mu$ = {0:<.3f}
$\alpha$ = {1:<.3f}
$\gamma$ = {2:<.3f}'.
format(mu_med, a_med, g_med)), fontsize = 16)
plt.legend(fontsize = 14)
plt.show()
リストC.2のプログラムMain.pyを実行すると、入力データファイル名の設定が求められる（図C.1）。
図C.1
図１の形式のデータファイル名を設定すると計算が始まる。Stanによるサンプリングが終了すると、パラメータαの事後分布が表示される（図C.2）。
図C.2
図C.2のフォームを閉じると、パラメータγの事後分布が表示される（図C.3）。
図C.3
図C.3のフォームを閉じると、パラメータμの事後分布が表示される（図C.4）。
図C.4
図C.4のフォームを閉じると、データのヒストグラムと事後分布の中央値をパラメータの点推定値とするex-Wald分布のグラフが表示される（図C.5）。
図C.5
図C.5のフォームを閉じると、プログラムの実行終了となる。