マグニチュード推定法
ベイズ分析（PyStan 3 対応）
サンプルスクリプトは、Ubuntu 22/WSL/Windows 10においてPython 3.12.3/Anaconda仮想環境で実行した（2024.07.04）。
pystan 3.10.0、arviz 0.18.0を用いた。Arviz 0.18.0は、次のコマンドでインストールした。
pip install arviz==0.18.0
プロセティックな感覚（量的な感覚）の場合、物理刺激量
と感覚量
の間の関係を表す式として、べき法則と呼ばれている次式
あるいは、これをより一般化した次式
がある。ただし、式番号は＜岡本安晴（2014）「心理学データ分析と測定」勁草書房＞におけるものである。
心理物理学実験において、第
試行で提示された物理刺激の物理量を
で、その感覚量のマグニチュード推定量（Magnitude Estimation: ME）を
で表すとき、以下のモデルを設定する（Stevens, 1975; p. 169）。
モデル（10.2.1）に対しては、
モデル（10.2.2）に対しては、
ここで、
は平均
、分散
の正規分布を表す。
上記の確率モデルをデータに適用して、ベイズ分析によりパラメータ値を求める。ベイズ分析については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の題３部、第９章〜第15章、で説明している。ベイズ分析は、PyStan 3を用いているが、WindowsにおいてPyStan 3を利用する場合は、WindowsにUbuntuをインストールしてUbuntu(WSL: Windows Subsysytem for Linux)で使用すればよい。UbuntuのWindowsへのインストールなどの手続きは簡単である。説明のウェブサイトを用意しているので、参照されたい。
まず、一般化されたモデル（10.2.4）による分析について説明する。その後、モデル（10.2.4）において
と置いた場合のモデル（10.2.3）による分析を示す。
一般化されたモデル（10.2.4）の場合
モデル（10.2.4）に基づいてベイズ分析を行うためのStanスクリプトをリスト１のように用意した。
リスト１ モデル（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> v_a;
real<lower = 0.0> v_b;
real<lower = 0.0, upper =
1.0> v_c;
real<lower = 0.0>
v_sgm;
}
transformed parameters {
real a;
real b;
real c;
real sgm;
a = 10.0 * v_a + 0.001;
b = 10.0 * v_b;
c = v_c * minSt;
sgm = v_sgm + 0.01;
}
model {
v_a ~ exponential(0.1);
v_b ~ exponential(0.1);
v_c ~ beta(1.0, 1.0);
v_sgm ~ exponential(0.1);
for (i in 1:N) {
LogMEs[i] ~ normal(log(a * pow(Sts[i] - c, b)), sgm);
}
}
モデル（10.2.2）では、３つのパラメータの間に強い連関が予想される。パラメータ
が小さくなったときは、パラメータ
の値を大きく、あるいはパラメータ
の値を負の方向に大きくすればモデル（10.2.2）とデータの関係の維持が期待される。このパラメータ間の連関は推定の不安定性として現れる。推定の不安定さを抑える方法として、事前分布による制約が考えられる。例えば、パラメータが絶対閾値を反映するように事前分布を設定することが考えられる。リスト１では、パラメータの範囲を、0以上、提示刺激値の最小値以下としている。パラメータの推定が不安定なときに事前分布による制約を用いる方法は、Wichmann & Hill（2001）により推奨されている。,
リスト１のStanスクリプトを用いるPythonスクリプトをリスト２に示す。リスト１およびリスト２のファイルなどは、圧縮ファイルme_c_PS3files.zipにまとめた。ダウンロード展開して自由に利用されたい。
リスト２のPythonスクリプトをWindowsにインストールされたUbuntu 22上のAnacondaのPython 3.12仮想環境で実行した。
次のコマンド
(py312) ・・・/me_c_PS3files$ python
me_c.py
を実行すると、入力データファイル名の設定が求められる。
(py312) ・・・/me_c_PS3files$ python me_c.py
Input Data File (*.xlsx) =
入力データファイルは、Excelファイルとして図１の形式で用意する。
・
・
・
１行目は、データの変数名を入力する。データは２行目から入力する。１列目はデータの識別用で、図１では試行が入力されている。この第１列目は、データの識別用であり、数値でも文字でも何でもよい。また、識別用であるが、同じ値を繰り返して用いてもよい。データの分析においては無視される。２列目に提示刺激の物理量、３列目にマグニチュード推定値を入力する。
入力データファイル名を設定入力すると（次の例では、図１のファイルData.xlsxが設定されている）、データが読み込まれ、次に出力用テキストファイル名の入力が求められる。
(py312) ・・・/me_c_PS3files$ python me_c.py
Input Data File (*.xlsx) = Data.xlsx
Output Text File (*.txt) =
適当なテキストファイル名を設定する。
Output Text
File (*.txt) = Results.txt
出力ファイル名を設定すると、計算が始まり、MCMCサンプリング後、まずトレース図が表示される（図２）。
図２
図２のフォームを閉じると、パラメータ、、の３D散布図が表示される（図３）。
図３
図３の３D散布図は、図上をマウスでドラッグすると回転することができる。
図４のフォームを閉じると、パラメータの事後分布が表示される（図４）。
図４
図４のフォームを閉じると、パラメータの事後分布が表示される（図５）。
図５
図５のフォームを閉じると、パラメータの事後分布が表示される（図６）。
図６
図６のフォームを閉じると、図７のフォームが表示される。
図７
MCMCサンプリングの結果から500個のサンプル値
をランダムに選び、各サンプル値に対して点線でグラフ
を描いたものである。第１四分位数Q!と第３四分位数Q3に対応する値が青色と緑色の破線で描かれている。
中央の赤線は、各パラメータの中央値を推定値としてグラフ（10.2.2）を描いたものである。
図７のフォームを閉じると、スクリプトの実行終了となり、端末には以下のようなメッセージが出力されている。
Results.txt was saved.
(py312) ・・・/me_c_PS3files$
出力ファイルを開くと、以下のような内容である。
Input Data File = Data.xlsx
Data...
1049
5
1888 10
1049
7
180 4.5
・
・
・
583 3.1
324 2.9
180 2.7
Summary:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
a 0.176 0.104 0.043 0.362 0.002 0.001 4988.0 5462.0 1.0
b 0.540 0.075 0.392 0.677 0.001 0.001 5094.0 5407.0 1.0
c 53.070 35.620 0.009 115.363 0.501 0.354 4734.0 5885.0 1.0
sgm 0.329 0.035 0.265 0.394 0.000 0.000 7590.0 7972.0 1.0
読み込まれたデータと、MCMCサンプリングの統計量が出力されている。
パラメータｃの事後分布（図６）において、「c=0」の可能性が示されたので、次に「c=0」であるモデル（10.2.1）すなわちモデル（10.2.3）による分析を行う。
確率モデル
のStanスクリプトをリスト３のように用意した。
リスト３ 確率モデル（10.2.3）のStanスクリプト（ME.stan）
data {
int N;
array[N] real Sts;
array[N] real MEs;
}
transformed data {
array[N] real LogMEs;
LogMEs = log(MEs);
}
parameters {
real<lower = 0.0> v_a;
real<lower = 0.0> v_b;
real<lower = 0.0>
v_sgm;
}
transformed parameters {
real<lower = 0.0> a;
real<lower = 0.0> b;
real<lower = 0.0> sgm;
a = 10.0 * v_a + 0.001;
b = 10.0 * v_b;
sgm = 10.0 * v_sgm + 0.001;
}
model {
v_a ~ exponential(0.1);
v_b ~ exponential(0.1);
v_sgm ~ exponential(0.1);
for (i in 1:N) {
LogMEs[i] ~ normal(log(a * pow(Sts[i], b)), sgm); // Eq. (10.2.3)
}
}
リスト３を用いるPythonスクリプトは、リスト４のように用意した。使用法は、リスト２と同じである。スクリプトファイルなどは、まとめて圧縮ファイルme_PS3files.zipとした。ダウンロード展開すれば自由に利用できる。Windowsで実行するときは、Anaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境を利用すればよい。
以下のように、pythonコマンドでスクリプトファイルme.pyを実行すると入力データファイル名を聞いてくる。
(py311) ・・・/me_PS3files$ python me.py
Input Data File (*.xlsx) =
入力データファイルはリスト２の場合と同様で、図１の形式で用意する。
入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。適当なテキストファイル名を設定する。
(py311) ・・・/me_PS3files$ python me.py
Input Data File (*.xlsx) = Data.xlsx
Output Text File (*.txt) = Results.txt
出力ファイル名を設定すると、計算が始まる。
MCMCサンプリングが終わると、まずトレースが表示される（図８）。
図８
図８のフォームを閉じると、パラメータとの散布図が表示される（図９）。
図９
図９のフォームを閉じると、パラメータの事後分布が表示される（図10）。
図10
図10のフォームを閉じると、パラメータの事後分布が表示される（図11）。
図11
図11のフォームを閉じると、図12のグラフが表示される。
図12
MCMCサンプリングの結果から500個のサンプル値
をランダムに選び、各サンプル値に対して点線でグラフ
を描いたものである。第１四分位数Q!と第３四分位数Q3に対応する値が青色と緑色の破線で描かれている。
中央の赤線は、各パラメータの中央値を推定値としてグラフ（10.2.1）を描いたものである。
図12のフォームを閉じると、スクリプトの実行終了である。
端末画面には、下記のようなメッセージが表示される。
Results.txt was saved.
(py311) ・・・/me_PS3files$
出力ファイルを開くと、下記のような内容である。
Input Data File = Data.xlsx
Data...
1049 5
1888 10
1049 7
・
・
・
583 3.1
324 2.9
180 2.7
Summary:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
a 0.093 0.034 0.038 0.154 0.000 0.000 5237.0 6214.0 1.0
b 0.616 0.055 0.515 0.720 0.001 0.001 5244.0 6029.0 1.0
sgm 0.324 0.035 0.261 0.390 0.000 0.000 7370.0 6850.0 1.0
岡本安晴（2019）いまさら聞けないPOythonでデータ分析．丸善出版
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.
リスト２ リスト１を利用するPythonスクリプト（me_c.py）
import numpy as np
import pandas as pd
import stan
import matplotlib as mplt
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
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))
fout.write('\nData...\n')
for st, me in zip(Sts, MEs):
print(st, me)
fout.write(f'{st:>10}{me:>10}\n')
Data = {'N':N, 'Sts':Sts, 'MEs':MEs}
with open('MEc.stan', 'r') as fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_samples = 5000)
print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['a', 'b', 'c', 'sgm'])
print('Summary =\n',smry)
fout.write(f'\nSummary:\n{smry}\n')
az.plot_trace(i_data, var_names =
['a', 'b', 'c'])
plt.tight_layout()
plt.show()
d_fm = fit.to_frame()
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(111, projection =
'3d')
ax.scatter(d_fm['a'], d_fm['b'],
d_fm['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()
mean_a = d_fm['a'].mean()
P025_a, P50_a, P975_a =
np.percentile(fit['a'], [2.5, 50, 97.5])
plt.title('Posterior Distribution of
a', fontsize = 20)
sb.kdeplot(d_fm['a'], clip = (0.0,
np.max(d_fm['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.show()
mean_b = d_fm['b'].mean()
P025_b, P50_b, P975_b =
np.percentile(fit['b'], [2.5, 50, 97.5])
plt.title('Posterior Distribution of
b', fontsize = 20)
sb.kdeplot(d_fm['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.show()
mean_c = d_fm['c'].mean()
P025_c, P50_c, P975_c =
np.percentile(fit['c'], [2.5, 50, 97.5])
plt.title('Posterior Distribution of
c', fontsize = 20)
sb.kdeplot(d_fm['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()
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)
St_grid = np.linspace(min_St, max_St,
100)
tn_samples = len(d_fm)
n_samples = 500
idx_samples =
np.random.choice(range(tn_samples), n_samples, replace = False)
curves_smpl = []
for s in idx_samples:
va = d_fm['a'][s]
vb = d_fm['b'][s]
vc = d_fm['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 = 'b', ls
= '--', label = 'Q1')
plt.plot(St_grid, c_s_q3, c = 'g', ls
= '--', label = 'Q3')
plt.plot(St_Values, [PowerLaw(v,
P50_a, P50_b, P50_c) for v in St_Values],
'r-', 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.legend()
plt.show()
print()
print(foutnm + ' was saved.')
print()
リスト４ リスト３を用いるPythonスクリプト（me.py）
import numpy as np
import pandas as pd
import stan
import matplotlib as mplt
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
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))
fout.write('\nData...\n')
for st, me in zip(Sts, MEs):
print(st, me)
fout.write(f'{st} {me}\n')
Data = {'N':N, 'Sts':Sts, 'MEs':MEs}
with open('ME.stan', 'r') as fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_samples = 5000)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['a', 'b', 'sgm'])
print('Summary =\n',smry)
fout.write(f'\nSummary:\n{smry}\n')
az.plot_trace(i_data, var_names =
['a', 'b', 'sgm'])
plt.show()
d_fm = fit.to_frame()
plt.title('Scattergram of (a,b)',
fontsize = 20)
plt.plot(d_fm['a'], d_fm['b'], 'bo')
plt.xlabel('a', fontsize = 18)
plt.ylabel('b', fontsize = 18)
plt.show()
mean_a = np.mean(d_fm['a'])
P025_a, P50_a, P975_a =
np.percentile(d_fm['a'], [2.5, 50, 97.5])
plt.title('Posterior Distribution of
a', fontsize = 20)
sb.kdeplot(d_fm['a'], clip = (0.0,
max(d_fm['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.show()
mean_b = np.mean(d_fm['b'])
P025_b, P50_b, P975_b =
np.percentile(d_fm['b'], [2.5, 50, 97.5])
plt.title('Posterior Distribution of
b', fontsize = 20)
sb.kdeplot(d_fm['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.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)
St_grid = np.linspace(min_St, max_St,
100)
tn_samples = len(d_fm)
n_samples = 500
idx_samples =
np.random.choice(range(tn_samples), n_samples, replace = False)
curves_smpl = []
for s in idx_samples:
va = d_fm['a'][s]
vb = d_fm['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 = 'b', ls
= '--', label = 'Q1')
plt.plot(St_grid, c_s_q3, c = 'g', ls
= '--', label = 'Q3')
plt.plot(St_Values, [PowerLaw(v, P50_a,
P50_b) for v in St_Values],
'r-', 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.legend()
plt.show()
fout.close()
print()
print(foutnm + ' was saved.')
print()