Stanでpsychophysics──階層ベイズモデルで恒常法データを分析する──

00/22 Osaka.Stan #5 – LT資料 2017年7月1日 サ イ コ フ ィ ジ ッ ク ス Stanでpsychophysics ──階層ベイズモデルで恒常法データを分析する── 大阪大学大学院人間科学研究科D2・日本学術振興会 武藤 拓之 (Hiroyuki Muto) Twitter: Web: @mutopsy http://kiso.hus.osaka-u.ac.jp/muto/

• http://kiso.hus.osaka-u.ac.jp/muto/

01/22 自己紹介  武藤 拓之 (むとう ひろゆき) • 大阪大学大学院人間科学研究科D2  研究分野 • 認知心理学 • メインテーマは空間的思考の身体性 • 【宣伝】つい最近論文がJEP:HPPにアクセプトされました！ Muto, H., Matsushita, S., & Morikawa, K. (in press). Spatial perspective taking mediated by whole-body motor simulation. Journal of Experimental Psychology: Human Perception and Performance.  ベイズ歴 • 2017/4/25に自分がベイジアンであることに気が付きました。

02/22 錯視を測りたい  3本の線分は同じ長さ。

03/22 錯視を測りたい  矢羽を付けると長さが変わって見える＝ミュラー・リヤー錯視  実際のところ何倍に見えているのか？ →心理物理学的測定法で測定できる！そのうちの1つが恒常法。 (psychophysical methods)

04/22 恒常法の例 左右の線分，どっちが長い？ 1st 標準刺激 比較刺激 (知覚される大きさを測定したい対象) (ものさし) みぎ！ ひだり！ 2 nd 3 rd なやむ… … このような試行をランダムな順序で繰り返して0/1型のデータを収集する。

05/22 得られるデータの例 (ある個人のデータ)  参加者ごとに プロビット関数をあてはめる プロビット関数は 累積正規分布の ことだよ！

06/22 得られるデータの例 120 (ある個人のデータ)  参加者ごとに プロビット関数をあてはめる  選択率が50%となる長さ ＝主観的等価点 (正規分布の平均値)  この人の場合， は の 120%の長さに見えたと 解釈される  正規分布の標準偏差は ノイズの大きさを表す

07/22 3種類の推定方法を比較してみる  3種類の推定方法で結果を比較 1. 参加者ごとに最尤推定 2. 参加者ごとにベイズ推定 3. 階層ベイズでまとめて推定  使用するデータ Web実験で収集したミュラー・リヤー錯視測定データ  参加者：Twitter・Facebookで募集した延べ26名 (うち1名は逆に回答していたと思われるため分析から除外)  標準刺激：主線が200 px (=100%) の外向図形と内向図形 (今回の分析では外向図形条件のデータのみを用いる)  比較刺激：長さ50─200% (10%刻み) の線分11種類  提示位置をカウンターバランスしランダム順に提示。各セル8試行。 実際の刺激の例→

08/22

1.参加者ごとに最尤推定 (1/1)
Rコード
M <- 1:n
SD <- M
for (i in 1:n){
di <- data.frame(
x = seq(50,150,10),
Rate = Rates[i]
)

#比較刺激の幅
#比較刺激の選択率データのベクトル

fit <- glm(formula = di$Rate ~ di$x, family = binomial(probit))
c1 <- coefficients((fit))
p1 <- c(-c1[1]/c1[2],1/c1[2])
M[i] <- p1[1]
SD[i] <- p1[2]
}
M
SD

#参加者ごとの平均値（等価点）ベクトル
#参加者ごとの標準偏差ベクトル

#glm()で最尤推定

09/22

2.参加者ごとにベイズ推定 (1/2)
Stanコード (1/2)
data {
int N;
int NP;
real Length[N];
int P[N];
int <lower=0, upper=1> SelectC[N];
}

#全部で何試行か。今回は2200 (= 88試行 × 25人)。
#参加者の数。今回は25。
#比較刺激の水準 (単位は%) = {50, 60, …, 140, 150}
#参加者番号 = {1, 2, …, 24, 25}
#比較刺激を選んだか否かのベクトル = {0, 1, 1, …}

parameters{
real<lower=0, upper=200> mu[NP];
real<lower=0, upper=50> sigma[NP];
}

#参加者一人ひとりの平均
#参加者一人ひとりの標準偏差

母数の範囲を指定しないと
うまく収束しないよ！

transformed parameters{
real<lower=0, upper=1> p[N];
for (n in 1:N){
p[n] = normal_cdf(Length[n], mu[P[n]], sigma[P[n]]);
}
}

#累積正規分布の確率を取得

10/22

2.参加者ごとにベイズ推定 (2/2)
Stanコード (2/2)
model{
for (n in 1:N){
SelectC[n] ~ bernoulli(p[n]);
}
}
generated quantities{
real<lower=0, upper=200> mu0;
real<lower=0, upper=50> sigma0;
mu0 = 0;
sigma0 = 0;
for (np in 1:NP){
mu0 = mu0 + mu[np];
sigma0 = sigma0 + sigma[np];
}
mu0 = mu0 / NP;
sigma0 = sigma0 / NP;
}

#ベルヌーイ分布から反応データ (0/1) が生じるモデル

#参加者間の平均 (mu0とsigma0) を生成

11/22

3.階層ベイズでまとめて推定 (1/2)
Stanコード (1/2)
data {
int N;
int NP;
real Length[N];
int P[N];
int <lower=0, upper=1> SelectC[N];
}

parameters{
real<lower=0, upper=200> mu[NP];
real<lower=0, upper=50> sigma[NP];
real<lower=0, upper=200> mu0;
real<lower=0, upper=50> sigma0;
real<lower=0, upper=200> s_mu0;
real<lower=0, upper=200> s_sigma0;
}

#全体平均
#全体標準偏差
#参加者による平均のばらつき
#参加者による標準偏差のばらつき

12/22

3.階層ベイズでまとめて推定 (2/2)
Stanコード (2/2)
transformed parameters{
real<lower=0, upper=1> p[N];
for (n in 1:N){
p[n] = normal_cdf(Length[n], mu[P[n]], sigma[P[n]]);
}
}
model{
for (np in 1:NP){
mu[np] ~ normal(mu0,s_mu0);
#参加者ごとの平均を得る
sigma[np] ~ normal(sigma0,s_sigma0);
#参加者ごとの標準偏差を得る
}
for (n in 1:N){
SelectC[n] ~ bernoulli(p[n]);
参加者一人ひとりの平均が
}
Normal(全体平均, s_mu) に
}

従うと仮定したよ！
標準偏差も同様！

13/22 ある参加者の結果 階層ベイズ 個別にベイズ推定 個別に最尤推定 ※推定値としてEAPを使用。

14/22 みんなの結果 (1/2) 個別に最尤推定 個別にベイズ推定 階層ベイズ

15/22 みんなの結果 (2/2) 個別に最尤推定 個別にベイズ推定 階層ベイズ

16/22 推定法による母平均の違い  点推定値はほとんど同じ  区間推定幅が最も狭いのは 個別にベイズ推定したとき  ベイズ推定で過大推定傾向 (分布の歪みが原因？）  区間推定幅が最も狭いのは 個別にベイズ推定したとき エラーバーは95%信頼区間または95%確信区間。

17/22 データを減らしたらどうなるか  データを半分にしてみる Rコード #ランダムに欠損を作る V <- 1:nrow(d) Smpl <- sample(V,nrow(d)*0.5) d_sps <- d[Smpl,] #何パーセントにするか ↓ある参加者の結果 ※条件によって分母が異なる 階層ベイズ 個別にベイズ推定 個別に最尤推定

18/22 欠損データのみんなの結果 (1/2) 個別に最尤推定 個別にベイズ推定 階層ベイズ

19/22 欠損データのみんなの結果 (2/2) 個別に最尤推定 個別にベイズ推定 階層ベイズ

20/22 推定法による母平均の違い(欠損データ)  どの推定法でも あまり変わらない。  データを削っても 意外と正確に推定できている。  ベイズ推定で過大推定傾向  区間推定幅が最も狭いのは 個別にベイズ推定したとき エラーバーは95%信頼区間または95%確信区間。

21/22 まとめと考察  全体平均の点推定値はどの推定方法でもあまり変わらなかった。  データ数を減らしても割と正確に全体平均を推定できた。 →ミュラー・リヤー錯視が頑健すぎるから？  全体平均の区間推定幅は個別にベイズ推定した時に最も狭くなった。 →なぜかは不明。 →今回のデータに特有の結果なのか一般的な傾向なのか要検討。  全体標準偏差の推定値は推定方法に影響された (ベイズで過大推定傾向)。 →分布の歪みが原因？参加者ごとのMAP推定値の平均をとるべき？  全体標準偏差の区間推定幅は個別にベイズ推定した時に最も狭くなった。 →なぜかは不明。  結局どの推定方法が一番いいの？ →見かけ上は個別にベイズ推定した時に精度が高かった。 →最初に想定するモデル（分析者の信念）とデータの特徴次第？

22/22 どんなにくるしくても おいしいものたべて ベイズしたらなおるよ！ なおるよ！