r - 如何在 R 中加速这个简单的函数

Question

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

但是，对于大型数据集，它非常慢。因此，我尝试分析和检查源代码。

大部分时间 (>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近似于无穷大)。

这个问题真的不需要特殊的背景统计知识，但是 link或 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 > 10000，lambda ^ 100 肯定是 NAN。如果我用其他语言编程，也许我会使用 reduce ，但我担心 R 中的 reduce 很慢。

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

您可能可以进一步优化它，但这应该会获得大部分 yield ，除非您可以利用某种内置函数而不是显式求和。

输入长度为 10000 时，在我的机器上需要 0.148 秒，而 test 需要 6.850 秒:

> 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