ミュラーリヤー錯視
実験（調整法）とデータ分析のPythonプログラム例
心理学実験とデータ分析をPythonプログラミングで試みた。心理学実験として「ミュラーリヤー錯視」における調整法を取り上げ、プログラミングはライブラリPyGameを利用した。データの分析は、最小2乗法をライブラリSciPyを利用して行い、ベイズ分析はPyStanを利用した。以下に、実験、最小2乗法、ベイズ分析の順に説明する。
ミュラーリヤー実験におけるキー操作よる長さの調節などは、キーイベントなどののイベント処理として行った。Pythonによるデータ分析については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
本サイトのプログラムファイルなどは、圧縮ファイルmuler_lyer_files.zipにまとめた。ダウンロード展開すれば、自由に利用できる。
ミュラーリヤー錯視の提示図形として図E1に示すものを使用した。
図E1
実験のためのプログラムは、リストE1に示すものである。観察者は、図E1のような刺激図形が提示されたとき、左右の矢軸が同じ長さに見えるように調整することが求められる（調整法）。矢軸はキーボードの右矢印→キーを押すと長くなり左矢印キー←を押すと短くなる。左右の矢軸が同じ長さに見えたところで「OK」ボタンをクリックする。「OK」ボタンをクリックすると、少しの時間の後、次の刺激図形が提示されるので、左右が同じ長さに見るように調整して「OK」ボタンをクリックする。この試行を20試行繰り返して実験終了である。
キー操作などの処理は、イベントのチェックとそれに対応する処理という方法で、PyGameを利用して以下のスクリプトで行っている。
for
e in pygame.event.get():
if e.type == pygame.QUIT: # Windowsのcloseボタンがクリックされたか
pygame.quit()
sys.exit()
elif e.type == pygame.KEYDOWN:
if e.key == pygame.K_RIGHT:
# 右矢印キーが押下されたか
ck_go_L = True
if e.key == pygame.K_LEFT: # 左矢印キーが押下されたか
ck_go_S = True
elif e.type == pygame.KEYUP: # キーが離されたか
ck_go_L = False
ck_go_S = False
elif e.type == pygame.MOUSEBUTTONDOWN:
ms_point = pygame.mouse.get_pos()
if ok_button.rect.collidepoint(ms_point): # ボタンがクリックされたか
screen.fill((255,255,255))
# 画面を白一色に消去
pygame.display.flip()
ck_loop = False
break
プログラムでは、ライブラリPyGameの他、データをExcelファイルとして保存するためにライブラリPandasを用いている。PandasによりExcelファイルに出力するときに、ライブライopenpyxlも必要である。すなわち、
pip install
pygame pandas openpyxl
とインストールしておけばよい。
実験が終わると、データの出力ファイル名を聞いてくる。
Output File
Name(*.xlsx) =
適当なExcelファイル名を設定すればよい。
Output File Name(*.xlsx) = sample.xlsx
sample.xlsx was saved.
ファイル名を設定すると、ファイルに保存された旨、表示されてプログラムの終了である。
実験終了後、リストE2のプログラムでデータを見ることができる。リストE2のプログラムの実行には、ライブラリPandasとMatplotlibが必要である。
コマンド
pip install
matplotlib
を実行して、ライブラリMatplotlibを追加インストールしておけばよい。
リストE2のプログラムを実行すると、入力ファイル名を聞いてくるので、表示したいデータファイル名を設定する。次の例では、上の実験で保存したファイルsample.xlsxが設定されている。
Input
File(*.xlsx) = sample.xlsx
ファイル名を設定すると、データが読み込まれて、まずヒストグラムが表示される（図E2）。
図E2
試行の始まりにおいて、右の矢軸が長い刺激図形が提示されるのがU条件で、下図の刺激提示画面が表示され、
右の矢軸が短い刺激図形が提示されるのがL条件で、下図の刺激提示画面が表示される。
観察者は提示された刺激画面において、左右の矢軸の長さが同じに見えるように右側の矢軸の長さを調整することが求められる。
図E２において、U条件およびL条件の平均値と標準偏差がグラフ上部２行目に表示されている。表示１行目に示されている平均値は、U条件の平均値とL条件の平均値の２つの平均値の平均値である。バイアスはこの平均値とU条件における平均値との差である。
調整法において、どの刺激値から始めるのかという要因がPSEの設定に強く影響していることが判る。図E2において、赤の半円は左側の矢軸の物理的長さのピクセル値である。左右の矢軸が物理的に同じ長さであるように右側の矢軸の長さが設定されたとすれば、この長さになる。
図E2のウィンドウを閉じると、図E3の折れ線グラフが表示される。
図E3
奇数試行がL条件、偶数試行がU条件である。U条件でのPSEがＬ条件より大きいということが安定した傾向であることが読み取れる。左の矢軸の長さ（ピクセル値）が、赤色の鎖線で示されている。
リストE１ ミュラーリヤー錯視実験Pythonプログラム例（MullerLyer.py）
import pygame
import sys
import pandas as pd
pygame.init()
screen = pygame.display.set_mode((1200, 800))
pygame.display.set_caption('Muller Lyer')
screen.fill((255,255,255))
pygame.display.flip()
xL = 600 - 350
xR = 600 + 350
L = 80
class Button:
def __init__(self):
self.screen_rect = screen.get_rect()
self.width, self.height = 200, 70
self.button_color = (0,0,255)
self.text_color = (255,255,255)
self.font = pygame.font.SysFont(None, 58)
self.rect = pygame.Rect(600-(self.width//2), 800 - (self.height * 2),
self.width, self.height)
self.msg_image = self.font.render('OK', True, self.text_color,
self.button_color)
self.msg_image_rect = self.msg_image.get_rect()
self.msg_image_rect.center = self.rect.center
def draw_button(self):
screen.fill(self.button_color,
self.rect)
screen.blit(self.msg_image, self.msg_image_rect)
def disp_pair(dL):
screen.fill((255,255,255))
#
#
左右の矢軸の描画
#
pygame.draw.line(screen,
(0,0,0), (xL,400), (xR+dL,400), width = 5)
#
#
左の矢羽根の描画
#
pygame.draw.line(screen,
(0,0,0), (xL+L,400-L), (xL,400), width = 5)
pygame.draw.line(screen,
(0,0,0), (xL+L,400+L), (xL, 400), width = 5)
#
#
中央の矢羽根の描画
#
pygame.draw.line(screen, (0,0,0),
(600-L, 400+L), (600, 400), width = 5)
pygame.draw.line(screen,
(0,0,0), (600-L, 400-L), (600, 400), width = 5)
#
#
右の矢羽根の描画
#
pygame.draw.line(screen,
(0,0,0), (xR+L+dL, 400-L), (xR+dL, 400), width = 5)
pygame.draw.line(screen,
(0,0,0), (xR+L+dL, 400+L), (xR+dL, 400), width = 5)
#
#
ボタンの設置
#
global ok_button
ok_button = Button()
ok_button.draw_button()
#
#
ディスプレイ画面への表示
#
pygame.display.flip()
pygame.time.wait(50)
# 50msecのブレーキ
def disp_longer():
global dL
if xR + L + dL + 1 <
1200:
dL
+= 1
# 右の矢軸の伸長
disp_pair(dL)
def disp_shorter():
global dL
if xR + dL - 1 > 600:
dL
-= 1
# 右の矢軸の収縮
disp_pair(dL)
n_trials = 20 # 総試行数
trials = [] # 試行数
conds = [] # 提示条件
resps = [] # 反応
for i_trial in range(n_trials):
trials.append(i_trial+1)
if i_trial % 2 == 1:
dL =
100
# 最初に十分長い矢軸を提示する
conds.append('U')
else:
dL =
-200
# 最初に十分短い矢軸を提示する
conds.append('L')
disp_pair(dL)
ck_go_L = False
# 矢軸を延ばすかどうかのフラグ
ck_go_S = False
# 矢軸を縮めるかどうかのフラグ
ck_loop = True
# While文のループを継続するかどうかのフラグ
while True:
#
# イベントのチェック
#
for
e in pygame.event.get():
if e.type == pygame.QUIT: # Windowsのcloseボタンがクリックされたか
pygame.quit()
sys.exit()
elif e.type == pygame.KEYDOWN:
if e.key == pygame.K_RIGHT:
# 右矢印キーが押下されたか
ck_go_L = True
if e.key == pygame.K_LEFT: # 左矢印キーが押下されたか
ck_go_S = True
elif e.type == pygame.KEYUP: # キーが離されたか
ck_go_L = False
ck_go_S = False
elif e.type == pygame.MOUSEBUTTONDOWN:
ms_point = pygame.mouse.get_pos()
if ok_button.rect.collidepoint(ms_point): # ボタンがクリックされたか
screen.fill((255,255,255))
# 画面を白一色に消去
pygame.display.flip()
ck_loop = False
break
if not
ck_loop:
break
# 当該試行の終了
if
ck_go_L:
disp_longer() # 右矢軸を伸ばす
if
ck_go_S:
disp_shorter()
# 右矢軸を縮める
resps.append(350 + dL) # 右矢軸の長さの記録
pygame.time.wait(1000) # １秒待ち
pygame.quit()
print('trials =\n', trials)
print('conds =\n', conds)
print('resps =\n', resps)
data = pd.DataFrame({'trial':trials,
'cond':conds, 'resp':resps})
fl_nm = input('Output File
Name(*.xlsx) = ')
with pd.ExcelWriter(fl_nm) as
ex_writer:
data.to_excel(ex_writer,
index=False)
# データをExcelファイルに保存
print(fl_nm, 'was saved.')
リストE2 データ表示用プログラム例（ShowData.py）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
fl_nm = input('Input File(*.xlsx) = ')
with pd.ExcelFile(fl_nm) as
f_data:
# ファイル名の設定
data_xlsx =
pd.read_excel(f_data)
# DataFrame型への読み込み
print(data_xlsx)
# 書き出し
trials = data_xlsx['trial']
conds = data_xlsx['cond']
resps = data_xlsx['resp']
c_idx = conds == 'U'
#
# U条件のデータの取り出し
#
trials_U = trials[c_idx]
conds_U = conds[c_idx]
resps_U = resps[c_idx]
U_mean = np.mean(resps_U)
U_sd = np.std(resps_U)
#
# L条件、すなわち非U条件の取り出し
#
c_idx_L = [not v for v in c_idx]
trials_L = trials[c_idx_L]
conds_L = conds[c_idx_L]
resps_L = resps[c_idx_L]
L_mean = np.mean(resps_L)
L_sd = np.std(resps_L)
g_mean = (U_mean + L_mean) / 2
bias = U_mean - g_mean
#
# ヒストグラム
#
plt.hist(resps_U, label = 'U')
plt.hist(resps_L, label = 'L')
plt.plot([350], [0], c='r', marker =
'o', markersize=20, label = '350')
plt.title(f'Mean = {g_mean:.2f}, Bias = {bias:.2f}\n' +
f'L: mean={L_mean:.2f}, sd={L_sd:.2f}' +
f' U: mean={U_mean:.2f}, sd={U_sd:.2f}')
plt.legend()
plt.show()
#
# 折れ線グラフ
#
plt.plot(trials_U, resps_U, label =
'U')
plt.plot(trials_L, resps_L, label =
'L')
plt.xticks(range(1,21))
plt.plot([1,20], [350,350], c='r',
ls='--', label='350')
plt.legend()
plt.show()
上の実験で得られたデータを最小二乗法を用いて以下のように分析した。
調整法において２つの条件、最初の比較刺激値を標準刺激値より大きい（長い）条件（U条件）と小さい（短い）条件（L条件）、で主観的等価点（PSE）が異なる結果が得られたとき、真の主観的等価点が最初の提示刺激値の値によりバイアスを受けたと考える。すなわち、U条件でのPSEは平均
、分散
の正規分布に従い、L条件でのPSEは平均
、分散の正規分布に従うとする。ここで、
はバイアスが無いときの平均値である。
パラメータ推定を行うとき、データ値の分布を直接対象とするとデータ値の変動の影響を大きく受けるので、累積比率を対象とする。最小二乗法によりパタメータ値を求めるプログラムをリストL1のように作成した。
リストL1では、最小二乗法を適用するための目的関数を、次のクラス型ss_cumのメソッドを用意して、目的関数opt_fを設定している。
class ss_cum:
def __init__(self,
resps_U_s, cum_U, resps_L_s, cum_L):
self.resps_U_s = resps_U_s
self.cum_U = cum_U
self.resps_L_s = resps_L_s
self.cum_L = cum_L
def f_ss(self, m, s, b):
v_ss
= 0.0
for
i in range(len(self.resps_U_s)):
v_ss += (self.cum_U[i] - ss.norm.cdf(resps_U_s[i],
m + b, s)) ** 2
for
i in range(len(self.resps_L_s)):
v_ss += (self.cum_L[i] - ss.norm.cdf(resps_L_s[i],
m - b, s)) ** 2
return v_ss
上のクラス型ss_cumを用いて次の目的関数opt_fが設定されている。
func_ss = ss_cum(resps_U_s, cum_U,
resps_L_s, cum_L)
#
# 最小2乗法のための目的関数
#
def opt_f(x):
m = x[0]
s = x[1]**2
b = x[2]
return func_ss.f_ss(m, s, b)
リストL1を実行すると、入力データファイル名を聞いてくるので、上の実験で保存したデータファイル名を設定する。
Input
File(*.xlsx) = sample.xlsx
入力データファイル名を設定すると、データが読み込まれて、まずヒストグラムが表示される（図L1）。
図L1
図L1のフォームを閉じると、次に累積比率グラフが表示される（図L2）。
図L2
リストL1のプログラムでは、この累積比率グラフに対して求める累積分布関数の2乗誤差が最小になるようにパラメータ値が求められる。
図L2のウィンドウを閉じると計算が始まり、計算結果のパラメータによる累積比率のグラフが表示される（図L3）。
図L3
グラフの上部に、最小2乗法で求められたパラメータ推定値が表示されている。３本の累積比率グラフの中央のグラフは、バイアスbを０と置いときのグラフである。PSEの平均値が
、標準偏差が
、提示条件によるバイアスが
という結果である。
図L3のウィンドウを閉じると、図L4のグラフが表示される。
図L4
図L4は、図L3の累積比率のグラフを、密度関数のグラフで表したものである。
リストL1 最小2乗法によるパラメータ推定プログラム（lsm.py）
import pandas as pd
import numpy as np
import scipy.stats as ss
import scipy.optimize as so
import matplotlib.pyplot as plt
fl_nm = input('Input File(*.xlsx) = ')
with pd.ExcelFile(fl_nm) as
f_data:
# ファイル名の設定
data_xlsx =
pd.read_excel(f_data)
# DataFrame型への読み込み
print(data_xlsx)
# 書き出し
trials = data_xlsx['trial']
conds = data_xlsx['cond']
resps = data_xlsx['resp']
#
# U条件のデータの取り出し
#
idx_U = conds == 'U'
trials_U = trials[idx_U]
conds_U = conds[idx_U]
resps_U = resps[idx_U]
#
# L条件のデータの取り出し
#
idx_L = conds == 'L'
trials_L = trials[idx_L]
conds_L = conds[idx_L]
resps_L = resps[idx_L]
#
# ヒストグラム
#
plt.hist(resps_U, label = 'U')
plt.hist(resps_L, label = 'L')
plt.plot([350], [0], c='c', marker =
'o', markersize=20, label = '350')
plt.legend()
plt.show()
#
# 累積比率:U条件
#
resps_U_s = np.sort(resps_U)
cum_U = []
for i in range(len(resps_U_s)):
cum_U.append((i+1)/len(resps_U_s))
#
# 累積比率：L条件
#
resps_L_s = np.sort(resps_L)
cum_L = []
for i in range(len(resps_L_s)):
cum_L.append((i+1)/len(resps_L_s))
plt.plot(resps_U_s, cum_U, label =
'U')
plt.plot(resps_L_s, cum_L, label = 'L')
plt.title('Cumulative graph')
plt.legend()
plt.show()
#
# 2乗誤差の和を求めるクラス型
#
class ss_cum:
def __init__(self,
resps_U_s, cum_U, resps_L_s, cum_L):
self.resps_U_s = resps_U_s
self.cum_U = cum_U
self.resps_L_s = resps_L_s
self.cum_L = cum_L
def f_ss(self, m, s, b):
v_ss
= 0.0
for
i in range(len(self.resps_U_s)):
v_ss += (self.cum_U[i] - ss.norm.cdf(resps_U_s[i],
m + b, s)) ** 2
for
i in range(len(self.resps_L_s)):
v_ss += (self.cum_L[i] - ss.norm.cdf(resps_L_s[i],
m - b, s)) ** 2
return
v_ss
func_ss = ss_cum(resps_U_s, cum_U,
resps_L_s, cum_L)
#
# 最小2乗法のための目的関数
#
def opt_f(x):
m = x[0]
s = x[1]**2
b = x[2]
return func_ss.f_ss(m, s, b)
#
# 初期値
#
x0 = np.empty(3)
x0[0] = np.median(resps)
x0[1] = np.std(resps) ** 0.5
x0[2] = 0
results = so.minimize(opt_f, x0)
m = results.x[0]
s = results.x[1] ** 2
b = results.x[2]
print('\nm =', m, ' s =', s, ' b =', b)
#
# 求めた結果のグラフ表示：累積分布
#
r_min = np.min(resps)
r_max = np.max(resps)
r_range = np.arange(r_min, r_max,
0.001)
Model_cum_U = ss.norm.cdf(r_range,
m-b, s)
Model_cum_L = ss.norm.cdf(r_range,
m+b, s)
Model_cum_NoBias =
ss.norm.cdf(r_range, m, s)
plt.plot(r_range, Model_cum_U, label =
'U: model')
plt.plot(r_range, Model_cum_L, label =
'L: model')
plt.plot(r_range, Model_cum_NoBias,
label = 'b = 0')
plt.plot(resps_U_s, cum_U, label = 'U:
data')
plt.plot(resps_L_s, cum_L, label = 'L:
data')
plt.title(f'mu = {m:.2f}, sgm =
{s:.2f}, b = {b:.2f}')
plt.legend()
plt.show()
#
# 求めた結果のグラフ表示：比率分布
#
dns_U = ss.norm.pdf(r_range, m+b, s)
dns_L = ss.norm.pdf(r_range, m-b, s)
plt.plot(r_range, dns_U, label = 'U:
model')
plt.plot(r_range, dns_L, label = 'L:
model')
plt.hist(resps_U, density = True,
label = 'U: data')
plt.hist(resps_L, density = True,
label = 'L: data')
plt.plot([350], [0], c='c', marker =
'o', markersize=20, label = '350')
plt.legend()
plt.show()
上の実験で得られたデータに対してベイズ分析を行うために、以下のように確率モデルを設定した。
調整法において２つの条件、最初の比較刺激値を標準刺激値より大きい（長い）条件（U条件）と小さい（短い）条件（L条件）、で主観的等価点（PSE）が異なる結果が得られたとき、真の主観的等価点が最初の提示刺激値の値によりバイアスを受けたと考える。すなわち、U条件でのPSEは平均μ+b、分散の正規分布に従い、L条件でのPSEは平均μ-b、分散の正規分布に従うとする。ここで、μはバイアスが無いときの平均値である。この仮定の下で、ベイズ分析を行うためのStanスクリプトをリストB1のように用意した。
リストB1のスクリプトはPyStan 3対応である。岡本安晴（2019）ではPyStan 2について説明しているが、Pystan 2とPyStan 3の違いについてはこのウェブサイトなどで説明している。
リストB1のStanスクリプトを用いてベイズ分析を行うプログラムをリストB2のように作成した。リストB2のプログラムをWindowsで実行するときは、WSLでWindowsにUbuntuをインストールして、Ubuntu上で実行する。このウェブサイトの説明を参照されたい。
以下の実行例（2024.01.22）は、WSL上のUbuntuでのものである。
リストB2のプログラムを実行すると、入力データファイル（Excelファイル形式）の名前の設定が求められる。
(py311) *****/muler_lyer_files$ python
ba.py
Input File(*.xlsx) = sample.xlsx
上の例では、本ウェブサイトの実験データ例の名前が設定されている。
入力データファイル名を設定してEnterキーを押すと、データが読み込まれて処理が始まり、図B1のウィンドウが表示される。
図B１
グラフは、コード
plt.savefig('Fig_hist.png')
により、ファイルに保存されている。
図B１のウィンドウを閉じると、ベイズ分析のMCMCサンプリングが始まり、サンプリング後、サンプリング結果のグラフが表示される（図B２）。
図B２
図B２のウィンドウを閉じると、パラメータの事後分布のグラフが表示される（図B３）。
図B３
タイトルに点推定値としてMAP推定値（モード）と4分位数が表示されている。Gelmanら（2021）Regression and other stories.では、分布の代表値として平均値は不安定であり、中央値が勧められている（p. 73）が、モードは最良の全体的フィットを与えると説明されている（p. 107）。岡本（2022）は、モード（MAP推定値）がよいと報告している。MAP推定値を算出するときは、MCMCサンプリングのサンプル数を多く取って分布がモード付近で崩れないように気を付ける必要がある。
図B３のウィンドウを閉じると、パラメータｓの事後分布のグラフが表示される（図Ｂ４）。
図B４
図B４のウィンドウを閉じると、パラメータｂの事後分布のグラフが表示される（図Ｂ５）。
図B５
図B５のウィンドウを閉じると、図B6のグラフが表示される。
図B6
ヒストグラムに重ねて、パラメータm、s、bのMAP推定値によるU条件とL条件におけるPSEの分布が描かれている。
図B6のウィンドウを閉じると、プログラムの実行終了となる。
端末画面上には、以下のような出力が表示されている。
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
vm 0.001 0.123 -0.243 0.225 0.001 0.001 31018.0 24625.0 1.0
vs 0.540 0.099 0.373 0.725 0.001 0.000 25101.0 23874.0 1.0
vb 0.879 0.122 0.651 1.112 0.001 0.000 32494.0 25429.0 1.0
m 303.117 2.681 297.828 307.993 0.015 0.011 31018.0 24625.0 1.0
s 11.757 2.156 8.115 15.777 0.014 0.010 25101.0 23874.0 1.0
b 19.096 2.647 14.144 24.179 0.015 0.010 32494.0 25429.0 1.0
parameter m
map = 303.03
quantiles = 301.40, 303.11, 304.86
parameter s
map = 10.87
quantiles = 10.22, 11.46, 12.96
parameter b
map = 19.15
quantiles = 17.39, 19.09, 20.82
リストB1 ベイズ分析のためのStanスクリプト（model_adj.stan）
data {
int mU;
array[mU] real ResU;
int mL;
array[mL] real ResL;
real mean;
real sd;
}
parameters {
real vm;
real<lower=0.0> vs;
real vb;
}
transformed parameters {
//
// Rescaling of the
parameters
//
real m;
real s;
real b;
m = vm * sd + mean;
s = vs * sd + 0.01;
b = vb * sd;
}
model {
vm ~ normal(0.0, 10.0);
vs ~ exponential(0.1);
vb ~ normal(0.0, 10.0);
ResU ~ normal(m+b, s);
ResL ~ normal(m-b, s);
}
リストB2 リストB1のStanスクリプトを用いてベイズ分析を行うプログラム（ba.py）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import stan
import arviz as az
import scipy.stats as ss
fl_nm = input('Input File(*.xlsx) = ')
with pd.ExcelFile(fl_nm) as
f_data:
# ファイル名の設定
data_xlsx =
pd.read_excel(f_data)
# DataFrame型への読み込み
print(data_xlsx) # 書き出し
trials = data_xlsx['trial']
conds = data_xlsx['cond']
resps = data_xlsx['resp']
#
# U条件のデータの取り出し
#
idx_U = conds == 'U'
trials_U = trials[idx_U]
conds_U = conds[idx_U]
resps_U = resps[idx_U]
#
# L条件のデータの取り出し
#
idx_L = conds == 'L'
trials_L = trials[idx_L]
conds_L = conds[idx_L]
resps_L = resps[idx_L]
#
# ヒストグラム
#
plt.hist(resps_U, label = 'U')
plt.hist(resps_L, label = 'L')
plt.plot([350], [0], c='c', marker =
'o', markersize=20, label = '350')
plt.legend()
plt.tight_layout()
plt.savefig('Fig_hist.png')
plt.show()
resps_U = np.array(resps_U) # Stan用に型の再設定
resps_L = np.array(resps_L) # Stan用に型の再設定
v_mean = np.mean(resps)
v_sd = np.std(resps)
Data = {'mU':len(resps_U), 'ResU':resps_U,
'mL':len(resps_L),
'ResL':resps_L, 'mean':v_mean, 'sd':v_sd}
with open('model_adj.stan') as f_stan:
sm =
stan.build(f_stan.read(), data = Data)
fit = sm.sample(num_samples = 10000)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
az.plot_trace(i_data, var_names =
('m', 's', 'b'))
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
def CalcMAPEst(samples, a = 0.05,
n_points = 2000):
"""
Calculate
a MAP estimate from a KDE graph on [Lp, Up]
Lp
and Up are 100*a/2 and 100(1-a/2) percentile points of samples
"""
import numpy as np
import scipy.stats as ss
Lp, Up =
np.percentile(samples, [100 * a/2, 100 * (1 - a/2)]) # import numpy as np
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
# import scipy.stats as ss
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
dframe =fit.to_frame()
m_paras = dframe['m']
s_paras = dframe['s']
b_paras = dframe['b']
m_map = CalcMAPEst(m_paras)[0]
m_q1, m_q2, m_q3 =
np.percentile(m_paras, [25, 50, 75])
print(f'parameter m\nmap =
{m_map:.2f}')
print(f'quantiles = {m_q1:.2f},
{m_q2:.2f}, {m_q3:.2f}')
s_map = CalcMAPEst(s_paras)[0]
s_q1, s_q2, s_q3 =
np.percentile(s_paras, [25, 50, 75])
print(f'\nparameter s\nmap =
{s_map:.2f}')
print(f'quantiles = {s_q1:.2f},
{s_q2:.2f}, {s_q3:.2f}')
b_map = CalcMAPEst(b_paras)[0]
b_q1, b_q2, b_q3 = np.percentile(b_paras,
[25, 50, 75])
print(f'\nparameter b\nmap =
{b_map:.2f}')
print(f'quantiles = {b_q1:.2f},
{b_q2:.2f}, {b_q3:.2f}')
az.plot_kde(m_paras)
plt.title(f'm:
map={m_map:.2f}\nq1={m_q1:.3f}, q2={m_q2:.2f}, q3={m_q3:.2f}')
plt.savefig('Fig_m.png')
plt.show()
az.plot_kde(s_paras)
plt.title(f's:
map={s_map:.2f}\nq1={s_q1:.3f}, q2={s_q2:.2f}, q3={s_q3:.2f}')
plt.savefig('Fig_s.png')
plt.show()
az.plot_kde(b_paras)
plt.title(f'b:
map={b_map:.2f}\nq1={b_q1:.3f}, q2={b_q2:.2f}, q3={b_q3:.2f}')
plt.savefig('Fig_b.png')
plt.show()
x_coord = np.arange(np.min(resps),
np.max(resps), 0.01)
U_model = ss.norm.pdf(x_coord, m_map +
b_map, s_map)
L_model = ss.norm.pdf(x_coord, m_map -
b_map, s_map)
plt.hist(resps_U, density = True,
label = 'U')
plt.hist(resps_L, density = True,
label = 'L')
plt.plot([350], [0], c='c', marker =
'o', markersize=20, label = '350')
plt.plot(x_coord, U_model, lw = 5,
label = 'U-model')
plt.plot(x_coord, L_model, lw = 5,
label = 'L-model')
plt.title(f'm = {m_map:.2f}, s = {s_map:.2f}, b = {b_map:.2f}', fontsize = 16)
plt.legend()
plt.tight_layout()
plt.savefig('Fig_distri_model.png')
plt.show()
なお、心理学データ分析と測定の説明を岡本（2021）で、Pythonデータ分析についての説明を岡本（2019）で行っている。