マグニチュード推定法
ベイズ分析
（PyStan 3 対応版のウェブサイトも用意した）
プロセティックな感覚（量的な感覚）の場合、物理刺激量
と感覚量
の間の関係を表す式として、べき法則と呼ばれている次式
あるいは、これをより一般化した次式
がある。ただし、式番号は＜岡本安晴（2014）「心理学データ分析と測定」勁草書房＞におけるものである。
心理物理学実験において、第
試行で提示された物理刺激の物理量を
で、その感覚量のマグニチュード推定量（Magnitude Estimation: ME）を
で表すとき、以下のモデルを設定する（Stevens, 1975; p. 169）。
モデル（10.2.1）に対しては、
モデル（10.2.2）に対しては、
上記の確率モデルをデータに適用して、ベイズ分析によりパラメータ値を求める。ベイズ分析については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の題３部、第９章〜第15章、で説明している。
まず、一般化されたモデル（10.2.4）による分析について説明する。その後、モデル（10.2.4）において
と置いた場合のモデル（10.2.3）による分析例を示す。
本ウェブサイトのプログラム例とサンプルデータのファイルは、圧縮ファイルfilesMagEstBysn.zipにまとめている。ダウンロード展開すれば、自由に利用できる。
モデル（10.2.4）をStanスクリプトとしてリスト１のように表した。ベイズ分析のためのライブラリStanのPython用モジュールPyStanのインストールについては、このウェブサイトで説明している。
リスト１ モデル（10.2.4）のStanスクリプト（MEc.stan）
data {
int N;
real Sts[N];
real MEs[N];
}
transformed data {
real LogMEs[N];
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)
}
}
リスト１のStanスクリプト（ファイル名MEc.stan）を利用してマグニチュード推定法データ
を分析するPythonスクリプトは、本ウェブサイトの後ろにリストA.1として示した。
リストA.1のファイルMainMagEst_cBysn.pyを実行すると、入力データファイル名の設定が求められる（図１）。なお、改訂版の実行は、Anaconda/Windows 10におけるPython 3.7の仮想環境においてpystan 2.17.1.0をインストールして行った（2024.07.03）。改訂版（リスト１）では、パラメータ
の範囲を、０以上、最小の刺激値以下と設定した。不安定なパラメータを追加したときは、制約を付けて対応することが心理物理学において推奨されている（Wichmann & Hill, 2001, p. 1295）。
図１
入力データファイルは、Excelファイルとして図２の形式で用意する。
・
・
・
図２
図２のデータは、＜岡本安晴「心理学データ分析と測定」勁草書房＞の表10.1.1のデータを入力したものである。１行目は、データの変数名を入力する。データは２行目から入力する。
１列目はデータの識別用で、図２では試行が入力されている。この第１列目は、データの識別用であり、数値でも文字でも何でもよい。また、識別用であるが、同じ値を繰り返して用いてもよい。データの分析においては無視される。２列目に提示刺激の物理量、３列目にマグニチュード推定値を入力する。
入力データファイル名を設定入力すると、データが読み込まれ、次に出力用テキストファイル名の入力が求められる（図１）。
出力用テキストファイル名として適当なテキストファイル名を設定する。図１では、Results.txtと設定している。
出力テキストファイル名の設定後、Stanスクリプトのコンパイルが始まり、コンパイル・ビルド後、MCMCによる事後分布のサンプリングが始まる。
サンプリングが終了すると、パラメータa,b,cの事後分布の散布図が表示される（図３）。
図３
サンプリングによるシミュレーションで推定された事後分布であるので、プログラムの実行ごとにランダムな変動があることに注意。
図３のグラフは、ドラッグして回転することができる（図４）。
図４
図４を見ると、パラメータ
は３次元空間の点であるが、平らなネクタイ状の帯が曲がった形のような分布であることが分かる。このような分布の場合、パラメータを１つずつサンプリング更新する方法では、サンプリングは分布内を偏りなく行うことは難しいと言われている。これに対して、Stanで用いられているHMCという方法は、図４のような分布の場合も分布内を偏りなくサンプリングできることが期待されるとされている。
グラフは保存ボタンをクリックしてファイルに保存することができる。ウィンドウの右上の閉じるボタンをクリックするとウィンドウは閉じられる。
図５のウィンドウを閉じると、パラメータaの事後分布が表示される（図５）。
図５
図５のウィンドウを閉じると、パラメータｂの事後分布が表示される（図６）。
図６
図６のウィンドウを閉じると、パラメータｃの事後分布が表示される（図７）。
図７
パラメータｃの事後分布の95％CIが［2.21, 128.90］であり、
の事後分布が0の方に偏っているので、
とするモデル（10.2.1）およびモデル（10.2.3）の適用が考えられる。
図７のウィンドウを閉じると、パラメータの対（a,b）の散布図が表示される（図８）。
図８
図８のウィンドウを閉じると、パラメータの対（a,c）の散布図が表示される（図９）。
図９
図９のウィンドウを閉じると、パラメータの対（b,c）の散布図が表示される（図10）。
図10
図10のウィンドウを閉じると、図11のグラフが表示される。
図11
図12では、データの点、およびモデル（10.2.2）においてパラメータa,b,cの推定値として事後分布の中央値を採用したときのグラフが描かれている。
図12のウィンドウを閉じると、出力用テキストファイルが保存されたことを示すメッセージが表示されて、プログラムの実行終了となる（図12）。
図12
出力用テキストファイルは、プログラムの実行終了後、開くことができる。上の出力用テキストファイルを開くと、以下のようになっている。
Input Data File = Data.xlsx
Inference for Stan model:
anon_model_1da374e47b77e23ce5e9fb08a7797397.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 0.18 3.4e-3 0.1 0.06 0.11 0.15 0.22 0.46 927 1.0
c 51.99 1.14 36.06 2.2 21.64 46.18 76.74 128.93 1001 1.01
b 0.54 2.3e-3 0.08 0.38 0.49 0.54 0.59 0.68 1065 1.0
sgm 0.33 9.0e-4 0.04 0.27 0.3 0.32 0.35 0.41 1537 1.0
lp__ 30.94 0.04 1.48 27.34 30.17 31.26 32.04 32.84 1249 1.0
Samples were drawn using NUTS at Wed
Jul 3 18:27:07 2024.
For each parameter, n_eff is a crude
measure of effective sample size,
and Rhat is the potential scale
reduction factor on split chains (at
convergence, Rhat=1).
パラメータｃの事後分布（図７）において、「c=0」の可能性が示されたので、次に「c=0」であるモデル（10.2.1）すなわちモデル（10.2.3）による分析を行う。
確率モデル
のStanスクリプトをリスト２のように用意した。ベイズ分析のライブラリStanのPython用インターフェースであるPyStanのインストールについてはこのウェブサイトで説明している。
リスト２ 確率モデル（10.2.3）のスタンスクリプト（ME.stan）
data {
int N;
real Sts[N];
real MEs[N];
}
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)
}
}
リスト２のStanスクリプト（ME.stan）を利用するPythonスクリプトを、本ウェブサイトの後半にリストA.2として示す。スクリプトファイルとサンプルデータファイルは、圧縮ファイルfilesMagEstBysn.zipにまとめたので、ダウンロード展開すれば自由に利用できる。Pythonのインストールは、このウェブサイトで説明している。
リストA.2のプログラムMainMagEstBysn.pyを実行すると、まず入力データファイル名の設定が求められる（図14）。
図14
入力データファイルは、図２の形式のものである。図1４では、図２のファイルを指定している。
入力データファイル名を設定すると、データが読み込まれ、次に出力用テキストファイル名の設定が求められる（図15）。
図15
出力用テキストファイル名を設定すると、Stanスクリプトのコンパイル・ビルドと実行が始まり、MCMCサンプリングが行われる。MCMCサンプリングの終了後、まずパラメータaの事後分布が表示される（図16）。
図16
図16のウィンドウを閉じると、パラメータｂの事後分布が表示される（図17）。
図17
図17のウィンドウを閉じると、パラメータaとbの対（a,b）の散布図が表示される（図18）。
図18
図18のウィンドウを閉じると、データ点と関数
が表示される（図19）。
図19
関数のパラメータaとｂの値は、それぞれの事後分布の中央値である。
図19のウィンドウを閉じると、出力用テキストファイルが保存されたことが表示されて、プログラムの実行終了となる（図20）。
図20
プログラムの実行終了後、出力ファイルを開くと以下のようになっている。
Input Data File = Data.xlsx
Inference for Stan model:
anon_model_b305c1b80ad3836e8ae6a1d7ad789c19.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 0.09 1.1e-3 0.04 0.04 0.07 0.09 0.11 0.18 1037 1.01
b 0.62 1.6e-3 0.06 0.51 0.58 0.62 0.66 0.73 1326 1.01
sgm 0.32 8.9e-4 0.04 0.27 0.3 0.32 0.35 0.41 1573 1.0
lp__ 28.05 0.04 1.35 24.66 27.42 28.43 29.03 29.55 1151 1.0
Samples were drawn using NUTS at Sat
Oct 17 21:35:18 2020.
For each parameter, n_eff is a crude
measure of effective sample size,
and Rhat is the potential scale
reduction factor on split chains (at
convergence, Rhat=1).
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.
リストA.1 リスト１のStanスクリプトを利用するPythonスクリプト
import numpy as np
import pandas as pd
import pystan
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}
sm = pystan.StanModel(file =
'MEc.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fout.write('\n{}'.format(fit))
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(111, projection =
'3d')
ax.scatter(fit['a'], fit['b'],
fit['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 = fit['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(fit['a'], clip = (0.0,
max(fit['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 = fit['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(fit['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 = fit['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(fit['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.title('Scattergram of (a,b)',
fontsize = 20)
plt.plot(fit['a'], fit['b'], 'bo')
plt.xlabel('a', fontsize = 18)
plt.ylabel('b', fontsize = 18)
plt.show()
plt.title('Scattergram of (a,c)',
fontsize = 20)
plt.plot(fit['a'], fit['c'], 'bo')
plt.xlabel('a', fontsize = 18)
plt.ylabel('c', fontsize = 18)
plt.show()
plt.title('Scattergram of (b,c)',
fontsize = 20)
plt.plot(fit['b'], fit['c'], 'bo')
plt.xlabel('b', fontsize = 18)
plt.ylabel('c', fontsize = 18)
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.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.show()
fout.close()
print()
print(foutnm + ' was saved.')
print()
リストA.2 リスト２のStanスクリプトを利用するPythonスクリプト（MainMagEstBysn.py）
import numpy as np
import pandas as pd
import pystan
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}
sm = pystan.StanModel(file =
'ME.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fout.write('\n{}'.format(fit))
mean_a = fit['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(fit['a'], clip = (0.0,
max(fit['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 = fit['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(fit['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()
plt.title('Scattergram of (a,b)',
fontsize = 20)
plt.plot(fit['a'], fit['b'], 'bo')
plt.xlabel('a', fontsize = 18)
plt.ylabel('b', fontsize = 18)
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.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.show()
fout.close()
print()
print(foutnm + ' was saved.')
print()