反応時間の確率分布：実験と分析

反応時間の確率分布

実験と分析

反応時間の分布は正規分布ではないことが報告されている（Palmerら、2011）。Palmerらは、５つの確率分布Ex-Gaussian、Ex-Wald、(shifted) Gamma、(shifted) Weibull、Gaussianを用いてデータ分析を行い、Ex-Gaussian、Ex-Wald、(shifted) Gammaの３つは同じようにデータへの当てはまりがよかったが、他の２つ、(shifted) Weibull、Gaussian、はよくなかったことを報告している。

なお、Okamoto (1981, 1982, 1980) は、単純反応時間において２状態の区別を報告している。反応時間のメカニズムの分析では状態の区別が必要であるが、本ウェブサイトではこの区別を行わず、単純に１つの分布を当てはめている。

本ウェブサイトでは、単純反応時間測定実験のプログラムとその分析プログラムを用意したが、分析結果例を下図に示す。

図１

緑色のヒストグラムがデータであり、青の曲線がshifted gamma 分布（以後、簡単にガンマ分布という）、橙色の曲線がGauss分布（以後、正規分布という）である。ガンマ分布は、Stanで簡単にプログラムが書けるので、本ウェブサイトでは３つの分布、Ex-Gaussian、Ex-Wald、(shifted) Gamma、の代表としてガンマ分布を用いる。図１に示されているように、ガンマ分布の当てはまりはよいが、正規分布の当てはまりは悪い。Ex-GaussianとEx-WaldのStanスクリプトは、別のウェブサイトに用意している。

本ウェブサイトでは、まず単純反応時間測定実験の実行例を示し、次に得られたデータの分析例を示す。その後に、実験を行うプログラム例、およびデータ分析のプログラム例を示す。プログラムファイルおよびサンプルデータファイルは圧縮ファイルfilesRT.zipにまとめた。ダウンロード展開して、自由に利用されたい。

なお、WindowsにおけるPsychoPy3およびPyStanのインストールについての補足説明を最後に用意している。

上記「補足説明」の方法で2023.02.10にPsychoPyのインストールを行うとうまくいかなかった。仮想環境に簡単なpipコマンドでインストールするとできたので、その説明を用意した。

単純反応時間測定実験

単純反応時間を測定する実験のプログラムSimpleRT_Exp.pyをPsychoPy 3のCoderで作成した。PsychoPy 3は、スタンドアロン版をインストールしたので、他のPythonのインストールとは独立した環境でPsychoPy 3のPythonを利用することができる。SimpleRT_Exp.pyのリストは、このウェブサイトの後半にリスト１として載せた。

なお、Windows10においてPsychoPy3の再インストールのためにStandalone版をダウンロードすると「・・・はお使いのデバイスに問題を引き起こす可能性があるため、ブロックされました」という警告が出た（2020.09）。Windows10にAnacondaをインストールして（簡単にインストールできる）、AnacondaにおいてPsychopy3をインストールしてみるとできた。このAnacondaでPsychoPy3をインストールする方法についてはこのウェブサイトの最後に説明を用意した。注意：　この「最後に説明を用意した」とある方法では、2023.02.10にはうまくいかなかったが、別の方法ではインストールできた。

プログラムSimpleRT_Exp.pyをPsychoPy 3のCoderで開いて実行すると、図２のフォームが表示される。

図２、

キーZを押すと、図２のウィンドウ内の文字は消え、図３の状態になる。

図３

数秒後（毎試行ごとにランダムに変わる）、図４のように文字Mが表示される。

図４

文字Mが表示されたら、キーMを押す。文字の表示からキーMを押すまでの時間が反応時間として記録され、次の試行に移る。

全ての試行（サンプルプログラムでは100試行）が終了すると、データの保存用ファイルの名前の設定を求めるウィンドウが表示される。ファイルへは、テキストファイルとして出力されるので、ファイル拡張子は、「.txt」としておくとよい。図５では、「MyData.txt」と設定されている。

図５

保存されたデータの分析例を、次に説明する。

データの分析

データ分析プログラムAnalysis.pyは、上の実験制御プログラムの実行環境であるスタンドアロン版PsychoPy 3とは別にインストールされたPython3で実行した。PsychoPy 3は、ディレクトリC:\Program Files\PsychoPy3の下にインストールされているが（デフォルトでのインストール先）、データ分析用Python3は、ディレクトリC:\Python374の下にインストールしている（Python 3のインストール時に作成指定）。また、PyStanのインストールについてはこのウェブサイトで説明しているが、本ウェブサイトの最後にも説明を用意している。

データにガンマ分布を適用してベイズ分析を行った。ガンマ分布は次式で与えられる。

ガンマ分布に対するStanスクリプト（shiftedGamma.stan）は以下に示すように簡単である。

data {

int N;

real RT[N];

real max_a;

real max_b;

real max_x;

}

parameters {

real<lower = 0.0, upper = max_a> param_a;

real<lower = 0.0, upper = max_b>param_b;

real<lower = 0.0, upper = max_x> param_x;

}

model {

param_a ~ uniform(0.0, max_a);

param_b ~ uniform(0.0, max_b);

param_x ~ uniform(0.0, max_x);

for(i in 1:N){

RT[i] - param_x ~ gamma(param_a, param_b);

}

正規分布による分析用のStanスクリプト（Normal.stan）は以下のようになる。

data {

int N;

real RT[N];

real min_mu;

real max_mu;

real max_sgm;

}

parameters {

real<lower = min_mu, upper = max_mu> mu;

real<lower = 0.0, upper = max_sgm> sgm;

}

model {

mu ~ uniform(min_mu, max_mu);

sgm ~ uniform(0.0, max_sgm);

for(i in 1:N){

RT[i] ~ normal(mu, sgm);

}

これらのStanスクリプトを用いた分析プログラムAnalysis.pyは、リスト２に載せてある。PythonにおけるStan（PyStan）の使い方の説明は、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を参照されたい。

プログラムAnalysis.pyを実行すると、図６のように入力データファイル名を聞いてくる。

図６

図６では、図５で設定したデータ出力ファイル名を入力ファイルとして設定している。ファイル名の設定後、読み込んだデータの出力が行われ、続いてStanによるガンマ分布の分析が始まる。StanのMCMCサンプリングが終了すると、まずパラメータの事後分布が表示される（図７）。

図７

図中、タイトルのlen()は、サンプル数を表している。図７のウィンドウを閉じると、の事後分布が表示される（図８）。

図８

図８のウィンドウを閉じると、次にの事後分布が表示される（図９）。

図９

図９のウィンドウを閉じると、正規分布による分析が始まる。Stanのサンプリングが終了すると、の事後分布が表示される（図１０）。

図１０

図１０のウィンドウを閉じると、の事後分布が表示される（図１１）。

図１１

図１１のウィンドウを閉じると、反応時間（RT）データのヒストグラムとともに、上のサンプリングの中央値を各パラメータの点推定値とするガンマ分布と正規分布のグラフが表示される（図１２）。

図１２

図１２のウィンドウを閉じると、プログラムは終了する。

プログラムソースリスト（プログラムファイルは圧縮ファイルfilesRT.zipとしてまとめた）

リスト１．単純反応時間測定実験のプログラムSimpleRT_Exp.py（Coder/PsychoPy 3）

import math, numpy, random

from psychopy import core, event, visual, gui

import psychopy.sound.backend_ptb as pss

myWin = visual.Window(size=(500, 300), pos=None, color='black',

colorSpace='rgb',

rgb=None, dkl=None, lms=None, fullscr=False, allowGUI=True,

monitor=None, bitsMode=None, winType=None, units='pix',

gamma=None, blendMode='avg', screen=0, viewScale=None,

viewPos=None, viewOri=0.0, waitBlanking=True, allowStencil=False,

multiSample=False, numSamples=2, stereo=False, name='window1',

checkTiming=True, useFBO=False, useRetina=True, autoLog=True,

gammaErrorPolicy='raise', bpc=(8, 8, 8), depthBits=8,

stencilBits=8, backendConf=None)

myClock = core.Clock()

print('size = ', myWin.size)

mySound = pss.SoundPTB(value = 2000, secs = 0.5, volume = 1.0)

latencies = []

for i in range(5):

t = 2.0 * (1.2**i)

latencies.append(t)

g_latencies = []

for i in range(40):

numpy.random.shuffle(latencies)

g_latencies = g_latencies + latencies

print('g_latnecies = ', g_latencies)

n_trials = 100 # Total number of trials

RTs = numpy.empty(n_trials)

Intvls = numpy.empty(n_trials)

information = visual.TextStim(myWin, text = '', height = 24, color = 'white')

t = 1

pos_g_latencies = 0

while True:

# 第t試行

mySound.play()

msg = '{0}/{1}'.format(t, n_trials) + \

'\nPress the Z key, \nwhen youv are ready.' + \

'\n\nAfter some seconds, \nthe character M will be presented.' + \

'\n\nPress the M key as fast as possible \nwhen M appears.'

information.setText(msg)

information.setHeight(24)

information.draw()

myWin.flip()

event.clearEvents()

while True:

pressedList = event.getKeys(keyList = ['z'])

if len(pressedList) > 0:

if pressedList[0] == 'z': # 試行開始

event.clearEvents()

break

myWin.flip()

t1 = myClock.getTime()

flying = False

while True:

pressedList = event.getKeys(keyList = ['m'])

if len(pressedList) > 0:

if pressedList[0] == 'm': # フライング

flying = True

information.setText('Flying')

information.draw()

myWin.flip()

core.wait(2)

event.clearEvents()

pos_g_latencies += 1

break

event.clearEvents()

if myClock.getTime() > t1 + g_latencies[pos_g_latencies]:

break

if flying == True:

continue

information.setText('M')

information.setHeight(48)

information.draw()

myWin.flip()

event.clearEvents()

t1 = myClock.getTime() # 刺激提示時刻

t2 = None

while True:

pressedList = event.getKeys(keyList = ['m'])

if len(pressedList) > 0:

if pressedList[0] == 'm':

t2 = myClock.getTime() # 反応時刻

break

myWin.flip()

core.wait(2)

RTs[t-1] = t2 - t1 # 反応時間

Intvls[t-1] = g_latencies[pos_g_latencies]

t += 1

pos_g_latencies += 1

if pos_g_latencies >= len(g_latencies):

pos_g_latnecies = 0

if t > n_trials:

break

# 反応時間の出力と保存

myWin.close()

for i in range(n_trials):

print('RT: ', RTs[i], ' Intvl: ', Intvls[i])

flnm = gui.fileSaveDlg(prompt = u'出力ファイル名(*.txt)')

print('file name = ', flnm)

fout = open(flnm, 'w')

fout.write('/\n')

fout.write('{0:>10s}{1:>10s}\n'.format('RT', 'Interval'))

for i in range(n_trials):

fout.write('{0:>10.5f}{1:>10.5f}\n'.format(RTs[i], Intvls[i]))

fout.write('/\n')

fout.close()

print(flnm, ' was saved.')

core.quit()

リスト２．データ分析用PythonスクリプトAnalysis.py

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import scipy.stats as ss

import pystan

# データの読み込み

flnm = input('Data file name = ')

fin = open(flnm, 'r')

data = fin.readlines()

fin.close()

pos = 0

while True:

if len(data[pos]) > 0:

if data[pos][0] == '/':

break

pos += 1

RTs = []

Intvls = []

pos += 2

while True:

if data[pos][0] == '/':

break

slist = data[pos].split()

RTs.append(float(slist[0]))

Intvls.append(float(slist[1]))

pos += 1

for x, y in zip(RTs, Intvls):

print(x, ' ', y)

RTs = np.array(RTs)

xi_hat = RTs.min() # シフト値の最大値

print('xi_hat = ', xi_hat)

RTs_xi = RTs - xi_hat

# モーメント法に依るパラメータ値の推定

beta_hat = RTs_xi.mean() / RTs_xi.var()

print('beta_hat = ', beta_hat)

alpha_hat = beta_hat * RTs_xi.mean()

print('alpha_hat = ', alpha_hat)

max_a = alpha_hat * 100.0

max_b = beta_hat * 100.0

max_x = xi_hat

Data = {'N': len(RTs), 'RT': RTs, 'max_a': max_a, 'max_b': max_b, 'max_x': max_x}

sm = pystan.StanModel(file = 'shiftedGamma.stan') # ガンマ分布

fit = sm.sampling(data = Data, n_jobs = 1)

param_a = fit['param_a']

param_b = fit['param_b']

param_x = fit['param_x']

plt.hist(param_a)

plt.xlabel(r'$\alpha$', fontsize = 14)

plt.title(r'len($\alpha$) = {}'.format(len(param_a)), fontsize = 18)

plt.show()

plt.hist(param_b)

plt.xlabel(r'$\beta$', fontsize = 14)

plt.title(r'len($\beta$) = {}'.format(len(param_b)), fontsize = 18)

plt.show()

plt.hist(param_x)

plt.xlabel(r'$\xi$', fontsize = 14)

plt.title(r'len($\xi$) = {}'.format(len(param_x)), fontsize = 18)

plt.show()

# 中央値

med_a = np.percentile(param_a, 50)

med_b = np.percentile(param_b, 50)

med_xi = np.percentile(param_x, 50)

print('med_a = ', med_a, ' med_b = ', med_b, ' med_xi = ', med_xi)

max_sgm = (RTs.var()**0.5) * 10.0

min_mu = RTs.mean() - max_sgm

max_mu = RTs.mean() + max_sgm

Data = {'N': len(RTs), 'RT': RTs, 'min_mu': min_mu, 'max_mu': max_mu, 'max_sgm': max_sgm}

sm1 = pystan.StanModel(file = 'Normal.stan') # 正規分布

fit1 = sm1.sampling(data = Data, n_jobs = 1)

param_mu = fit1['mu']

param_sgm = fit1['sgm']

plt.hist(param_mu)

plt.xlabel('$\mu$', fontsize = 14)

plt.title(r'len($\mu$) = {}'.format(len(param_mu)))

plt.show()

plt.hist(param_sgm)

plt.xlabel('$\sigma$', fontsize = 14)

plt.title(r'len($\sigma$) = {}'.format(len(param_sgm)))

plt.show()

# 中央値

med_mu = np.percentile(param_mu, 50)

med_sgm = np.percentile(param_sgm, 50)

print('med_mu = ', med_mu, ' med_sgm = ', med_sgm)

rt = np.arange(0, 0.5, 0.001)

# ガンマ分布

pdf_values = ss.gamma.pdf(rt, med_a, med_xi, 1/med_b)

plt.plot(rt, pdf_values, linewidth = 3, label = 'gamma(ξ,α,β)')

# 正規分布

pdf1_values = ss.norm.pdf(rt, med_mu, med_sgm)

plt.plot(rt, pdf1_values, linewidth = 3, label = 'normal(μ,σ)')

# データ

plt.hist(RTs, density = True, label = 'RT')

plt.xlabel('Time', fontsize = 14)

plt.title('ξ = {0:<.3f}, α = {1:<.3f}, β = {2:<.3f}\n'.format(med_xi, med_a, med_b) +\

'μ = {0:<.3f}, σ = {1:<.3f}'.format(med_mu, med_sgm),

fontsize = 16)

plt.legend(fontsize = 14)

plt.show()

引用文献

Palmer, E. M., and Horowitz, T. S., Torralba, A., and Wolfe, J. M. (2011). What are the shapes of response time distributions in visual search? Journal of Experimental Psychology: Human Perception and Performance, 37, 58-71.

WindowsにおけるPsychoPy3とPyStanのインストールについて

Windows10においてMicrosoft EdgeでPsychoPy3のStandalone版をインストールするためにダウンロードを試みると「StandalonePsychoPy….exeはお使いのデバイスに問題を引き起こす可能性があるため、ブロックされました。」という警告が出た（2020.09）。そのため、以下のようにAnacondaでインストールを試みるとインストールできた（2023.02.10に行うと、うまくいかなかったので、pipコマンドで直接インストールするとできた。この説明も後に示した）。また、PyStanもAnacondaでは簡単にインストールできるが、これはPsychoPyの説明の後で解説する。