是否有可能加快我创建相关矩阵的功能？

Question

我已经编写了以下函数来估计使用所谓的CramérV的多项式变量的成对相关性.我vcd为此目的使用了包,但据我所知,没有现有的函数会从矩阵或data.frame类似的函数创建V的对称相关矩阵cor.

功能是:

require(vcd)
get.V<-function(y){
  col.y<-ncol(y)
  V<-matrix(ncol=col.y,nrow=col.y)
  for(i in 1:col.y){
    for(j in 1:col.y){
      V[i,j]<-assocstats(table(y[,i],y[,j]))$cramer
    }
  }
  return(V)
}

但是,对于大量变量,它变得相对较慢.

no.var<-5
y<-matrix(ncol=no.var,sample(1:5,100*no.var,TRUE))
get.V(y)

随着您的增加,no.var计算时间可能会爆炸.由于我需要将其应用于data.frame100和更高的长度,我的问题是,是否有可能通过更优雅的编程来"加速"我的功能.谢谢.

Answer 1

除了减少执行的测试次数,或者以其他方式优化整个功能的运行,我们也可以 assocstats加快速度.我们首先建立一个测试用例,以确保我们不会意外地创建一个不正确的更快的函数.

x <- vcd::Arthritis$Improved
y <- vcd::Arthritis$Treatment
correct <- vcd::assocstats(table(x, y))$cramer
correct

## [1] 0.3942

is_ok <- function(x) stopifnot(all.equal(x, correct))

我们首先要制作一个assocstats非常接近原版的版本.

cramer1 <- function (x, y) {
  mat <- table(x, y)

  tab <- summary(MASS::loglm(~1 + 2, mat))$tests

  phi <- sqrt(tab[2, 1] / sum(mat))
  cont <- sqrt(phi ^ 2 / (1 + phi ^ 2))

  sqrt(phi ^ 2 / min(dim(mat) - 1))
}
is_ok(cramer1(x, y))

这里最慢的操作将是loglm,所以在我们尝试更快地进行之前,值得寻找替代方法.一个小的谷歌搜索找到一个有用的博客文章.我们也试试:

cramer2 <- function(x, y) {
  chi <- chisq.test(x, y, correct=FALSE)$statistic[[1]]

  ulength_x <- length(unique(x))
  ulength_y <- length(unique(y))

  sqrt(chi / (length(x) * (min(ulength_x, ulength_y) - 1)))
}
is_ok(cramer2(x, y))

性能如何叠加:

library(microbenchmark)

microbenchmark(
  cramer1(x, y),
  cramer2(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1080.0 1149.3 1182.0 1222.1 2598   100
##  cramer2(x, y)  800.7  850.6  881.9  934.6 1866   100

cramer2()是比较快的.chisq.test()可能是瓶颈,所以让我们看看我们是否可以通过做更少的事情来更快地完成这个功能: chisq.test()除了计算测试统计数据之外还有很多,所以我们可能会让它更快.几分钟的细心工作将功能降低到:

chisq_test <- function (x, y) {
  O <- table(x, y)
  n <- sum(O)

  E <- outer(rowSums(O), colSums(O), "*")/n

  sum((abs(O - E))^2 / E)
}

然后,我们可以创建一个新的cramer3()使用chisq.test().

cramer3 <- function(x, y) {
  chi <- chisq_test(x, y)

  ulength_x <- length(unique(x))
  ulength_y <- length(unique(y))

  sqrt(chi / (length(x) * (min(ulength_x, ulength_y) - 1)))
}
is_ok(cramer3(x, y))
microbenchmark(
  cramer1(x, y),
  cramer2(x, y),
  cramer3(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1088.6 1138.9 1169.6 1221.5 2534   100
##  cramer2(x, y)  796.1  840.6  865.0  906.6 1893   100
##  cramer3(x, y)  334.6  358.7  373.5  390.4 1409   100

而现在,我们有我们自己的简单版本chisq.test(),我们可以利用的结果伊克出多一点的速度table()找出独特的元素的数量x和y:

cramer4 <- function(x, y) {
  O <- table(x, y)
  n <- length(x)
  E <- outer(rowSums(O), colSums(O), "*")/n

  chi <- sum((abs(O - E))^2 / E)
  sqrt(chi / (length(x) * (min(dim(O)) - 1)))
}
is_ok(cramer4(x, y))
microbenchmark(
  cramer1(x, y),
  cramer2(x, y),
  cramer3(x, y),
  cramer4(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1097.6 1145.8 1183.3 1233.3 2318   100
##  cramer2(x, y)  800.7  840.5  860.7  895.5 2079   100
##  cramer3(x, y)  334.4  353.1  365.7  384.1 1654   100
##  cramer4(x, y)  248.0  263.3  273.2  283.5 1342   100

不错 - 我们使用R代码的速度提高了4倍.从这里开始,您可以通过以下方式获得更快的速度:

• 使用tcrossprod()而不是outer()
• table()为这个特殊(2d)案例制作更快的版本
• 使用Rcpp从表格数据计算测试统计量