×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
Andrew Yonelinasの再認記憶2過程信号検出理論 (dual-process signal-detection model) の推定をRでやってみる
本人のホームページ: http://psychology.ucdavis.edu/labs/Yonelinas/index.html
ExcelソルバーでのサンプルをRでやっただけ
# データはExcelファイル (DPSDSSE) に従う
hit <- c(0.87, 0.75, 0.61, 0.47, 0.34) # recognition課題のヒット率
fa <- c(0.71, 0.5, 0.31, 0.14, 0.05) # フォールス・アラーム
est.v.start <- c(rep(0.5, 7)) # 推定するのは7つ。ro, d', cが5つ。初期値を全部0.5にする
# 推定式から出るヒット率、フォールス・アラーム率の予測値と、それぞれの実測値との二乗誤差を最小化するパラメータを出す
# 最適化のための関数定義
fc <- function(x){
ro <- x[1]
rn <- 0
dp <- x[2]
cv <- x[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
difhit <- (hit-prehit)^2
diffa <- (fa-prefa)^2
sum(difhit+diffa)
}
# optim関数で最適化
(opres <- optim(par=est.v.start, fn=fc, method="BFGS"))
est.v <- opres$par
names(est.v) <- c("ro", "dprime", "c1","c2","c3","c4","c5")
est.v # 最終的な推定値
# 改めて予測を出す
ro <- est.v[1]
rn <- 0
dp <- est.v[2]
cv <- est.v[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
pre.o <- cbind(prehit=prehit, hit, prefa=prefa, fa) # 予測、実測をまとめる
rownames(pre.o) <- NULL
pre.o
fc(est.v) # 最終的な二乗誤差
# プロット
plot(fa, hit, xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i")
# 曲線を描く。
# 適当なx軸データ (fa) をつくり、推定値に基づいてhitを予測する。その上でラインを引く。
plfa <- c(0.001,0.01,0.02,0.035,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99999)
plcv <- qnorm(plfa, mean=0, sd=1)
plhit <- plfa+ro+(((1-ro)*(pnorm(plcv, -dp, oldvar)))-((1-rn)*(pnorm(plcv, 0, 1))))
sum(plhit>1) # 1を超えるかどうか調べる
lines(spline(plfa, plhit))
abline(0,1) # 対角線
# もっとうまい描きかたがあるだろうな…
本人のホームページ: http://psychology.ucdavis.edu/labs/Yonelinas/index.html
ExcelソルバーでのサンプルをRでやっただけ
# データはExcelファイル (DPSDSSE) に従う
hit <- c(0.87, 0.75, 0.61, 0.47, 0.34) # recognition課題のヒット率
fa <- c(0.71, 0.5, 0.31, 0.14, 0.05) # フォールス・アラーム
est.v.start <- c(rep(0.5, 7)) # 推定するのは7つ。ro, d', cが5つ。初期値を全部0.5にする
# 推定式から出るヒット率、フォールス・アラーム率の予測値と、それぞれの実測値との二乗誤差を最小化するパラメータを出す
# 最適化のための関数定義
fc <- function(x){
ro <- x[1]
rn <- 0
dp <- x[2]
cv <- x[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
difhit <- (hit-prehit)^2
diffa <- (fa-prefa)^2
sum(difhit+diffa)
}
# optim関数で最適化
(opres <- optim(par=est.v.start, fn=fc, method="BFGS"))
est.v <- opres$par
names(est.v) <- c("ro", "dprime", "c1","c2","c3","c4","c5")
est.v # 最終的な推定値
# 改めて予測を出す
ro <- est.v[1]
rn <- 0
dp <- est.v[2]
cv <- est.v[3:7]
oldvar <- 1
prehit <- ro+((1-ro)*(pnorm(cv, -dp, oldvar)))
prefa <- (1-rn)*(pnorm(cv, 0,1))
pre.o <- cbind(prehit=prehit, hit, prefa=prefa, fa) # 予測、実測をまとめる
rownames(pre.o) <- NULL
pre.o
fc(est.v) # 最終的な二乗誤差
# プロット
plot(fa, hit, xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i")
# 曲線を描く。
# 適当なx軸データ (fa) をつくり、推定値に基づいてhitを予測する。その上でラインを引く。
plfa <- c(0.001,0.01,0.02,0.035,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99999)
plcv <- qnorm(plfa, mean=0, sd=1)
plhit <- plfa+ro+(((1-ro)*(pnorm(plcv, -dp, oldvar)))-((1-rn)*(pnorm(plcv, 0, 1))))
sum(plhit>1) # 1を超えるかどうか調べる
lines(spline(plfa, plhit))
abline(0,1) # 対角線
# もっとうまい描きかたがあるだろうな…
PR
Comment
Trackback
Trackback URL
Comment form