グループ変数が沢山ある場合の統計量の計算

グループ変数が複数セットあり、測定値が 1 セットある様なデータに対して、標本平均、総和、標本分散、偏差平方和をたぶん高速に計算できる関数を作りました。
以前作った高速 tapply(X, INDEX, var) もどき関数 (グループ変数が1セットで、測定値が複数セットあるようなデータ用) に似ている関数です。

MultipleGroupStat <- function(x, g, type) {

  m <- ncol(g)
  n <- nrow(g)
  uni <- unique(factor(g))
  nuni <- length(uni)

  c <- NULL
  for(i in 1:nuni) {
    c <- cbind(c, g/i)
  }
  c[c!=1] <- 0
  p <- rep(seq(1, 1+m*(nuni-1), by = m), m) + rep(1:m-1, each=nuni)
  c <- c[,p]
  r <- NULL
  
  # Mean
  if (type == 1 || type == "mean") {
    r <- matrix(t(c) %*% x / colSums(c), nrow=nuni)
  }
  # Summation
  else if (type == 2 || type == "sum") {
    r <- matrix(t(c) %*% x, nrow=nuni)
  }
  # Variance (sample variance)
  else if (type == 3 || type == "var") {
    r <- matrix((t(c) %*% x^2 - ((t(c) %*% x)^2 / colSums(c))) /
           (colSums(c) - 1), nrow=nuni)
  }
  # Sum of squared deviation
  else if (type == 4 || type == "ss") {
    r <- matrix(t(c) %*% x^2 - ((t(c) %*% x)^2 / colSums(c)), nrow=nuni)
  }

  return(r)

}

適当にデータを生成し、

set.seed(392784)

n    <- 20     # データの数
m    <- 50000   # グループ変数セットの数
nuni <- 3      # グループの数 (例では 3 水準)

g <- matrix(sample(seq(1, nuni), n*m, replace=TRUE), ncol=m)
x <- matrix(rnorm(n), ncol=1)

apply と by を駆使した場合と比べてみます。

Rprof()
applyby <- apply(g, 2, function(g, x) by(x, g, mean), x)
Rprof(NULL)
# apply() & by()
summaryRprof()$sampling.time

Rprof()
MGS <- MultipleGroupStat(x, g, "mean")
Rprof(NULL)
# MultipleGroupStat()
summaryRprof()$sampling.time

require(compiler)
MultipleGroupStatc <- cmpfun(MultipleGroupStat)
Rprof()
MGS <- MultipleGroupStatc(x, g, "mean")
Rprof(NULL)
# MultipleGroupStatc() (compiler package)
summaryRprof()$sampling.time

> # apply() & by()
> summaryRprof()$sampling.time (2 回分)
[1] 60.6
[1] 59.44

# MultipleGroupStat()
> summaryRprof()$sampling.time (2 回分)
[1] 3.12
[1] 3.16

# MultipleGroupStatc() (compiler package)
> summaryRprof()$sampling.time (2 回分)
[1] 2.82
[1] 2.86

今流行の compiler パッケージの cmpfun 関数を使ってみたところ、若干スピードが上がりました。

履歴

2011/05/16
高速 tapply(X, INDEX, var) もどき関数でご指摘いただいた誤りを同様に修正