×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
Van Selst, M., & Jolicoeur, P. (1994). A solution to the effect of sample size on outlier elimination. The quarterly journal of experimental psychology, 47A, 631-650.
より。
結論としては基準変動型のいずれかがお勧めらしい
## 仮想反応時間データ生成
rt <- round(sample(rnorm(10000, mean=300, sd=10)+(300*rexp(10000)), 20))
hist(rt)
win.graph(); plot(rt)
mean(rt); sd(rt)
range(rt)
# 基準変動型手続きに使うデータ
ssz <- c(4,5,6,7,8,9,10,11,12,13,14,15,20,25,30,35,50,100) # 抽出するサンプルの数
nrv <- c(1.458,1.68,1.841,1.961,2.050,2.12,2.173,2.22,2.246,2.274,2.31,2.326,2.391,2.410,2.4305,2.45,2.48,2.50) # 基準変動型非再帰手続き用
mrv <- c(8, 6.20, 5.30, 4.80, 4.475,4.25,4.11,4.00,3.92,3.85,3.80,3.75,3.64,3.595,3.55,3.54,3.51,3.50) # 基準変動型修正再帰手続き用
## 表にしてみる
tb <- rbind(nrv, mrv)
colnames(tb) <- ssz
tb
## 線形補間で予測値を得る関数を生成
nrvfun <- approxfun(ssz, nrv) # 基準変動型非再帰手続き用
mrvfun <- approxfun(ssz, mrv) # 基準変動型修正再帰手続き用
## 基準変動型非再帰手続き (Non-recursive with moving criterion)
mrt <- mean(rt)
sdrt <- sd(rt)
sszrt <- length(rt)
crco <- nrvfun(sszrt) # 線形補間関数による基準値
cv <- c(mrt-(crco*sdrt), mrt+(crco*sdrt))
rt2 <- rt[cv[1]<=rt & rt<=cv[2]]
mean(rt); sd(rt); length(rt)
mean(rt2); sd(rt2); length(rt2) # 最終出力
range(rt); range(rt2); cv; crco
## 基準変動型修正再帰手続き (Modified recursive with moving criterion)
rttmp <- rt[rt<max(rt)]
mrttmp <- mean(rttmp)
sdrttmp <- sd(rttmp)
sszrt <- length(rt)
crco <- mrvfun(sszrt) # 線形補間関数による基準値
cv <- c(mrttmp-(crco*sdrttmp), mrttmp+(crco*sdrttmp))
minv <- min(rt); maxv <- max(rt)
ucv.l <- cv[1]; ucv.u <- cv[2]
svct <- rt
while((minv < ucv.l | maxv > ucv.u)&(length(svct)>3)) {
if (minv < ucv.l) {svct<- svct[svct>minv]}
if (maxv > ucv.u) {svct<- svct[svct<maxv]}
cat(ucv.l, ucv.u, crco, "\n") # 最後に使われた基準値とカットオフポイントを返す
if (length(svct) < 2) break
maxv <- max(svct)
minv <- min(svct)
svct2 <-svct[svct<maxv]
if (length(svct2) < 2) break
ucv.l <- mean(svct2)-(crco*sd(svct2))
ucv.u <- mean(svct2)+(crco*sd(svct2))
}
rt2 <- svct
mean(rt); sd(rt)
mean(rt2); sd(rt2) # 最終出力
cv; crco
range(rt); length(rt)
range(rt2); length(rt2)
## ハイブリッド手続き (Hybrid procedure)
crco <- 2.5; crco.nr <- crco # 基準値は勘で
## まず非再帰手続きを行う
mrt <- mean(rt)
sdrt <- sd(rt)
sszrt <- length(rt)
cv <- c(mrt-(crco*sdrt), mrt+(crco*sdrt)); cv.nr <- cv
rt.nrc <- rt[cv[1]<=rt & rt<=cv[2]]
## 修正再帰手続きを行う。この際、非再帰手続きの基準値に1プラスする
crco <- crco.nr+1; crco.mr <- crco
rttmp <- rt[rt<max(rt)]
mrttmp <- mean(rttmp)
sdrttmp <- sd(rttmp)
cv <- c(mrttmp-(crco*sdrttmp), mrttmp+(crco*sdrttmp)); cv.mr <- cv
minv <- min(rt); maxv <- max(rt)
ucv.l <- cv[1]; ucv.u <- cv[2]
svct <- rt
while((minv < ucv.l | maxv > ucv.u)&(length(svct)>3)) {
if (minv < ucv.l) {svct<- svct[svct>minv]}; cv.mr[1] <- ucv.l
if (maxv > ucv.u) {svct<- svct[svct<maxv]}; cv.mr[2] <- ucv.u
cat(ucv.l, ucv.u, crco, "\n") # 最後に使われた基準値とカットオフポイントを返す
if (length(svct) < 2) break
maxv <- max(svct)
minv <- min(svct)
svct2 <-svct[svct<maxv]
if (length(svct2) < 2) break
ucv.l <- mean(svct2)-(crco*sd(svct2))
ucv.u <- mean(svct2)+(crco*sd(svct2))
}
rt.mrc <- svct
## 非再帰手続の平均と修正再帰手続の平均を平均する
mean(c(mean(rt.nrc),mean(rt.mrc))) # 最終の代表値
mean(rt); sd(rt); length(rt)
mean(rt.nrc); sd(rt.nrc); length(rt.nrc); cv.nr; crco.nr
mean(rt.mrc); sd(rt.mrc); length(rt.mrc); cv.mr; crco.mr
range(rt); range(rt.nrc); range(rt.mrc)
より。
結論としては基準変動型のいずれかがお勧めらしい
## 仮想反応時間データ生成
rt <- round(sample(rnorm(10000, mean=300, sd=10)+(300*rexp(10000)), 20))
hist(rt)
win.graph(); plot(rt)
mean(rt); sd(rt)
range(rt)
# 基準変動型手続きに使うデータ
ssz <- c(4,5,6,7,8,9,10,11,12,13,14,15,20,25,30,35,50,100) # 抽出するサンプルの数
nrv <- c(1.458,1.68,1.841,1.961,2.050,2.12,2.173,2.22,2.246,2.274,2.31,2.326,2.391,2.410,2.4305,2.45,2.48,2.50) # 基準変動型非再帰手続き用
mrv <- c(8, 6.20, 5.30, 4.80, 4.475,4.25,4.11,4.00,3.92,3.85,3.80,3.75,3.64,3.595,3.55,3.54,3.51,3.50) # 基準変動型修正再帰手続き用
## 表にしてみる
tb <- rbind(nrv, mrv)
colnames(tb) <- ssz
tb
## 線形補間で予測値を得る関数を生成
nrvfun <- approxfun(ssz, nrv) # 基準変動型非再帰手続き用
mrvfun <- approxfun(ssz, mrv) # 基準変動型修正再帰手続き用
## 基準変動型非再帰手続き (Non-recursive with moving criterion)
mrt <- mean(rt)
sdrt <- sd(rt)
sszrt <- length(rt)
crco <- nrvfun(sszrt) # 線形補間関数による基準値
cv <- c(mrt-(crco*sdrt), mrt+(crco*sdrt))
rt2 <- rt[cv[1]<=rt & rt<=cv[2]]
mean(rt); sd(rt); length(rt)
mean(rt2); sd(rt2); length(rt2) # 最終出力
range(rt); range(rt2); cv; crco
## 基準変動型修正再帰手続き (Modified recursive with moving criterion)
rttmp <- rt[rt<max(rt)]
mrttmp <- mean(rttmp)
sdrttmp <- sd(rttmp)
sszrt <- length(rt)
crco <- mrvfun(sszrt) # 線形補間関数による基準値
cv <- c(mrttmp-(crco*sdrttmp), mrttmp+(crco*sdrttmp))
minv <- min(rt); maxv <- max(rt)
ucv.l <- cv[1]; ucv.u <- cv[2]
svct <- rt
while((minv < ucv.l | maxv > ucv.u)&(length(svct)>3)) {
if (minv < ucv.l) {svct<- svct[svct>minv]}
if (maxv > ucv.u) {svct<- svct[svct<maxv]}
cat(ucv.l, ucv.u, crco, "\n") # 最後に使われた基準値とカットオフポイントを返す
if (length(svct) < 2) break
maxv <- max(svct)
minv <- min(svct)
svct2 <-svct[svct<maxv]
if (length(svct2) < 2) break
ucv.l <- mean(svct2)-(crco*sd(svct2))
ucv.u <- mean(svct2)+(crco*sd(svct2))
}
rt2 <- svct
mean(rt); sd(rt)
mean(rt2); sd(rt2) # 最終出力
cv; crco
range(rt); length(rt)
range(rt2); length(rt2)
## ハイブリッド手続き (Hybrid procedure)
crco <- 2.5; crco.nr <- crco # 基準値は勘で
## まず非再帰手続きを行う
mrt <- mean(rt)
sdrt <- sd(rt)
sszrt <- length(rt)
cv <- c(mrt-(crco*sdrt), mrt+(crco*sdrt)); cv.nr <- cv
rt.nrc <- rt[cv[1]<=rt & rt<=cv[2]]
## 修正再帰手続きを行う。この際、非再帰手続きの基準値に1プラスする
crco <- crco.nr+1; crco.mr <- crco
rttmp <- rt[rt<max(rt)]
mrttmp <- mean(rttmp)
sdrttmp <- sd(rttmp)
cv <- c(mrttmp-(crco*sdrttmp), mrttmp+(crco*sdrttmp)); cv.mr <- cv
minv <- min(rt); maxv <- max(rt)
ucv.l <- cv[1]; ucv.u <- cv[2]
svct <- rt
while((minv < ucv.l | maxv > ucv.u)&(length(svct)>3)) {
if (minv < ucv.l) {svct<- svct[svct>minv]}; cv.mr[1] <- ucv.l
if (maxv > ucv.u) {svct<- svct[svct<maxv]}; cv.mr[2] <- ucv.u
cat(ucv.l, ucv.u, crco, "\n") # 最後に使われた基準値とカットオフポイントを返す
if (length(svct) < 2) break
maxv <- max(svct)
minv <- min(svct)
svct2 <-svct[svct<maxv]
if (length(svct2) < 2) break
ucv.l <- mean(svct2)-(crco*sd(svct2))
ucv.u <- mean(svct2)+(crco*sd(svct2))
}
rt.mrc <- svct
## 非再帰手続の平均と修正再帰手続の平均を平均する
mean(c(mean(rt.nrc),mean(rt.mrc))) # 最終の代表値
mean(rt); sd(rt); length(rt)
mean(rt.nrc); sd(rt.nrc); length(rt.nrc); cv.nr; crco.nr
mean(rt.mrc); sd(rt.mrc); length(rt.mrc); cv.mr; crco.mr
range(rt); range(rt.nrc); range(rt.mrc)
PR
x<-1:10
y<-x # 同じもの
y10<-c(1:9,100)#最後だけ違う
y01<-c(100,2:10)#最初だけ違う
# 基本的に等号演算子で確認できる
x==y
x==y10
x==y01
## 同じなのはいくつかカウント
sum(x==y)
sum(x==y10)
#とりあえず同じかどうか調べる
identical(x,y10)#ベクトルだけではなくオブジェクトもわかる
#yの中にxのどの要素があるか調べる。yの中にあるxの要素を返し、そうでないものはNA。excelでいうとcountifみたいなものか
match(x,y10)
#countif的にマッチするものの数を数える
sum(complete.cases(match(x,y10)))
# どこが同じ調べる。要素の位置を返す
xx <- c(11:20, 1:10)
yy <- c(11:20, 101:110)
which(xx==yy) # 要素の1番目から10番目まで同じ
y<-x # 同じもの
y10<-c(1:9,100)#最後だけ違う
y01<-c(100,2:10)#最初だけ違う
# 基本的に等号演算子で確認できる
x==y
x==y10
x==y01
## 同じなのはいくつかカウント
sum(x==y)
sum(x==y10)
#とりあえず同じかどうか調べる
identical(x,y10)#ベクトルだけではなくオブジェクトもわかる
#yの中にxのどの要素があるか調べる。yの中にあるxの要素を返し、そうでないものはNA。excelでいうとcountifみたいなものか
match(x,y10)
#countif的にマッチするものの数を数える
sum(complete.cases(match(x,y10)))
# どこが同じ調べる。要素の位置を返す
xx <- c(11:20, 1:10)
yy <- c(11:20, 101:110)
which(xx==yy) # 要素の1番目から10番目まで同じ
書き溜めたものをQuick-Rのようにしたいと思ったがこのブログでは上手くいかないな…
階層のあるメニューにできればいいんだけど
Rjpwikiで楽しそうなことが始まっている
psych, sem, carパッケージなんかもあるといいな
階層のあるメニューにできればいいんだけど
Rjpwikiで楽しそうなことが始まっている
psych, sem, carパッケージなんかもあるといいな
SPSSのフリーのクローン、PSPPというのがあるらしい。
SPSSのsavというデータファイルも開くことができる
http://www.gnu.org/software/pspp/
RやOpenOfficeでもそうだけど、この世には限られた自分の時間を
こういうことに使う人がいて、有料のソフトやコンテンツはさらに高機能化が
求められるだろう。
SPSSのsavというデータファイルも開くことができる
http://www.gnu.org/software/pspp/
RやOpenOfficeでもそうだけど、この世には限られた自分の時間を
こういうことに使う人がいて、有料のソフトやコンテンツはさらに高機能化が
求められるだろう。
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) # 対角線
# もっとうまい描きかたがあるだろうな…
