myscratchpad

こちらのページを参考。参考というかコピペ。感謝
http://hosho.ees.hokudai.ac.jp/~kubo/log/2010/kubolog20100215.html

plot(1:10) # まず最初のプロット (x 軸と y 軸(左側)が書かれる)
print(par("usr")) # その時のユーザー座標系
par(usr=c(par("usr")[1:2], 100.8, 105.2)) # ユーザー座標系を変更 (x 座標系はそのまま)
points(1:5, 105:101, type = "b", col="red") # 追加の点を変更座標系に赤色で書く (軸は書かれない)
axis(4, col = "red") # 変更座標系の (y) 軸を右側に書く

この記事のとは違う
対数変換してから非再帰法+基準値置換でやる
Ratcliff, R. (1993). Methods for dealing with reaction time outliers. Psychological Bulletin, 114, 510-532.
っぽいかんじ

nn <- 50
rt <- round(c((rnorm(nn, mean=300, sd=20)+(300*rexp(nn)))))

## 関数定義。反応時間のベクトルrt, 基準cop (デフォルト2.5SD) を与える。ifm=Tにすると他の情報も出す。
elrt <- function(rt, cop=2.5, ifm=FALSE) {
rt <- na.omit(as.vector(rt))
rtl <- log(rt)
cop <- cop
meanv <- mean(rtl)
sdv <- sd(rtl)
cv.l <- meanv-(cop*sdv)
cv.u <- meanv+(cop*sdv)
rtl2 <- rtl
tf.l <- rtl<cv.l
tf.u <- rtl>cv.u
rtl2[tf.l] <- cv.l # 外れ値をカットオフ値で置換
rtl2[tf.u] <- cv.u
    nart <- rt
    nart[tf.l] <- NA
    nart[tf.u] <- NA
    rplrt <- rt
    rplrt[tf.l] <- exp(cv.l) # 生rtの外れ値に対しカットオフ値を指数変換して置換
    rplrt[tf.u] <- exp(cv.u)
    logrt <- rtl2
    reslist <- list()
    reslist$logrt <- logrt
    reslist$nart <- nart
    reslist$rplrt <- rplrt
    reslist$outliers.l <- rt[tf.l]
    reslist$outliers.u <- rt[tf.u]
    reslist$cop <- paste(cop, "SD", sep="")
    reslist$cv <- c(cv.l, sum(tf.l), cv.u, sum(tf.u))
    names(reslist$cv) <-  c("lower", "nofl", "upper",  "nofu")
ifelse(ifm==FALSE, return(rtl2), return(reslist))
}
(logrt <- elrt(rt)) # 対数変換・カットオフ置換のみ
(elres <- elrt(rt, ifm=T)) # logrtに加え、外れ値をNAにしたもの、指数変換で置き換えたもの、カットオフ値とその個数

## 基本的に、解析等が終了しその結果をhtmlファイルにまとめておくためのもの
dat <- data.frame(s=paste("s", 1:10), val=rnorm(10))
library(psych)

library(R2HTML)
HTMLStart(outdir=getwd(), file="report",extension="html", echo=FALSE, HTMLframe=TRUE)
HTML.title("My Report", HR=1)

HTML.title("Description of my data", HR=3)
describe(dat)

HTMLhr()

HTML.title("X Y Scatter Plot", HR=2)
plot(dat$val)
HTMLplot()

HTMLStop()


Quick-R: http://www.statmethods.net/interface/output.html
CRAN: http://cran.r-project.org/web/packages/R2HTML/index.html
開発者 (?) : http://www.stat.ucl.ac.be/R2HTML/

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)

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番目まで同じ