如何在R中加速这个简单的功能

Question

我正在尝试使用R中的COMPoissonReg进行Conway-Maxwell-Poisson回归

但是,对于大型数据集来说,它非常慢.因此,我试图剖析并检查源代码.

大部分时间(> 95%)花在一个函数上COMPoissonReg:::computez,相当于:

test <- function (lambda, nu, max=100) 
{
    forans <- matrix(0, ncol = max + 1, nrow = length(lambda))
    for (j in 1:max) {
        temp <- matrix(0, ncol = j, nrow = length(lambda))
        for (i in 1:j) {
            temp[, i] <- lambda/(i^nu)
        }
        for (k in 1:length(lambda)) {
            forans[k, j + 1] <- prod(temp[k, ])
        }
    }
    forans[, 1] <- rep(1, length(lambda))
    ans <- rowSums(forans)
    return(ans)
}

v这里是nu,lambda是向量,max是上限s(这里它设置为100作为近似无穷大).

并不真正需要的特殊背景的问题统计知识,但链接或LINK2是在这里,以防万一.

一个简单的脚本来测试性能,这需要8秒,如果我懒洋洋地cmpfun编译它,它需要4秒.我相信它有可能得到进一步改善.(没有在C中重写,我的目标是约0.05秒,这样我就不必重构迭代调用此函数的包中的代码.)

lambda <- rnorm(10000, 1.5, 0.3)
Rprof(tmp <- tempfile())
sum(log(test(lambda, 1.2)))
Rprof()
summaryRprof(tmp)

更新

我意识到另一个问题:浮点运算限制.做电源系列是危险的,它很快就会溢出,特别是如果我们必须进行矢量化.例如lambda ^ 100,当然> NAN如果lambda> 10000. reduce如果我用其他语言编程,我可能会使用,但我担心R减速很慢.

Answer 1

通过避免循环,您可以比使用的功能更快.例如:

test2<-function(lambda,nu,max=100){
  len<-length(lambda)
  mm<-matrix(rep(lambda,each=max+1),max+1,len)
  mm<-mm^(0:max)
  mm<-mm/factorial(0:max)^nu
  colSums(mm)
}

使用长度为100的lambda,运行速度提高约50倍:

> require(microbenchmark)
> lam<-rnorm(100)
> max(abs(test(lam,1.2)-test2(lam,1.2)))
[1] 4.510281e-16
> microbenchmark(test(lam,1.2),test2(lam,1.2),times=10)
Unit: milliseconds
            expr       min        lq    median        uq       max neval
  test(lam, 1.2) 77.124705 77.422619 78.241945 79.635746 81.260280    10
 test2(lam, 1.2)  1.335716  1.373116  1.401411  1.507765  1.562447    10

您可以更多地优化它,但这应该获得大部分收益,除非您可以利用某种内置函数而不是明确地进行求和.

在输入长度10000时,我的机器需要0.148秒,而以下是6.850秒test:

> lam<-rnorm(10000)
> max(abs(test(lam,1.2)-test2(lam,1.2)))
[1] 3.552714e-15
> system.time(test2(lam,1.2))
   user  system elapsed 
  0.132   0.016   0.148 
> system.time(test(lam,1.2))
   user  system elapsed 
  6.780   0.056   6.850