読者です 読者をやめる 読者になる 読者になる

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

R

関数を作ってみました。

HS.tapply <- function(x, index, type) {

  # create index indicate matrix
  index <- as.numeric(factor(index))
  mindex <- NULL
  for(i in 1:length(unique(index))) {
    mindex <- cbind(mindex, index-(i-1))
  }
  if (type == 5 || type == "max" || type == 6 || type == "min") {
    mindex[mindex!=1] <- NA
  } else {
    mindex[mindex!=1] <- 0
  }

  r <- NULL
  
  # Mean
  if (type == 1 || type == "mean") {
    r <- t(mindex) %*% x / as.numeric(table(index))
  }
  # Summation
  else if (type == 2 || type == "sum") {
    r <- t(mindex) %*% x
  }
  # Sum of squared deviation
  else if (type == 3 || type == "ss") {
    r <- t(mindex) %*% x^2 - (t(mindex) %*% x)^2 / as.numeric(table(index))
  }
  # Variance (sample variance)
  else if (type == 4 || type == "var") {
    r <- (t(mindex) %*% x^2 - (t(mindex) %*% x)^2 /
         as.numeric(table(index))) / (as.numeric(table(index))-1)
  }
  # Maximum
  else if (type == 5 || type == "max") {
    for(i in 1:length(unique(index))) {
      r <- rbind(r, x[cbind(
             max.col(ifelse(is.na(t(mindex[,i]*x)), -Inf, t(mindex[,i]*x))),
             1:(ncol(x))
           )])
    }
  }
  # Minimum
  else if (type == 6 || type == "min") {
    for(i in 1:length(unique(index))) {
      r <- rbind(r, x[cbind(
             max.col(ifelse(is.na(t(mindex[,i]*x)), -Inf, -t(mindex[,i]*x))),
             1:(ncol(x))
           )])
    }
  }

  return(r)
}

データが行列状になっていても計算してくれます(グループ変数が等しく、異なる測定値がたくさんあるようなデータ)。
データが多いとメモリが足りなくなるので、適当に分割して計算します。

n <- c(5, 4, 7)
g <- rep(1:length(n), n)
gs <- c("a", "a", "a", "a", "a", "b", "b", "b", "b", "c", "c", "c", "c", "c", "c", "c")
rep <- 20000
x <- matrix(rnorm(sum(n)*rep), ncol=rep, nrow=sum(n))

Rprof()
varA <- apply(x, 2, function(x, g) tapply(x, g, var), g)
Rprof(NULL)
print(a<-summaryRprof()$sampling.time)

Rprof()
varF <- HS.tapply(x, g, "var")
#varF <- HS.tapply(x, gs, "var")  # グループ変数は文字列等でもOK
Rprof(NULL)
print(f<-summaryRprof()$sampling.time)

max(varA-varF)
varA[,1:3]
varF[,1:3]
f/a


あまりの早さにびっくり。

> print(a<-summaryRprof()$sampling.time)
[1] 12.02
> print(f<-summaryRprof()$sampling.time)
[1] 0.02
> max(varA-varF)
[1] 1.776357e-15
> varA[,1:3]
        [,1]      [,2]      [,3]
1 0.70445164 2.2072017 0.6809683
2 0.09444263 0.8143882 0.4464633
3 0.28558027 0.9014923 0.9471029
> varF[,1:3]
           [,1]      [,2]      [,3]
[1,] 0.70445164 2.2072017 0.6809683
[2,] 0.09444263 0.8143882 0.4464633
[3,] 0.28558027 0.9014923 0.9471029
> f/a
[1] 0.001663894

条件演算子

R の条件演算子には &, |, &&, || などがありますが、前者二つはベクトル化されている演算子であり、後者二つはベクトル化されていないという違いがあります。
したがって、データフレームなどに subset 関数をつかって条件を満たす症例だけを抽出したい場合などは両者を混同しないように注意が必要です。

df<-data.frame(x=runif(10),y=rnorm(10))
subset(df, y<0 & x>0.5)   # (a) & を使う
subset(df, y<0 && x>0.5)  # (b) && を使う

(a) は正しく二つの条件を同時に満たす症例を抽出できますが、(b) はそうなりません。

履歴

2009/11/19
関数の定義をちょこっと修正

2009/11/27
Max, Minを追加

2011/05/16
ご指摘いただいた誤りを修正

2011/07/03
条件演算子についての記述を追加