2026.03改訂
単純反応時間の実験プログラムをPyGameを利用して書いてみた。PyGameの利用は簡単であり、コマンド
pip install pygame
でインストールできる（Windowsの場合）。
プログラムファイルなどは圧縮ファイルsrt_pygame_files.zipにまとめた。ダウンロード展開すれば、自由に利用できる。
プログラムは、リスト１に示す通りである。
実験において表示する日本語文は、予めお絵かきソフト（ペイント）で作成しておき、実験においては画像ファイルとして読み込んで表示した。
例えば、説明文は次のコードでイメージimg_explnとしている。
img_expln = pg.image.load('expln.png')
rect_expln = img_expln.get_rect()
これを、次のコードで表示している。
screen.blit(img_expln,
((800-rect_expln.width)//2,(500-rect_expln.height)//2))
pg.display.flip()
表示すると、図１のようになる。
図１
実験は、図１の状態で「Gキー」を押すと反応時間計測試行へ進む。
各試行は「ヨーイ」文の表示で始まる（図２）。
図２
画面左上には、現在の試行数と実験の総試行数が表示されている。
「Fキー」を押すと、「ヨーイ」文は消える（図３）。
図３
しばらく後、「ドン」文が表示される（図４）。
図４
「Jキー」を押すと、「ドン」文は消える。「ドン」文の表示から「Jキー」の押下までが反応時間として記録される。
「ドン」が消えてから暫く後に「ヨーイ」文が表示され、図２に戻る。「Fキー」を押して、上の手順を繰り返す。
総試行が終了すると、出力ファイル名の設定が求められる。
出力ファイル名(*.txt)
= Results.txt
出力ファイル名（テキストファイル名）を設定すると、プログラムの終了である。
出力ファイルに格納されたデータの分析プログラムをリスト２、リスト２で用いるStanスクリプトをリスト３のように用意した。スクリプトファイルなどはsrt_pygame_files.zipにまとめた。
リスト３では、３パラメータ対数正規分布を用いているが、３パメータ対数正規分布についてはこのウェブサイトで説明している。反応時間の分布に３パラメータ対数正規分布がよく当てはまることがこのウェブサイトで説明されている。
プログラム（スクリプト）は、CmdStanPy対応である。CmdStanPyの簡単な説明をこのウェブサイトで行った。
CmdStanPyのインストールされた環境において、リスト２のプログラムを実行すると、入力ファイル名を聞いてくる。
(stan) ****/srt_pygame_files$ python anal_srt.py
Data file
name = Results.txt
上の実験で保存したファイル名を設定すると、読み込まれ、リスト３のStanスクリプトによる分析が行われる。上の例では、入力データファイルは、リスト２およびリスト３と同じフォルダに置かれている。
MCMCサンプリングが終了すると、トレース図が表示される（図A.1）。
図A.1
図A.1のWindowsの右上のX印をクリックして閉じると、図A.2のグラフが表示される。
図A.2
３パラメータ対数正規分布のパラメータの事後分布のグラフである。平均値パラメータ、中央値パラメータ、標準偏差パラメータの中央値が各分布の代表値としてタイトル欄に表示さている。
図A.2のWindowを閉じると、図A.3のグラフが表示される。
図A.3
データのヒストグラムと、事後分布の中央値をパラメータ値とする３パラメータ対数正規分布のグラフである。タイトル欄に表示されている値は、３パラメータ対数正規分布に対応する正規分布のパラメータである。
図A.3のWindowを閉じると、スクリプトの実行終了である。
import pygame
as pg
from pygame
import font
from pygame.time
import get_ticks as time_ms
import numpy.random
as nr
pg.init()
font.init()
screen = pg.display.set_mode(size=(800,
500)) # 画面の大きさ
screen.fill((0,0,0))
# 黒で塗りつぶす
my_font = pg.font.SysFont(None,32)
# フォントオブジェクト
#
# 説明文のイメージ
#
img_expln = pg.image.load('expln.png')
rect_expln = img_expln.get_rect()
#
# 「ヨーイ」文のイメージ
#
img_rdy = pg.image.load('youi.png')
rect_rdy = img_rdy.get_rect()
#
#
「ドン」文のイメージ
#
img_don = pg.image.load('don.png')
rect_don = img_don.get_rect()
#
#
説明文の表示
#
screen.blit(img_expln,
((800-rect_expln.width)//2,(500-rect_expln.height)//2))
pg.display.flip()
ck = False
while not ck:
for e in pg.event.get():
if e.type == pg.KEYDOWN:
if e.key == pg.K_g: # キーGが押されたかのチェック
ck = True
break
#
#
１秒間のブラックアウト
#
screen.fill((0,0,0))
pg.display.flip()
pg.time.wait(1000)
w_times = [2, 3, 4.5, 6] # 「ドン」が提示されるまでの時間のリスト（秒）
n_trials = 50
# 総試行数
times = []
# 反応時間を格納するリスト
for i_trial
in range(n_trials):
#
# 試行数と「ヨーイ」文の表示
#
trial_str
= f'Trial {i_trial+1}/{n_trials}'
img_trial
= my_font.render(trial_str,
True, (255,255,255), (0, 0, 0))
screen.blit(img_trial, (10, 10))
screen.blit(img_rdy, ((800-rect_rdy.width)//2,
(500-rect_rdy.height)//2))
pg.display.flip()
ck = False
while not ck:
for
e in pg.event.get():
if e.type == pg.KEYDOWN:
if e.key == pg.K_f: # 「キーF」が押されたか？
ck = True
break
#
# ブラックアウト
#
screen.fill((0,0,0))
pg.display.flip()
w_t
= nr.choice(w_times) *
1000 # ブラックアウトの時間の選択
c_t
= time_ms()
while time_ms()
< c_t + w_t:
# w_tミリ秒経過したか？
pg.event.get()
# フライング対策
#
# 「ドン」文の提示
#
screen.blit(img_don, ((800-rect_don.width)//2,
(500-rect_don.height)//2))
pg.display.flip()
t0 = time_ms()
# 反応時間計測の始点
ck = False
while not ck:
for
e in pg.event.get():
if e.type == pg.KEYDOWN:
if e.key == pg.K_j: # 「キーJ」が押されたか？
ck = True
break
t1 = time_ms()
# 反応時間計測の終点
#
#
ブラックアウト
#
screen.fill((0,0,0))
pg.display.flip()
times.append(t1-t0) # 反応時間をリストtimesに追加
pg.time.wait(2000) # ２秒待ち
font.quit()
pg.quit()
fl_nm = input('\n出力ファイル名(*.txt) = ')
fout = open(fl_nm,
'w')
for v in times:
fout.write(f'{v}\n')
fout.close()
print(f'\n{fl_nm}
was saved.\n')
import numpy
as np
import matplotlib.pyplot
as plt
import seaborn as sb
from cmdstanpy
import CmdStanModel
import pandas as pd
import arviz
as az
fl_nm = input('Data file name = ')
f_in = open(fl_nm)
data = f_in.read().splitlines()
print(data)
data = [float(v) for v in data]
print(data)
Data = {'n': len(data),
'x':data}
sm = CmdStanModel(stan_file='lognormal3p.stan')
fit = sm.sample(data=Data)
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# inf_dataを用いてのトレースプロット
az.plot_trace(inf_data, var_names=['mean', 'mode', 'sd',
'theta'],
figsize=(12,8))
plt.tight_layout()
plt.show()
fit_dframe = fit.draws_pd()
# PandasのDataFrame型に変換
print(fit_dframe.keys())
plt.figure(figsize=(12,5))
plt.subplot(121)
modes = fit_dframe['mode']
mode_med = np.median(modes)
sb.kdeplot(modes, label='mode')
means = fit_dframe['mean']
means_med = np.median(means)
sb.kdeplot(means, label='mean')
plt.xlabel('mode and mean', fontsize=14)
plt.title(f'mode(med.)={mode_med:.1f},'
+
f'
mean(med.)={means_med:.1f}', fontsize=16)
plt.legend()
plt.subplot(1,2,2)
sds = fit_dframe['sd']
sds_med = np.median(sds)
sb.kdeplot(sds)
plt.xlabel('sd', fontsize=14)
plt.title(f'sd={sds_med:.3f}',
fontsize=16)
plt.show()
def lognormal_pdf(x,
mu, sgm):
v1 = sgm
* x * ((2*np.pi)**0.5)
v2 = np.exp(-((np.log(x)
- mu)**2) / (2 * (sgm**2)))
return v2 / v1
x_min = np.min(data)
x_max = np.max(data)
mu_med = np.median(fit_dframe['mu'])
sgm_med = np.median(fit_dframe['sgm'])
theta_med = np.median(fit_dframe['theta'])
xcoord = np.linspace(x_min, x_max, 1000)
ycoord = lognormal_pdf(xcoord - theta_med, mu_med, sgm_med)
plt.title(f'mu={mu_med:.2f}, sgm={sgm_med:.3f}, theta={theta_med:.3f}',
fontsize=14)
plt.hist(data, bins=15, density=True,
label='data')
plt.plot(xcoord, ycoord, lw=3, label='LogNormal')
plt.xlabel('RT(msec)', fontsize=14)
plt.legend()
plt.show()
リスト３ ベイズ分析用Stanスクリプト（lognormal3p.stan）
data {
int n;
array[n] real x;
}
transformed data {
real min_x;
min_x
= min(x);
}
parameters {
real mu;
real<lower=0.0001,
upper=10000.0> sgm;
real<lower=0.0, upper=min_x> theta;
}
model {
mu ~ normal(0.0, 10000);
sgm
~ uniform(0.0001, 10000.0);
for (i
in 1:n) {
x[i] - theta ~ lognormal(mu, sgm);
}
}
generated quantities {
real mean;
real mode;
real sd;
mean = theta + exp(mu +
0.5*(sgm^2));
mode = theta + exp(mu -
(sgm^2));
sd
= (exp(2*mu + (sgm^2)) * (exp(sgm^2) - 1.0))^0.5;
}