R 2.13.0のコンパイラをつかってみる (2) Cλ

R 2.13.0のコンパイラをつかってみるその2。類似度指数のC_λを計算してみる。

library(compiler)

# MorisitaのCλ
clambda <- function(c1, c2) {
  if (length(c1) != length(c2)) {
    stop("number of species not matched.")
  }
  n1 <- sum(c1)
  n2 <- sum(c2)
  if (n1 == 0 | n1 == 1 | n2 == 0 | n2 == 1) {
    stop("too few number of individuals.")
  }
  l1 <- sum(c1 * (c1 - 1)) / (n1 * (n1 - 1))
  l2 <- sum(c2 * (c2 - 1)) / (n2 * (n2 - 1))
  cc <- sum(c1 * c2)
  return (2 * cc / ((l1 + l2) * n1 * n2))
}

# Cλ行列 -- 各行が群集
mclambda <- function(mat) {
  n.cmm <- nrow(mat)
  cl <- matrix(0, ncol = n.cmm, nrow = n.cmm)
  for (i in 1:n.cmm) {
    for (j in 1:n.cmm) {
      cl[i, j] <- clambda(mat[i, ], mat[j, ])
    }
  }
  cl
}

## 架空データ
set.seed(123)
# 種数
n.spc <- 160
# 群集の数
n.cmm <- 480
# 個体数が0でない確率
p.nz <- 0.5
# 個体数が0でないところの平均個体数
m.ppl <- 8
# 群集ごと種ごとの個体数
cmm <- matrix(rbinom(n.spc * n.cmm, 1, p.nz) *
                rpois(n.spc * n.cmm, m.ppl),
              ncol = n.spc)


## コンパイル
c.clambda <- cmpfun(clambda)
mclambda2 <- function(mat) {
  n.cmm <- nrow(mat)
  cl <- matrix(0, ncol = n.cmm, nrow = n.cmm)
  for (i in 1:n.cmm) {
    for (j in 1:n.cmm) {
      cl[i, j] <- c.clambda(mat[i, ], mat[j, ])
    }
  }
  cl
}

c.mclambda <- cmpfun(mclambda2)

結果

> system.time(mclambda(cmm))
   ユーザ   システム       経過  
     9.068      0.034      9.043 
> system.time(mclambda2(cmm))
   ユーザ   システム       経過  
     6.744      0.023      6.714 
> system.time(c.mclambda(cmm))
   ユーザ   システム       経過  
     5.871      0.017      5.841 

タグ：R