第10章 感覚の強さを測る
マグネチュード推定法
べき法則
あるいは、その一般化された式
による分析プログラム例を用意した。式番号は＜岡本安晴「心理学データ分析と測定」勁草書房＞におけるものである。
図１に示されているデータを分析する。上掲書の表10.1.1のものである。
・
・
・
図１ マグネチュード推定法データ
Excelにおいて、1行目に変数名を書く。１列目にはデータの識別情報を設定する。通し番号でも文字列でも何でもよい。第２列目に提示刺激値、第３列目に刺激値に対するマグネチュード推定値を書く。図１のExcelデータファイルは次のスクリプトで読み込み、Stで刺激値、MEでマグネチュード推定値のデータを表すようにする。
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
VarLabels = data.columns[1:].values;
X = data.values
St = X[:, 1]
ME = X[:, 2]
読み込んだデータをモデル（10.2.4）に基づいて分析する。
Stanスクリプトは以下のようになる。
data {
int N;
real St[N];
real ME[N];
real min_c;
}
transformed data {
real logME[N];
for (i in 1:N)
logME[i] = log(ME[i]);
}
parameters {
real<lower = 0.0> a;
real<lower = 0.0> b;
real<lower = 0.0> pc;
real<lower = 0.0> sgm;
}
transformed parameters {
real c;
c = min_c - pc;
}
model {
pc ~ exponential(0.001);
sgm ~ exponential(0.001);
a ~ exponential(0.001);
b ~ exponential(0.001);
for (i in 1:N) {
logME[i] ~ normal(log(a) + b * log(St[i] - c), sgm);
}
}
上のStanスクリプトを用いる分析プログラムをリスト１としてこのウェブサイトの後半に示した。リスト１のファイルなどは圧縮ファイルfiles10.zipとしてまとめたので、ダウンロードして自由に利用できる。
リスト１のプログラムを実行すると、入力データファイル名が聞かれる（図２）。
図２
図２の場合は、図１のデータファイル名Data.xlsxを設定している。入力データファイル名を設定すると、データが読み込まれて計算が始まる。Stanによるサンプリングが終了すると、３つのパラメータの点推定値がそれぞれの事後分布の中央値として与えられて、式（10.2.2）のグラフがデータの散布図とともに表示される（図３）。
図３
図３のウィンドウを閉じると、図４のグラフが表示される。
図４
図４のウィンドウを閉じると、プログラムの実行終了となる。
実行終了後、端末には以下の出力が表示されている。
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 0.09 3.1e-3 0.09 2.9e-3 0.02 0.06 0.12 0.35 923 1.0
b 0.69 5.7e-3 0.16 0.42 0.57 0.67 0.79 1.05 796 1.01
pc 253.9 4.75 135.22 69.13 156.69 227.34 326.67 582.01 811 1.0
sgm 0.32 8.9e-4 0.03 0.27 0.3 0.32 0.34 0.4 1451 1.0
c -73.9 4.75 135.22 -402.0 -146.6
-47.34 23.31 110.87 811 1.0
lp__ 32.79 0.05 1.43 29.09 32.1 33.11 33.81 34.55 935 1.0
Samples were drawn using NUTS at Sun
May 24 15:34:09 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).
a(mean) = 0.0867
95%CI = [0.00289, 0.346]
a(Q1) = 0.0243, a(Med) = 0.0561, a(Q3) = 0.116
b(mean) = 0.69
95%CI = [0.422, 1.05]
b(Q1) = 0.574, b(Med) = 0.674, b(Q3) = 0.787
c(mean) = -73.9
95%CI = [-4.01e+02,
1.11e+02]
c(Q1) = -1.47e+02, c(Med) = -47.3, c(Q3) = 23.3
パラメータｃの事後分布の第１四分位数Q1と第３四分位数Q3の間に０があるので、ｃの値を０としたモデル（10.2.1）で分析してみる。このときのStanスクリプトは以下のようになる。
data {
int N;
real St[N];
real ME[N];
}
transformed data {
real logME[N];
for (i in 1:N)
logME[i] = log(ME[i]);
}
parameters {
real<lower = 0.0> a;
real<lower = 0.0> b;
real<lower = 0.0> sgm;
}
model {
sgm ~ exponential(0.001);
a ~ exponential(0.001);
b ~ exponential(0.001);
for (i in 1:N) {
logME[i] ~ normal(log(a) + b * log(St[i]), sgm);
}
}
このStanスクリプトを用いるPythonスクリプト例をリスト２としてウェブサイトの後半に示した。リスト２のファイルなどは、圧縮ファイルfiles10.zipとしてまとめたので、ダウンロードして自由に利用できる。
リスト２のプログラムを実行すると、入力データファイル名の設定が求められる（図５）。
図５
入力データファイルは図１の形式で用意する。入力データファイル名を設定すると計算が始まる。Stanによるサンプリングが終了すると、図６のグラフが表示される。
図６
図６のウィンドウを閉じると、図７のグラフが表示される。
図７
図７のウィンドウを閉じるとプログラムの実行終了となる。実行終了後の端末には、以下のような出力が表示されている。
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 0.09 1.0e-3 0.03 0.04 0.07 0.09 0.11 0.18 1123 1.0
b 0.62 1.7e-3 0.06 0.51 0.58 0.62 0.65 0.73 1079 1.0
sgm 0.32 9.0e-4 0.03 0.27 0.3 0.32 0.35 0.4 1467 1.0
lp__ 28.09 0.04 1.28 24.68 27.54 28.41 29.0 29.54 1215 1.0
Samples were drawn using NUTS at Sun
May 24 16:23:09 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).
a(mean) = 0.0932
95%CI = [0.0433, 0.176]
a(Q1) = 0.0686, a(Med) = 0.0867, a(Q3) = 0.113
b(mean) = 0.616
95%CI = [0.506, 0.727]
b(Q1) = 0.576, b(Med) = 0.617, b(Q3) = 0.655
本ウェブサイトのPythonの実行例は、Windows上でのAnacondaによる仮想環境において行った。簡単なAnacondaの説明のウェブサイトを用意している。
PythonにおけるStan（PyStan）の説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.optimize as so
import pystan
import seaborn as sb
import pickle
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values
print('X =\n', X)
St = X[:, 1]
ME = X[:, 2]
i = 1
for v in zip(St, ME):
print(i, v)
i += 1
N = len(St)
print('N = ', N)
min_St = St.min()
print('min_St = ', min_St)
Data = {'N': N, 'St': St, 'ME': ME,
'min_c': min_St}
sm = pystan.StanModel(file =
'MEabc.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
a_p025, a_Q1, a_med, a_Q3, a_p975 =
np.percentile(fit['a'], [2.5, 25, 50, 75, 97.5])
b_p025, b_Q1, b_med, b_Q3, b_p975 =
np.percentile(fit['b'], [2.5, 25, 50, 75, 97.5])
c_p025, c_Q1, c_med, c_Q3, c_p975 =
np.percentile(fit['c'], [2.5, 25, 50, 75, 97.5])
print('\na(mean) =
{0:.3}'.format(fit['a'].mean()))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(a_p025, a_p975))
print(' a(Q1) = {0:.3}, a(Med) = {1:.3}, a(Q3) = {2:.3}'.format(
a_Q1, a_med, a_Q3))
print('\nb(mean) =
{0:.3}'.format(fit['b'].mean()))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(b_p025, b_p975))
print(' b(Q1) = {0:.3}, b(Med) = {1:.3}, b(Q3) = {2:.3}'.format(
b_Q1, b_med, b_Q3))
print('\nc(mean) =
{0:.3}'.format(fit['c'].mean()))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(c_p025, c_p975))
print(' c(Q1) = {0:.3}, c(Med) = {1:.3}, c(Q3) = {2:.3}'.format(
c_Q1, c_med, c_Q3))
St_vs = np.arange(St.min(), St.max(),
(St.max() - St.min()) / 1000)
ME_hat = a_med * ((St_vs - c_med) **
b_med)
plt.plot(St, ME, 'bo')
plt.plot(St_vs, ME_hat, 'g', label =
'$\psi = a * (\phi - c)^b$')
plt.xlabel('St', fontsize = 14)
plt.ylabel('ME', fontsize = 14)
plt.title('a(med) = {0:.3} b(med) = {1:.3} c = {2:.3}'.format(
a_med, b_med, c_med), fontsize = 14)
plt.legend(fontsize = 14)
plt.show()
plt.figure(figsize = (12, 6))
plt.subplot(2, 3, 1)
sb.kdeplot(fit['a'])
plt.title('Post. Distri. of a',
fontsize = 14)
plt.subplot(2, 3, 2)
sb.kdeplot(fit['b'])
plt.title('Post. Distri. of b',
fontsize = 14)
plt.subplot(2, 3, 3)
sb.kdeplot(fit['c'])
plt.title('Post. Distri. of c', fontsize
= 14)
plt.subplot(2, 3, 4)
plt.scatter(fit['a'], fit['b'], alpha
= 0.3)
plt.xlabel('a', fontsize = 12)
plt.ylabel('b', fontsize = 12)
plt.title('Scatter plot (a, b)',
fontsize = 14)
plt.subplot(2, 3, 5)
plt.scatter(fit['b'], fit['c'], alpha
= 0.3)
plt.xlabel('b', fontsize = 12)
plt.ylabel('c', fontsize = 12)
plt.title('Scatter plot (b, c)',
fontsize = 14)
plt.subplot(2, 3, 6)
plt.scatter(fit['c'], fit['a'], alpha
= 0.3)
plt.xlabel('c', fontsize = 12)
plt.ylabel('a', fontsize = 12)
plt.title('Scatter plot (c, a)',
fontsize = 14)
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.optimize as so
import pystan
import seaborn as sb
import pickle
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values
print('X =\n', X)
St = X[:, 1]
ME = X[:, 2]
i = 1
for v in zip(St, ME):
print(i, v)
i += 1
N = len(St)
print('N = ', N)
Data = {'N': N, 'St': St, 'ME':
ME} #, 'min_c':
min_St}
sm = pystan.StanModel(file =
'MEab.stan')
fit = sm.sampling(data = Data, n_jobs =
1)
print(fit)
a_p025, a_Q1, a_med, a_Q3, a_p975 =
np.percentile(fit['a'], [2.5, 25, 50, 75, 97.5])
b_p025, b_Q1, b_med, b_Q3, b_p975 =
np.percentile(fit['b'], [2.5, 25, 50, 75, 97.5])
print('\na(mean) =
{0:.3}'.format(fit['a'].mean()))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(a_p025, a_p975))
print(' a(Q1) = {0:.3}, a(Med) = {1:.3}, a(Q3) = {2:.3}'.format(
a_Q1, a_med, a_Q3))
print('\nb(mean) =
{0:.3}'.format(fit['b'].mean()))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(b_p025, b_p975))
print(' b(Q1) = {0:.3}, b(Med) = {1:.3}, b(Q3) = {2:.3}'.format(
b_Q1, b_med, b_Q3))
St_vs = np.arange(St.min(), St.max(),
(St.max() - St.min()) / 1000)
ME_hat = a_med * (St_vs ** b_med)
plt.plot(St, ME, 'bo')
plt.plot(St_vs, ME_hat, 'g', label =
'$\psi = a * (\phi^b)$')
plt.xlabel('St', fontsize = 14)
plt.ylabel('ME', fontsize = 14)
plt.title('a(med) = {0:.3} b(med) = {1:.3}'.format(a_med,
b_med), fontsize = 14)
plt.legend(fontsize = 14)
plt.show()
plt.figure(figsize = (8, 6))
plt.subplot(2, 2, 1)
sb.kdeplot(fit['a'])
plt.title('Post. Distri. of a',
fontsize = 14)
plt.subplot(2, 2, 2)
sb.kdeplot(fit['b'])
plt.title('Post. Distri. of b',
fontsize = 14)
plt.subplot(2, 2, 3)
plt.scatter(fit['a'], fit['b'], alpha
= 0.3)
plt.xlabel('a', fontsize = 12)
plt.ylabel('b', fontsize = 12)
plt.title('Scatter plot (a, b)',
fontsize = 14)
plt.tight_layout()
plt.show()