myscratchpad

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. を参考

source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1270317393")

ssz <- c(4,5,6,7,8,9,10,11,12,13,14,15,20,25,30,35,50,100)
nos <- 10000
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)
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)

cop <- 2.5
system.time(rawres <- rawfun(ssz, nos, cop)) # 生データ平均値
system.time(medres <- medfun(ssz, nos, cop)) # 生データ中央値

system.time(lores <- lofun(ssz, nos, cop)) # 全体対数変換のみ
system.time(rores <- rofun(ssz, nos, cop)) # 全体対数変換、全体カットオフ
system.time(nrres.lo <- nrcs.lo(ssz, nos, cop)) # 全体対数変換、個別カットオフ
system.time(nrres <- nrcs(ssz, nos, cop)) # 非再帰法
cop <- 3.5
system.time(mrres <- mrcs(ssz, nos, cop)) # 修正再帰法
hybres <- (nrres+mrres)/2 # ハイブリッド法
cop <- nrv
system.time(nrreswmc <- nrcs(ssz, nos, cop)) # 非再帰基準変動
cop <- mrv
system.time(mrreswmc <- mrcs(ssz, nos, cop)) # 修正再帰基準変動

win.graph()
ylmt <- c(500, 650)
plot(ssz, rawres, type="l", xaxt="n", xlim=c(5,100), ylim=ylmt, xlab="Sample Size", yaxt="n", pch=1, col=1, lty=1) # 生データ平均値
axis(1, at=ssz, labels=as.character(ssz))
axis(2, ylim=ylmt)
#axis(2, at=seq(from=ylmt[1], to=ylmt[2], 5), labels=as.character(seq(from=ylmt[1], to=ylmt[2], 5)))
points(ssz, medres, type="l", col=1, pch=1, lty=3) # 生データ中央値
points(ssz, nrres, type="b", col=2, pch=2) # 非再帰
points(ssz, mrres, type="b", col=3, pch=3) # 修正再帰
points(ssz, hybres, type="b", col=4, pch=4) # ハイブリッド
points(ssz, nrreswmc, type="b", col=5, pch=5) # 非再帰基準変動
points(ssz, mrreswmc, type="b", col=6, pch=6) # 修正再帰基準変動
legend("topleft", lty=c(1,3,2), col=c(1,1,2), pch=c(NA, NA, 2), c("rawmean", "rawmedian", "non-recursive"), box.lty=0)
legend("top", lty=c(3:6), col=c(3:6), pch=c(3:6), c("modified recursive", "hybrid", "n.r.moving criterion", "m.r.moving criterion"), box.lty=0)

ylmt <- log(ylmt)
par(usr=c(par("usr")[1:2], ylmt[1], ylmt[2]))
axis(4, ylim=ylmt)
points(ssz, lores, type="b", col=1, ylim=ylmt, lty=2, pch=1)
points(ssz, rores, type="b", col=2, ylim=ylmt, lty=2, pch=2) # 全体カットオフ対数変換
points(ssz, nrres.lo, type="b", col=3, ylim=ylmt, lty=2, pch=3) # 修正再帰法対数変換
legend("topright", lty=c(2,2,2), col=1:3, pch=1:3, c("all.log", "log grand", "log individual"), box.lty=0)
#対数ってこう書いていいんかなー

## 出力pdf

## 表にする
tb1 <- rbind(rawres, medres, nrres, mrres, hybres, nrreswmc, mrreswmc)
tb2 <- rbind(rawres, explo=exp(lores), exprores=exp(rores), expnreslo=exp(nrres.lo), lores, rores, nrres.lo)
colnames(tb1)  <- ssz
colnames(tb2) <- ssz
options(digits=4)
tb1
tb2