マグニチュード推定法
ベイズ分析(CmdStanPy 対応)
CmdStanPyの準備の簡単な説明のウェブサイトを用意した。
プロセティックな感覚(量的な感覚)の場合、物理刺激量と感覚量
の間の関係を表す式として、べき法則と呼ばれている次式
あるいは、これをより一般化した次式
がある。ただし、式番号は<岡本安晴(2014)「心理学データ分析と測定」勁草書房>におけるものである。
心理物理学実験において、第試行で提示された物理刺激の物理量を
で、その感覚量のマグニチュード推定量(Magnitude Estimation: ME)を
で表すとき、以下のモデルを設定する(Stevens, 1975; p. 169)。
モデル(10.2.1)に対しては、
モデル(10.2.2)に対しては、
ここで、は平均
、分散
の正規分布を表す。
上記の確率モデルをデータに適用して、ベイズ分析によりパラメータ値を求める。ベイズ分析については、<岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版>の題3部、第9章〜第15章、で説明している。
まず、一般化されたモデル(10.2.4)による分析について説明する。その後、モデル(10.2.4)においてと置いた場合のモデル(10.2.3)による分析を示す。
なお、ファイルはmefiles.zipにまとめている。自由に利用されたい。
モデル(10.2.4)に基づいてベイズ分析を行うためのStanスクリプトをリストG1のように用意した。
モデル(10.2.2)では、3つのパラメータの間に強い連関が予想される。パラメータが小さくなったときは、パラメータ
の値を大きく、あるいはパラメータ
の値を負の方向に大きくすればモデル(10.2.2)とデータの関係の維持が期待される。このパラメータ間の連関は推定の不安定性として現れる。推定の不安定さを抑える方法として、事前分布による制約が考えられる。例えば、パラメータ
が絶対閾値を反映するように事前分布を設定することが考えられる。リストG1では、パラメータ
の範囲を、0以上、提示刺激値の最小値以下としている。パラメータの推定が不安定なときに事前分布による制約を用いる方法は、Wichmann & Hill(2001)により推奨されている。,
リストG1のStanスクリプトを用いるPythonスクリプトをリストG2に示す。リストG1およびリストG2のファイルなどは、圧縮ファイルmefiles.zipにまとめた。ダウンロード展開して自由に利用されたい。
リストG1のファイル(MEc.stan)、リストG2のファイル(ME_cBysn.py)と入力データファイルを同じフォルダ内に置き、そのフォルダをカレントディレクトリとする。CmdStanPyのインストールされた環境で次のコマンドを実行すればよい。
(stan) *****/mefiles$ python ME_cBysn.py
上のコマンドを実行すると、入力データファイル名と出力ファイル名の設定が求められる。次の例では、入力データファイルとしてData.xlsxが、出力ファイルとしてResults.txtが設定されている。
(stan) *****/mefiles$
python ME_cBysn.py
Input Data File (*.xlsx) = Data.xlsx
Output Text File (*.txt) = Results.txt
出力ファイルは適当なテキストファイル名を設定すればよい。
入力データファイルは、Excelファイルとして図G1の形式で用意する。
・
・
・
図G1
1行目は、データの変数名を入力する。データは2行目から入力する。1列目はデータの識別用で、図1では試行が入力されている。この第1列目は、データの識別用であり、数値でも文字でも何でもよい。また、識別用であるが、同じ値を繰り返して用いてもよい。データの分析においては無視される。2列目に提示刺激の物理量、3列目にマグニチュード推定値
を入力する。
出力ファイル名を設定すると、入力データが読み込まれ、リストG1のStanスクリプトのコンパイルが始まる。このコンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが行われ、サンプリングが終了すると、トレース図が表示される(図G2)。
図G2
図G2のフォームを閉じると、パラメータ、
、
の3次元散布図が表示される(図G3)。
図G3
図G3の3D散布図は、図上をマウスでドラッグすると回転することができる(図G4)。
図G4
図G4のWindowを閉じると、図G5のグラフが表示される。
図G5
パラメータa、b、cの事後分布のグラフである。
図G5のWindowを閉じると、図G6のグラフが表示される。
図G6
パラメータ対の散布図である。
図G6のWindowを閉じると、図G7のグラフが表示される。
図G7
左側のグラフでは、データの点と、事後分布の中央値をパラメータの推定値として式(10.2.2)のグラフが描かれている。Gelmanら(2021)は、事後分布の代表値として、平均値より安定した中央値を勧めている。
右側のグラフでは、MCMCサンプリングの結果から500個のサンプル値をランダムに選び、各サンプル値に対して点線でグラフ
が描かれている。第1四分位数Q!と第3四分位数Q3に対応する値が水色とライムの破線で描かれている。
中央の赤線は、各パラメータの中央値を推定値としてグラフ(10.2.2)を描いたものである。
図G7のフォームを閉じると、スクリプトの実行終了である。
リストG1 モデル(10.2.4)のStanスクリプト(MEc.stan)
data {
int N;
array[N] real Sts;
array[N] real MEs;
}
transformed data {
array[N] real LogMEs;
real minSt;
for (i
in 1:N) {
LogMEs[i] = log(MEs[i]);
}
minSt
= Sts[1];
for (i
in 2:N) {
if (minSt > Sts[i]) {
minSt = Sts[i];
}
}
}
parameters {
real<lower = 0.0> a;
real<lower = 0.0, upper =
minSt> c;
real<lower = 0.0> b;
real<lower = 0.0> sgm;
}
model {
a ~ exponential(0.001);
b ~ exponential(0.001);
c ~ uniform(0.0, minSt);
sgm
~ exponential(0.001);
for (i
in 1:N) {
LogMEs[i] ~ normal(log(a * pow(Sts[i] - c, b)), sgm);
// Eq. (10.2.4)
}
}
リストG2 リストG1を利用するPythonスクリプト(ME_cBysn.py)
from cmdstanpy
import CmdStanModel
import numpy
as np
import pandas as pd
import arviz
as az
import matplotlib as mplt
import matplotlib.pyplot
as plt
import seaborn as sb
finnm = input('Input Data File (*.xlsx) =
')
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
foutnm = input('Output Text File (*.txt) =
')
fout = open(foutnm,
'w')
fout.write('Input Data File = {}\n'.format(finnm))
N = len(data)
- 1
Sts = [v for v in data[1:, 1]]
MEs = [v for v in data[1:, 2]]
print('N =', N)
print(len(Sts))
print(len(MEs))
for st, me
in zip(Sts, MEs):
print(st,
me)
Data = {'N':N, 'Sts':Sts, 'MEs':MEs}
model = CmdStanModel(stan_file="MEc.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
print(fit.summary())
fout.write('\n{}'.format(fit.summary()))
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data, figsize=(10,7)) # トレースプロット
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(111,
projection = '3d')
ax.scatter(dframe['a'],
dframe['b'], dframe['c'])
ax.set_xlabel('a', fontsize
= 16)
ax.set_ylabel('b', fontsize
= 16)
ax.set_zlabel('c', fontsize
= 16)
plt.title('Scatter plot of (a,b,c)',
fontsize = 20)
plt.show()
plt.figure(figsize=(10,8))
plt.subplot(2,4,(1,2))
mean_a = dframe['a'].mean()
P025_a, P50_a, P975_a = np.percentile(dframe['a'], [2.5,
50, 97.5])
plt.title('Posterior Distribution of a', fontsize = 20)
sb.kdeplot(dframe['a'],
clip = (0.0, max(dframe['a'])),
label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI =
[{2:.2f}, {3:.2f}]'.
format(mean_a, P50_a, P025_a, P975_a))
plt.xlabel('a', fontsize
= 18)
plt.ylabel('Density', fontsize
= 14)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.subplot(2,4,(3,4))
mean_b = dframe['b'].mean()
P025_b, P50_b, P975_b = np.percentile(dframe['b'], [2.5,
50, 97.5])
plt.title('Posterior Distribution of b', fontsize = 20)
sb.kdeplot(dframe['b'],
label = 'Mean = {0:.2f}\nMed = {1:.2f}\n95%CI =
[{2:.2f}, {3:.2f}]'.
format(mean_b, P50_b, P025_b, P975_b))
plt.xlabel('b', fontsize
= 18)
plt.ylabel('Density', fontsize
= 14)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.subplot(2,4,(6,7))
mean_c = dframe['c'].mean()
P025_c, P50_c, P975_c = np.percentile(dframe['c'], [2.5,
50, 97.5])
plt.title('Posterior Distribution of c', fontsize = 20)
sb.kdeplot(dframe['c'],
label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI =
[{2:.2f}, {3:.2f}]'.
format(mean_c, P50_c, P025_c, P975_c))
plt.xlabel('c', fontsize
= 18)
plt.ylabel('Density', fontsize
= 14)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.show()
plt.figure(figsize=(10,
8))
plt.subplot(2,4,(1,2))
plt.title('Scattergram of (a,b)',
fontsize = 20)
plt.plot(dframe['a'],
dframe['b'], 'bo')
plt.xlabel('a', fontsize
= 18)
plt.ylabel('b', fontsize
= 18)
plt.tight_layout()
plt.subplot(2,4,(3,4))
plt.title('Scattergram of (a,c)',
fontsize = 20)
plt.plot(dframe['a'],
dframe['c'], 'bo')
plt.xlabel('a', fontsize
= 18)
plt.ylabel('c', fontsize
= 18)
plt.tight_layout()
plt.subplot(2,4,(6,7))
plt.title('Scattergram of (b,c)',
fontsize = 20)
plt.plot(dframe['b'],
dframe['c'], 'bo')
plt.xlabel('b', fontsize
= 18)
plt.ylabel('c', fontsize
= 18)
plt.tight_layout()
plt.show()
min_St = min(Sts)
max_St = max(Sts)
St_Values = np.arange(min_St, max_St, (max_St - min_St) / 1000.0)
def PowerLaw(x,
a, b, c):
return a * ((x - c) ** b)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.plot(Sts, MEs, 'bo', alpha = 0.5, label = 'Data')
plt.plot(St_Values,
[PowerLaw(v, P50_a, P50_b, P50_c) for v in St_Values], 'k-',
label = r'$\psi=a(\varphi-c)^b$' + ('\na = {0:.3f}'.format(P50_a)) +
('\nb = {0:.3f}'.format(P50_b)) +
('\nc = {0:.3f}'.format(P50_c)))
plt.xlabel(r'$\varphi$',
fontsize = 18)
plt.ylabel(r'$\psi$', fontsize
= 18)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.subplot(1,2,2)
St_grid = np.linspace(min_St, max_St, 100)
tn_samples = len(dframe)
n_samples = 500
idx_samples = np.random.choice(range(tn_samples), n_samples, replace =
False)
curves_smpl = []
for s in idx_samples:
va
= dframe['a'][s]
vb
= dframe['b'][s]
vc
= dframe['c'][s]
phis = PowerLaw(St_grid, va, vb,
vc)
plt.plot(St_grid, phis, ls = ':', lw = 1,
alpha = 0.5)
curves_smpl.append(phis)
curves_smpl = np.array(curves_smpl)
c_s_q1, c_s_q3 = np.percentile(curves_smpl, [25, 75], axis = 0)
plt.plot(St_grid,
c_s_q1, c = 'aqua', ls = '--', label = 'Q1')
plt.plot(St_grid,
c_s_q3, c = 'lime', ls = '--', label = 'Q3')
plt.plot(St_Values,
[PowerLaw(v, P50_a, P50_b, P50_c) for v in St_Values],
'r-', lw=3, label = r'$\psi=a(\varphi-c)^b$' + ('\na =
{0:.3f}'.format\
(P50_a)) + ('\nb = {0:.3f}'.format(P50_b)) +
('\nc = {0:.3f}'.format(P50_c)))
plt.plot(Sts, MEs,
'o', c = 'b', alpha = 0.5, label = 'Data')
plt.xlabel(r'$\varphi$',
fontsize = 18)
plt.ylabel(r'$\psi$', fontsize
= 18)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.show()
fout.close()
print()
print(foutnm
+ ' was saved.')
print()
確率モデル
のStanスクリプトをリストS1のように用意した。
リストS1を用いるPythonスクリプトは、リストS2のように用意した。使用法は、リストG2と同じである。スクリプトファイルなどは、まとめて圧縮ファイルmefiles.zipとした。ダウンロード展開すれば自由に利用できる。
リストS1のファイル、リストS2のファイル、および入力データファイルを同じフォルダに置いてカレントディレクトリをそのフォルダに移動して、CmdStanPyのインストールされた環境において、次のようなコマンドを実行すればよい。CmdStanPyの準備の簡単な説明は、このウェブサイトで行った。
(stan) *****/mefiles$ python MEBysn.py
上のコマンドを実行すると、入力データファイル名の設定が求められる。
(stan) *****/mefiles$
python MEBysn.py
Input Data File (*.xlsx) = Data.xlsx
Output Text File (*.txt) = Results.txt
入力データファイルは、Excelファイルとして図G1に示す様式で準備する。入力データファイル名を入力すると、次に出力ファイル名の設定が求められる。出力ファイルは、テキストファイル(ファイル拡張子が「.txt」)であれば任意の名前でよい。
出力ファイル名を設定すると、データが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルには多少時間が掛かる。コンパイル後、MCMCサンプリングが実行され、サンプリング、MCMCサンプリングのトレース図(図S1)が表示される。
図S1
図S1のWindowを閉じると、図S2のグラフが表示される。
図S2
パラメータaとbの事後分布と、aとbの散布図である。
図S2のWindowを閉じると、図S3のグラフが表示される。
図S3
左側のグラフは、事後分布の中央値を点推定値として次式モデル
を描いたものである。中央値を点推定値とすることは、Gelmanら(2021)が平均値より安定しているとして勧めている。
右のグラフは、MCMCサンプリングの結果から500個のサンプル値をランダムに選び、各サンプル値に対して点線でグラフ
を描いたものである。第1四分位数Q!と第3四分位数Q3に対応する値が水色とライムの破線で描かれている。中央の赤線は、各パラメータの中央値を推定値としてグラフ
を描いたものである。
図S3のWindowを閉じると、スクリプトの実行終了である。
端末には、以下のようなメッセージが表示されている。
Results.txt was saved.
(stan) *****/mefiles$
出力ファイルResults.txtを開くと、以下のような内容である。
Input Data File = Data.xlsx
Mean
MCSE StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ 28.077800 0.043047 1.348590 ... 1266.14 1898.61 1.00371
a 0.093302 0.001147 0.034230 ... 1313.06 1480.39 1.00503
b 0.615792 0.001877 0.055588 ... 1281.25 1445.08 1.00501
sgm 0.324856 0.000961 0.035652 ... 1559.21 2333.04 1.00312
[4 rows x 11 columns]
リストS1 モデル(10.2.3)のStanスクリプト(ME.stan)
data {
int N;
array[N] real Sts;
array[N] real MEs;
}
parameters {
real<lower = 0.0> a;
real<lower = 0.0> b;
real<lower = 0.0> sgm;
}
model {
a ~ exponential(0.001);
b ~ exponential(0.001);
sgm
~ exponential(0.001);
for (i
in 1:N) {
log(MEs[i]) ~ normal(log(a * pow(Sts[i], b)), sgm); // Eq. (10.2.3)
}
}
リストS2 リストS1のStanスクリプトを用いるPythonスクリプト(MEBysn.py)
from cmdstanpy
import CmdStanModel
import numpy
as np
import pandas as pd
import arviz
as az
import matplotlib as mplt
import matplotlib.pyplot
as plt
import seaborn as sb
finnm = input('Input Data File (*.xlsx) =
')
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
foutnm = input('Output Text File (*.txt) =
')
fout = open(foutnm,
'w')
fout.write('Input Data File = {}\n'.format(finnm))
N = len(data)
- 1
Sts = [v for v in data[1:, 1]]
MEs = [v for v in data[1:, 2]]
print('N =', N)
print(len(Sts))
print(len(MEs))
for st, me
in zip(Sts, MEs):
print(st,
me)
Data = {'N':N, 'Sts':Sts, 'MEs':MEs}
model = CmdStanModel(stan_file="ME.stan") # コンパイル
fit = model.sample(data=Data)
# サンうリング
print(fit.summary())
fout.write('\n{}'.format(fit.summary()))
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data, figsize=(10,7)) # トレースプロット
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
plt.figure(figsize=(10,7))
mean_a = dframe['a'].mean()
P025_a, P50_a, P975_a = np.percentile(dframe['a'], [2.5,
50, 97.5])
plt.subplot(2,4,(1,2))
plt.title('Posterior Distribution of a', fontsize = 20)
sb.kdeplot(dframe['a'],
clip = (0.0, max(dframe['a'])),
label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI =
[{2:.2f}, {3:.2f}]'.
format(mean_a, P50_a, P025_a, P975_a))
plt.xlabel('a', fontsize
= 18)
plt.ylabel('Density', fontsize
= 14)
plt.legend(fontsize =
12)
plt.tight_layout()
mean_b = dframe['b'].mean()
P025_b, P50_b, P975_b = np.percentile(dframe['b'], [2.5,
50, 97.5])
plt.subplot(2,4,(3,4))
plt.title('Posterior Distribution of b', fontsize = 20)
sb.kdeplot(dframe['b'],
label = 'Mean = {0:.2f}\nMed = {1:.2f}\n95%CI =
[{2:.2f}, {3:.2f}]'.
format(mean_b, P50_b, P025_b, P975_b))
plt.xlabel('b', fontsize
= 18)
plt.ylabel('Density', fontsize
= 14)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.subplot(2,4,(6,7))
plt.title('Scattergram of (a,b)',
fontsize = 20)
plt.plot(dframe['a'],
dframe['b'], 'bo')
plt.xlabel('a', fontsize
= 18)
plt.ylabel('b', fontsize
= 18)
plt.tight_layout()
plt.show()
min_St = min(Sts)
max_St = max(Sts)
St_Values = np.arange(min_St, max_St, (max_St - min_St) / 1000.0)
def PowerLaw(x,
a, b):
return a * (x ** b)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(Sts, MEs, 'bo', alpha = 0.5, label = 'Data')
plt.plot(St_Values,
[PowerLaw(v, P50_a, P50_b) for v in St_Values], 'k-',
label = r'$\psi=a\varphi^b$' + ('\na = {0:.3f}'.format(P50_a)) +
('\nb = {0:.3f}'.format(P50_b)))
plt.xlabel(r'$\varphi$',
fontsize = 18)
plt.ylabel(r'$\psi$', fontsize
= 18)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.subplot(1,2,2)
St_grid = np.linspace(min_St, max_St, 100)
tn_samples = len(dframe)
n_samples = 500
idx_samples = np.random.choice(range(tn_samples), n_samples, replace =
False)
curves_smpl = []
for s in idx_samples:
va
= dframe['a'][s]
vb
= dframe['b'][s]
phis = PowerLaw(St_grid, va, vb)
plt.plot(St_grid, phis, ls = ':', lw = 1,
alpha = 0.5)
curves_smpl.append(phis)
curves_smpl = np.array(curves_smpl)
c_s_q1, c_s_q3 = np.percentile(curves_smpl, [25, 75], axis = 0)
plt.plot(St_grid,
c_s_q1, c = 'aqua', ls = '--', label = 'Q1')
plt.plot(St_grid,
c_s_q3, c = 'lime', ls = '--', label = 'Q3')
plt.plot(St_Values,
[PowerLaw(v, P50_a, P50_b) for v in St_Values],
'r-', lw = 3, label = r'$\psi=a\varphi^b$' + ('\na =
{0:.3f}'.format\
(P50_a)) + ('\nb = {0:.3f}'.format(P50_b)))
plt.plot(Sts, MEs,
'o', c = 'b', alpha = 0.5, label = 'Data')
plt.xlabel(r'$\varphi$',
fontsize = 18)
plt.ylabel(r'$\psi$', fontsize
= 18)
plt.legend(fontsize =
12)
plt.tight_layout()
plt.show()
fout.close()
print()
print(foutnm
+ ' was saved.')
print()
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
岡本安晴(2014)心理学データ分析と測定.勁草書房
岡本安晴(2019)いまさら聞けないPythonでデータ分析.丸善出版
岡本安晴(2025)感覚・知覚測定法.和氣典二・重野純・村上郁也(編)感覚・知覚心理学ハンドブック 第三版(第2章)、誠信書房
Stevens, S. S. (1975). Psychophysics: Introduction to its Perceptual, Neural, and Social Prospects. Wiley.
Wichmann, F. A. & Hill, N. J. (2001). The psychometric function: I. Fitting, sampling, and goodness of fit. Perception & Psychophysics, 63, 1293-1313.