反応時間の分布
Ex-Gaussian分布、対数正規分布、３パラメータ対数正規分布
刺激に対する反応が複数のプロセスによって行われるとき、それらのプロセスの処理時間が反応時間に反映すると考えると、実験条件の違いにおける反応時間の差はプロセスの実験条件の違いにおける処理時間の違いを反映していると考えられ、心理学研究における有用な指標の1つとして用いられる。反応時間に関わる心理学的情報処理モデルはいろいろ提案されてきている（上田・椎名・浅川、2003）。反応時間の分析において、心理学的プロセスを具体的に設定するのではなく、反応時間の記述モデルにより分析することも広く行われている。対数正規分布が、上田ら（2003）、伊関（2020）で論じられており、指数付き正規分布（Exponentially modified Gaussian distribution: Ex-Gaussian distribution）が、伊関（2020）、三ヶ尻（2017）で論じられている。三ヶ尻（2017）では、Ex-Gaussianが不適切である場合のあることを示すグラフが提示されている。また、三ヶ尻（2017）では、対数正規分布による分析結果は示されていない。Ex-Gaussian分布および対数正規分布は、心理学における反応時間の分析における有用性の評価を慎重に行う必要があると考え、いくつかのグラフを選び、Ex-Gaussian分布、対数正規分布による再分析を行った。確率分布としては、３パラメータ対数正規分布も取り上げたが、これは反応時間がある値以上の値をとるときに対数正規分布より適切と考えられるためである。３パラメータ対数正規分布による平均値の差の分析例は、このウェブサイトに用意されている。また、比較のために再分析においては正規分布による分析結果も示した。
再分析の対象とするグラフは、三ヶ尻（2017）においてEx-Gaussian分布の当てはまりの良いものとして図3.1、悪いものとして図3.8および図3.10を選んで再分析した。グラフの値は、このウェブサイトのプログラムを再分析用に改編したものを利用した。いずれの再分析においても、データへの適合は３パラメータ対数正規分布（Johnson, Kotz, & Barakrishnan, 1994）がもっともよかった。
ここで、
Ex-Gaussian分布に従う確率変数
は、正規分布（normal distribution; Gaussian distribution）に従う確率変数
と指数分布に従う確率変数
の和
として表すことができる。とは独立である。
（２パラメータ）対数正規分布に従う確率変数は、正規分布に従う確率変数の指数関数
として表すことができる。
３パラメータ対数正規分布に従う確率変数は、対数正規分布に従う確率変数と定数
との和
として表すことができる。
なお、Dual System Model（二重砂山モデル）が、三ヶ尻（2017）で提案されている。Dual System Modelと３パラメータ対数正規分布の比較も行った。
図3.1（三ヶ尻、2017）の再分析を、リスト１のスクリプトにより行い、Figure 1の結果を得た。
Figure 1 図3.1の再分析
図3.1から読み取った該当区間における確率密度を縦棒グラフで示した。それぞれの分布のパラメータは最小2乗法により求めた。
データの分布が非対称であるので、緑の点線で示されている正規分布の当てはまりは良くない。橙色の曲線で示されている対数正規分布のグラフは、正規分布よりはデータへの当てはまりが良いが右裾の分布への当てはまりは良くない。赤色の1点破線で描かれているEx-Gaussian分布は全体としてデータの分布に良く当てはまっているように見えるが、モードでの値が高すぎる。これらに対して緑色破線で示されている３パラメータ対数正規分布は、全体としてデータに良く当てはまっているようである。
図3.8（三ヶ尻、2017）の再分析を、リスト２のスクリプトにより行い、Figure２の結果を得た。
Figure ２ 図3.８の再分析
図3.８から読み取った該当区間における確率密度を縦棒グラフで示した。それぞれの分布のパラメータは最小2乗法により求めた。
データの分布が非対称であるので、緑の点線で示されている正規分布の当てはまりは良くない。橙色の曲線で示されている対数正規分布のグラフは、正規分布よりはデータへの当てはまりが良いが右裾の分布への当てはまりは良くない。赤色の1点破線で描かれているEx-Gaussian分布は全体としてデータの分布に良く当てはまっているように見えるが、モードでの値が高すぎる。これらに対して緑色破線で示されている３パラメータ対数正規分布は全体としてデータに良く当てはまっているようである。
図3.10（三ヶ尻、2017）の再分析を、リスト３のスクリプトにより行い、Figure３の結果を得た。
Figure ３ 図3.10の再分析
図3.10から読み取った該当区間における確率密度を縦棒グラフで示した。それぞれの分布のパラメータは最小2乗法により求めた。
データの分布が非対称であるので、緑の点線で示されている正規分布の当てはまりは良くない。橙色の曲線で示されている対数正規分布のグラフは、正規分布よりはデータへの当てはまりが良いが右裾の分布への当てはまりは良くない。赤色の1点破線で描かれているEx-Gaussian分布は左裾が短く急峻であり、モードでの値が突出している。緑色破線で示されている３パラメータ対数正規分布は対数正規分布に似ているが、より全体をよく表しているようである。
全体として、３パラメータ対数正規分がよさそうである。３パラメータ対数正規分布は、ベイズ分析において対数正規分布と同じように用いることができることが、このウェブサイトに示されている。
本ウェブサイトの分析では、Ex-Gaussian分布の適用も試みられたが、最小2乗基準を用いた極値探索法では、関数exp()におけるオーバーフローを避けるために独立変数の単位を変える必要があった。スクリプト中の次の部分がそれに該当する。
X_eg = X / 1000
# 極値探索において関数exp()の
# オーバーフローを避けるためのスケーリング
Ydensity_eg = Ydensity
* 1000 # Xのスケーリングに伴う調整
Ex-Gaussian分布を用いるときは、注意が必要である。
３パラメータ対数正規分布とDual
Systemモデルの比較
３パラメータ対数正規分布とDual Systemモデルの比較を行う。
Dual Systemモデルでは、２つのsystem、System AとSystem B、考え、２つのsystemの処理が終わった時に反応が生起すると考えることができる。２つのsystemの処理が独立で処理時間が指数分布に従うとすれば、処理がt時間内に終わる確率は次式で与えられる。
は反応時間の閾値であり、式番号は、三ヶ尻（2017）におけるものである。
したがって、確率密度関数は次式で与えられる。
式番号(6.5rev)は、三ヶ尻の式(6.5)を修正したものであることを表す。
Dual bSystemモデルと３パラメータ対数正規分布の比較を、三ヶ尻（2017）の図6.3（実験参加者１のデータ）と図6.4（実験参加者２のデータ）に対して行った。グラフの値は、このウェブサイトのプログラムを変更して用いた。
図6.3のグラフの再分析を、リスト４により行った。結果をFigure 4に示す。
Figure 4 図6.3の再分析
Dual Systemモデルのパラメータは、三ヶ尻（2017）の表6.1の値を用いたものを緑の破線で、最小2乗の基準で極値探索法（Powell法）により求めた値によるものを青の1点鎖線で示した。表6.1の値による緑色の破線はデータに対して少し左寄りであるが、Powell法による値の青の1点鎖線と３パラメータ対数正規分布のグラフはほぼ重なっていて、データの分布をよく表していると言える。
図6.4のグラフの再分析を、リスト５により行った。結果をFigure 5に示す。
Figure 5 図6.4の再分析
Dual Systemモデルのパラメータは、三ヶ尻（2017）の表6.1の値を用いたものを緑の破線で、最小2乗の基準で極値探索法（Powell法）により求めた値によるものを青の1点鎖線で示した。表6.1の値による緑色の破線はデータに対して少し左寄りであるが、Powell法による値の青の1点鎖線と３パラメータ対数正規分布のグラフはほぼ重なっていて、データの分布をよく表していると言える。
井関龍太（2020）心理学者は反応時間をどう分析するか． 基礎心理学研究、38、243ー249．
Johnson, N. L., Kotz, S., & Barakrishnan, N. (1994) Continuous univariate distributions, Vol.1, 2nd Ed. John Wiley & Sons, Inc.
三ヶ尻陽一（2017）反応が競合する場合の反応時間分布と誤りを表現するモデルとその解釈．東京大学学位論文（12601甲第33686号）
上田卓司・椎名乾平・浅川伸一（2003）反応時間の確率モデル． 心理学評論、46、249ー273．
附記
リスト１ 図3.1の再分析スクリプト（F3_1.zip）
import pandas as pd
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.optimize
as so
import scipy.stats
as ss
pddata = pd.read_excel('MikajiriF3_1adj.xlsx')
print(pddata.keys())
X = pddata['Xadj']
Y = pddata['Yadj']
X = np.array(X)
Y = np.array(Y)
Ydensity = (Y / np.sum(Y))
/ 100
"""
対数正規分布による分析
"""
def lognormal_pdf(x,
mu, sgm):
v1 = sgm
* x * ((2*np.pi)**0.5)
v2 = np.exp(-((np.log(x)
- mu)**2) / (2 * (sgm**2)))
return v2 / v1
class sqr_diff():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
"""
最小２乗法を用いるための目的関数
"""
mu =
params[0]
sgm = 0.0001 + params[1]**2
v =
0.0
v +=
np.sum((self.y - lognormal_pdf(self.x, mu, sgm))**2)
return v
s_d = sqr_diff(X,
Ydensity)
x0 = [0, 1]
results = so.minimize(s_d.f, x0, method='Powell')
print(results)
mu = results.x[0]
sgm = 0.0001 + results.x[1]**2
xcoord = np.linspace(1,
3000, 1000)
y_log_norm = lognormal_pdf(xcoord, mu, sgm)
"""
３パラメータ対数正規分布
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
class sqr_diff_shift():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
shift = params[2]
v =
0.0
for i in range(self.n):
v += (self.y[i] -
lognormal_shift_pdf(self.x[i], mu, sgm, shift))**2
return v
s_d_t = sqr_diff_shift(X,
Ydensity)
x0 = [0, 1, 0]
results_t = so.minimize(s_d_t.f, x0, method='Powell')
print(results_t)
mu_t = results_t.x[0]
sgm_t = 0.0001 + results_t.x[1]**2
shift = results_t.x[2]
y_log_norm_t = []
for vx in xcoord:
y_log_norm_t.append(lognormal_shift_pdf(vx, mu_t, sgm_t, shift))
"""
指数付き正規分布
Exponentially Modified Gaussian distribution
"""
def exp_gaussian(t,
mu, sgm,tau):
v = (1/tau)*np.exp(-(t/tau) + (mu/tau) + ((sgm**2)/(2
* tau**2))) * \
ss.Normal(mu=0, sigma=1).cdf((t -
mu - ((sgm**2)/tau)) / sgm)
return v
class sqr_diff_exG():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
tau
= 0.01 + params[2]**2
v =
0.0
for i in range(self.n):
v += (self.y[i] - exp_gaussian(self.x[i], mu, sgm, tau))**2
return v
X_eg = X / 1000
# 極値探索において関数exp()の
# オーバーフローを避けるためのスケーリング
Ydensity_eg = Ydensity
* 1000 # Xのスケーリングに伴う調整
sqr_exg = sqr_diff_exG(X_eg, Ydensity_eg)
x0 = [10, 1, 1]
results_exg = so.minimize(sqr_exg.f, x0, method='Powell')
print(results_exg)
mu_exg = results_exg.x[0]
sgm_exg = 0.0001 + results_exg.x[1]**2
tau_exg = 0.01 + results_exg.x[2]**2
xcoord_eg = xcoord/1000 # Xのスケーリングに伴う調整
y_ex_g = exp_gaussian(xcoord_eg, mu_exg, sgm_exg, tau_exg) / 1000
"""
正規分布
"""
class sqr_diff_G():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
v = ss.Normal(mu=mu, sigma=sgm).pdf(self.x)
return np.sum((self.y
- v)**2)
sqr_g = sqr_diff_G(X,
Ydensity)
x0 = [500, 10]
results_g = so.minimize(sqr_g.f, x0, method='Powell')
print(results_g)
mu_g = results_g.x[0]
sgm_g = 0.0001 + results_g.x[1]**2
y_g = ss.Normal(mu=mu_g, sigma=sgm_g).pdf(xcoord)
plt.title('Data fron
Figure 3.1', fontsize=16)
plt.bar(X, Ydensity,
width=80, alpha=0.5, label='Data')
plt.plot(xcoord, y_log_norm_t, label='3 param lognormal', ls='--', lw=5)
plt.plot(xcoord, y_log_norm, label='lognormal', lw=3)
plt.plot(xcoord, y_g, label='Gausian', ls=':', lw=5)
plt.plot(xcoord, y_ex_g, label='Ex-Gaussian', ls='-.', lw=3)
plt.legend(handlelength=5)
plt.show()
リスト２ 図3.8の再分析用スクリプト（F3_8.zip）
import pandas as pd
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.optimize
as so
import scipy.stats
as ss
pddata = pd.read_excel('MikajiriF3_8adj.xlsx')
print(pddata.keys())
X = pddata['Xadj']
Y = pddata['Yadj']
X = np.array(X)
Y = np.array(Y)
Ydensity = (Y / np.sum(Y))
/ 100
"""
対数正規分布による分析
"""
def lognormal_pdf(x,
mu, sgm):
v1 = sgm
* x * ((2*np.pi)**0.5)
v2 = np.exp(-((np.log(x)
- mu)**2) / (2 * (sgm**2)))
return v2 / v1
class sqr_diff():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
"""
最小２乗法を用いるための目的関数
"""
mu =
params[0]
sgm = 0.0001 + params[1]**2
v =
0.0
v +=
np.sum((self.y - lognormal_pdf(self.x, mu, sgm))**2)
return v
s_d = sqr_diff(X,
Ydensity)
x0 = [0, 1]
results = so.minimize(s_d.f, x0, method='Powell')
print(results)
mu = results.x[0]
sgm = 0.0001 + results.x[1]**2
xcoord = np.linspace(1,
3000, 1000)
y_log_norm = lognormal_pdf(xcoord, mu, sgm)
"""
３パラメータ対数正規分布
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
class sqr_diff_shift():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
shift = params[2]
v =
0.0
for i in range(self.n):
v += (self.y[i] -
lognormal_shift_pdf(self.x[i], mu, sgm, shift))**2
return v
s_d_t = sqr_diff_shift(X,
Ydensity)
x0 = [0, 1, 0]
results_t = so.minimize(s_d_t.f, x0, method='Powell')
print(results_t)
mu_t = results_t.x[0]
sgm_t = 0.0001 + results_t.x[1]**2
shift = results_t.x[2]
y_log_norm_t = []
for vx in xcoord:
y_log_norm_t.append(lognormal_shift_pdf(vx, mu_t, sgm_t, shift))
"""
指数付き正規分布
Exponentially Modified Gaussian distribution
"""
def exp_gaussian(t,
mu, sgm,tau):
v = (1/tau)*np.exp(-(t/tau) + (mu/tau) + ((sgm**2)/(2
* tau**2))) * \
ss.Normal(mu=0, sigma=1).cdf((t -
mu - ((sgm**2)/tau)) / sgm)
return v
class sqr_diff_exG():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
tau
= 0.0001 + params[2]**2
v =
0.0
for i in range(self.n):
v += (self.y[i] - exp_gaussian(self.x[i], mu, sgm, tau))**2
return v
X_eg = X/1000
# 極値探索において関数exp()の
# オーバーフローを避けるためのスケーリング
Ydensity_eg = Ydensity
* 1000 # Xのスケーリングに伴う調整
sqr_exg = sqr_diff_exG(X_eg, Ydensity_eg)
x0 = [0.1, 0.1, 10]
results_exg = so.minimize(sqr_exg.f, x0, method='Powell')
print(results_exg)
mu_exg = results_exg.x[0]
sgm_exg = 0.0001 + results_exg.x[1]**2
tau = 0.0001 + results_exg.x[2]**2
xcoord_eg = xcoord/1000 # Xのスケーリングに伴う調整
y_ex_g = exp_gaussian(xcoord_eg, mu_exg, sgm_exg, tau)/1000
"""
正規分布
"""
class sqr_diff_G():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
v = ss.Normal(mu=mu, sigma=sgm).pdf(self.x)
return np.sum((self.y
- v)**2)
sqr_g = sqr_diff_G(X,
Ydensity)
x0 = [500, 10]
results_g = so.minimize(sqr_g.f, x0, method='Powell')
print(results_g)
mu_g = results_g.x[0]
sgm_g = 0.0001 + results_g.x[1]**2
y_g = ss.Normal(mu=mu_g, sigma=sgm_g).pdf(xcoord)
plt.title('Data from Figure 3.8', fontsize=16)
plt.bar(X, Ydensity,
width=80, alpha=0.5, label='Data')
plt.plot(xcoord, y_log_norm_t, label='3 param lognormal', ls='--', lw=5)
plt.plot(xcoord, y_log_norm, label='lognormal', lw=3)
plt.plot(xcoord, y_g, label='Gausian', ls=':', lw=5)
plt.plot(xcoord, y_ex_g, label='Ex-Gaussian', ls='-.', lw=3)
plt.legend(handlelength=5)
plt.show()
リスト３ 図3.10の再分析用スクリプト（F3_10.zip）
import pandas as pd
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.optimize
as so
import scipy.stats
as ss
pddata = pd.read_excel('MikajiriF3_10adj.xlsx')
print(pddata.keys())
X = pddata['Xadj']
Y = pddata['Yadj']
X = np.array(X)
Y = np.array(Y)
Ydensity = (Y / np.sum(Y))
/ 100
"""
対数正規分布による分析
"""
def lognormal_pdf(x,
mu, sgm):
v1 = sgm
* x * ((2*np.pi)**0.5)
v2 = np.exp(-((np.log(x)
- mu)**2) / (2 * (sgm**2)))
return v2 / v1
class sqr_diff():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
"""
最小２乗法を用いるための目的関数
"""
mu =
params[0]
sgm = 0.0001 + params[1]**2
v =
0.0
v +=
np.sum((self.y - lognormal_pdf(self.x, mu, sgm))**2)
return v
s_d = sqr_diff(X,
Ydensity)
x0 = [0, 1]
results = so.minimize(s_d.f, x0, method='Powell')
print(results)
mu = results.x[0]
sgm = 0.0001 + results.x[1]**2
xcoord = np.linspace(1,
3000, 1000)
y_log_norm = lognormal_pdf(xcoord, mu, sgm)
"""
３パラメータ対数正規分布
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
class sqr_diff_shift():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
shift = params[2]
v =
0.0
for i in range(self.n):
v += (self.y[i] -
lognormal_shift_pdf(self.x[i], mu, sgm, shift))**2
return v
s_d_t = sqr_diff_shift(X,
Ydensity)
x0 = [0, 1, 0]
results_t = so.minimize(s_d_t.f, x0, method='Powell')
print(results_t)
mu_t = results_t.x[0]
sgm_t = 0.0001 + results_t.x[1]**2
shift = results_t.x[2]
y_log_norm_t = []
for vx in xcoord:
y_log_norm_t.append(lognormal_shift_pdf(vx, mu_t, sgm_t, shift))
"""
指数付き正規分布
Exponentially Modified Gaussian distribution
"""
def exp_gaussian(t,
mu, sgm,tau):
v = (1/tau)*np.exp(-(t/tau) + (mu/tau) + ((sgm**2)/(2
* tau**2))) * \
ss.Normal(mu=0, sigma=1).cdf((t -
mu - ((sgm**2)/tau)) / sgm)
return v
class sqr_diff_exG():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
tau
= 0.0001 + params[2]**2
v =
0.0
for i in range(self.n):
v += (self.y[i] - exp_gaussian(self.x[i], mu, sgm, tau))**2
return v
X_eg = X/1000
# 極値探索において関数exp()の
# オーバーフローを避けるためのスケーリング
Ydensity_eg = Ydensity
* 1000 # Xのスケーリングに伴う調整
sqr_exg = sqr_diff_exG(X_eg, Ydensity_eg)
x0 = [0.1, 0.1, 10]
results_exg = so.minimize(sqr_exg.f, x0, method='Powell')
print(results_exg)
mu_exg = results_exg.x[0]
sgm_exg = 0.0001 + results_exg.x[1]**2
tau = 0.0001 + results_exg.x[2]**2
xcoord_eg = xcoord/1000
# Xのスケーリングに伴う調整
y_ex_g = exp_gaussian(xcoord_eg, mu_exg, sgm_exg, tau)/1000
"""
正規分布
"""
class sqr_diff_G():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
v = ss.Normal(mu=mu, sigma=sgm).pdf(self.x)
return np.sum((self.y
- v)**2)
sqr_g = sqr_diff_G(X,
Ydensity)
x0 = [500, 10]
results_g = so.minimize(sqr_g.f, x0, method='Powell')
print(results_g)
mu_g = results_g.x[0]
sgm_g = 0.0001 + results_g.x[1]**2
y_g = ss.Normal(mu=mu_g, sigma=sgm_g).pdf(xcoord)
plt.title('Data from Figure 3.10', fontsize=16)
plt.bar(X, Ydensity,
width=80, alpha=0.5, label='Data')
plt.plot(xcoord, y_log_norm_t, label='3 param lognormal', ls='--', lw=5)
plt.plot(xcoord, y_log_norm, label='lognormal', lw=3)
plt.plot(xcoord, y_g, label='Gausian', ls=':', lw=5)
plt.plot(xcoord, y_ex_g, label='Ex-Gaussian', ls='-.', lw=5)
plt.legend(handlelength=5)
plt.show()
リスト４ 図6.3の再分析（Fig6_3.zip）
import pandas as pd
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.optimize
as so
import scipy.stats
as ss
pddata = pd.read_excel('Fig6_3adj.xlsx')
print(pddata.keys())
X = pddata['Xadj']
Y = pddata['Yadj']
X = np.array(X)
Y = np.array(Y)
Ydensity = (Y / np.sum(Y))
/ 100
xcoord = np.linspace(1,
3000, 1000)
"""
３パラメータ対数正規分布
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
class sqr_diff_shift():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
shift = params[2]
v =
0.0
for i in range(self.n):
v += (self.y[i] -
lognormal_shift_pdf(self.x[i], mu, sgm, shift))**2
return v
s_d_t = sqr_diff_shift(X,
Ydensity)
x0 = [0, 1, 0]
results_t = so.minimize(s_d_t.f, x0, method='Powell')
print(results_t)
mu_t = results_t.x[0]
sgm_t = 0.0001 + results_t.x[1]**2
shift = results_t.x[2]
y_log_norm_t = []
for vx in xcoord:
y_log_norm_t.append(lognormal_shift_pdf(vx, mu_t, sgm_t, shift))
"""
Dual System model
"""
a = 0.09
b = 0.32
tc = 0.56
def p_ds(t,
a, b, tc):
if t <= tc:
return 0.0
v = (1.0/b)*np.exp(-(t - tc) / b) * (1.0 - np.exp(-(t - tc) / a)) + \
(1.0/a)*np.exp(-(t - tc)
/ a) * (1.0 - np.exp(-(t - tc)
/ b))
return v
xcoord_ds = xcoord/1000 # 単位をミリ秒から秒へ
y_ds = []
for vx in xcoord_ds:
y_ds.append(p_ds(vx, a, b, tc))
y_ds = np.array(y_ds) / 1000
class DualSystem():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
self.x_min = np.min(x)
def f(self, params):
a =
0.0001 + params[0]**2
b =
0.0001 + params[1]**2
tc = params[2]
// self.x_min - params[2]**2
v =
0.0
for i in range(self.n):
v += (self.y[i] - p_ds(self.x[i],
a, b, tc))**2
return v
X_ds = X / 1000
ydensity_ds = Ydensity
* 1000
d_sys = DualSystem(X_ds, ydensity_ds)
x0 = [a**0.5, b**0.5, tc]
results_ls = so.minimize(d_sys.f, x0, method='Powell')
print(results_ls)
a_ls = 0.0001 + results_ls.x[0]**2
b_ls = 0.0001 + results_ls.x[1]**2
tc_ls = results_ls.x[2]
print(a_ls, b_ls, tc_ls)
y_ls = []
for vx in xcoord_ds:
y_ls.append(p_ds(vx, a_ls,
b_ls, tc_ls))
y_ls = np.array(y_ls) / 1000
plt.title('Data fron
Figure 6.3', fontsize=16)
plt.bar(X, Ydensity,
width=80, alpha=0.5, label='Data')
plt.plot(xcoord, y_ds, label='Dual System Table 6.1', ls='--', lw=3, c='g', alpha=0.7)
plt.plot(xcoord, y_ls, label='Dual System LS', ls='-.', lw=3,
c='b', alpha=0.7)
plt.plot(xcoord, y_log_norm_t, c='orange',
label='3 param lognormal', ls='--', lw=5)
plt.legend(handlelength=5)
plt.show()
リスト５ 図6.4の再分析（Fig6_4.zip）
import pandas as pd
import matplotlib.pyplot
as plt
import numpy
as np
import scipy.optimize
as so
import scipy.stats
as ss
pddata = pd.read_excel('Fig6_4adj.xlsx')
print(pddata.keys())
X = pddata['Xadj']
Y = pddata['Yadj']
X = np.array(X)
Y = np.array(Y)
Ydensity = (Y / np.sum(Y))
/ 100
xcoord = np.linspace(1,
3000, 1000)
"""
３パラメータ対数正規分布
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
class sqr_diff_shift():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
def f(self, params):
mu =
params[0]
sgm = 0.0001 + params[1]**2
shift = params[2]
v =
0.0
for i in range(self.n):
v += (self.y[i] -
lognormal_shift_pdf(self.x[i], mu, sgm, shift))**2
return v
s_d_t = sqr_diff_shift(X,
Ydensity)
x0 = [0, 1, 0]
results_t = so.minimize(s_d_t.f, x0, method='Powell')
print(results_t)
mu_t = results_t.x[0]
sgm_t = 0.0001 + results_t.x[1]**2
shift = results_t.x[2]
y_log_norm_t = []
for vx in xcoord:
y_log_norm_t.append(lognormal_shift_pdf(vx, mu_t, sgm_t, shift))
"""
Dual System model
"""
a = 0.16
b = 0.48
tc = 0.43
def p_ds(t,
a, b, tc):
if t <= tc:
return 0.0
v = (1.0/b)*np.exp(-(t - tc) / b) * (1.0 - np.exp(-(t - tc) / a)) + \
(1.0/a)*np.exp(-(t - tc)
/ a) * (1.0 - np.exp(-(t - tc)
/ b))
return v
xcoord_ds = xcoord/1000 # 単位をミリ秒から秒へ
y_ds = []
for vx in xcoord_ds:
y_ds.append(p_ds(vx, a, b, tc))
y_ds = np.array(y_ds) / 1000
class DualSystem():
def __init__(self,
x, y):
self.x = np.array(x)
self.n = len(x)
self.y = np.array(y)
self.x_min = np.min(x)
def f(self, params):
a =
0.0001 + params[0]**2
b =
0.0001 + params[1]**2
tc = params[2]
#// self.x_min - params[2]**2
v =
0.0
for i in range(self.n):
v += (self.y[i] - p_ds(self.x[i],
a, b, tc))**2
return v
X_ds = X / 1000
ydensity_ds = Ydensity
* 1000
d_sys = DualSystem(X_ds, ydensity_ds)
x0 = [a**0.5, b**0.5, tc]
results_ls = so.minimize(d_sys.f, x0, method='Powell')
print(results_ls)
a_ls = 0.0001 + results_ls.x[0]**2
b_ls = 0.0001 + results_ls.x[1]**2
tc_ls = results_ls.x[2]
print(a_ls, b_ls, tc_ls)
y_ls = []
for vx in xcoord_ds:
y_ls.append(p_ds(vx, a_ls,
b_ls, tc_ls))
y_ls = np.array(y_ls) / 1000
plt.title('Data fron
Figure 6.4', fontsize=16)
plt.bar(X, Ydensity,
width=80, alpha=0.5, label='Data')
plt.plot(xcoord, y_ds, label='Dual System', ls='--', lw=3,
c='g', alpha=0.7)
plt.plot(xcoord, y_ls, label='Dual System', ls='-.', lw=3,
c='b', alpha=0.7)
plt.plot(xcoord, y_log_norm_t, c='orange', alpha=0.8,
label='3 param lognormal', ls='--', lw=5)
plt.legend(handlelength=5)
plt.show()