R で bootstrap method を実装してみよう

高速 tapply(X, INDEX, var) もどき関数 を作ったときに、簡単な実装を考えていたのでまとめてみます。

xxx.test <- function(d, g) {

  n      <- as.numeric(table(g))
  index  <- as.numeric(factor(g))
  mindex <- NULL
  for(i in 1:length(unique(index))) {
    mindex <- cbind(mindex, index-(i-1))
  }
  mindex[mindex!=1] <- 0

  p   <- length(unique(g))
  df  <- length(g) - p
  s   <- t(mindex) %*% d^2 - (t(mindex) %*% d)^2 / n
  v   <- t(rep(1, p)) %*% s/df

  T <- (c(1, -1) %*% (t(mindex) %*% d / n)) / (sqrt(v * sum(1/n)))

  return(T)

}

xxx.eval <- function(T, t0, tail) {

  if (tail=="two.sided") {
    q <- floor(abs(T)/as.vector(abs(t0)))
    q[q>1] <- 1
  }
  else if (tail=="greater" | tail=="less") {
    q <- floor(T/as.vector(t0))
    q[q<0] <- 0
    q[q>1] <- 1
  }

  pvalue <- sum(q)/r
  if (tail=="less"    & t0 > 0) pvalue <- 1 - pvalue
  else if (tail=="greater" & t0 < 0) pvalue <- 1 - pvalue

  return(pvalue)

}

bootmatrix <- function(d0, g, r, tail="two.sided") {

  n <- as.numeric(table(g))
  d <- do.call("rbind", 
         mapply(function(n) matrix(sample(d0, r * n, replace=T),
         nrow=n, ncol=r), n, SIMPLIFY = FALSE)
       )

  T  <- xxx.test( d, g)
  t0 <- xxx.test(d0, g)

  pvalue <- xxx.eval(T, t0, tail)

  return(list(T=T, d0=d0, t0=t0, pvalue=pvalue))

}

n  <- c(10, 10)
d0 <- c(rnorm(n[1], 0, 1), rnorm(n[2], 0, 1))
g  <- rep(1:length(n), n)

r  <- 100000

bm100K <- bootmatrix(d0, g, r)
tt     <- t.test(d0[1:10], d0[11:20], var.equal = TRUE)
bm100K$t0; bm100K$pvalue;
tt$statistic; tt$p.value;

統計量を求める関数 xxx.test が 行列演算 で書ける場合は高速化が期待できます。
評価関数 xxx.eval と標本抽出部分を、多変量化出来れば・・・

計算時間

Rprof()
bm1M <- bootmatrix(d0, g, 1000000)
Rprof(NULL)
print(a<-summaryRprof()$sampling.time)

$sampling.time
[1] 3.28

100万サンプル

関連記事

高速 tapply(X, INDEX, var) もどき関数