使Vectorize（）传递变暗或正确地向量化该函数

Question

exactci包中有一个函数，我想将参数作为矩阵传回并获取矩阵。因为它是，所有参数只能是长度为1的向量I挖成源，并发现这一块，我实际使用的功能（在这里与改性参数和还原的）：

exact.binom.minlike <- function(d1, d2, e1, e2){ 
    x   <- round(d1) 
    n   <- x + round(d2) 
    p   <- e1/(e1 + e2) 

    support  <- 0:n 
    f   <- dbinom(support, n, p) 
    d   <- f[support == x] 

    sum(f[f <= d * relErr]) 
} 

（这个返回p值对于使用minlike方法的泊松比相等的双面测试）

我看到我无法传入矩阵并取回矩阵的原因是由于内部创建了矢量support。我剥离下来的dbinom()部分如下：

f   <- exp(lfactorial(n) - 
        (lfactorial(support) + lfactorial(n - support)) + 
        support * log(p) + 
        (n - support) * log(1 - p) 
        ) 

这还给了相同的向量，f，罚款和花花公子，甚至有点快，但它似乎并没有解决我的problem-至少我不没有看到使用support作为一种载体。支持的长度会根据d1+d2的不同而有所不同，所以我一直在坚持一次比较。我已经能够做的最好的是坚持内Vectorize()整个事情，这需要矩阵就好作为参数，但返回一个向量，而不是一个矩阵：

exact.binom.minlike.stripped <- Vectorize(compiler:::cmpfun(function(d1, d2, e1, e2, relErr = 1 + 10^(-7)){ 
    x   <- round(d1) 
    n   <- x + round(d2) 
    p   <- e1/(e1 + e2) 

    support  <- 0:n 

    # where dbinom() is the prob mass function: 
    # n choose k * p^k * (1 - p)^(n - k) # log it to strip down, then exp it 
    f   <- exp(lfactorial(n) - 
         (lfactorial(support) + lfactorial(n - support)) + 
         support * log(p) + 
         (n - support) * log(1 - p) 
         ) 
    #f   <- dbinom(support,n,p) 
    d   <- f[support == x] 

    sum(f[f <= d * relErr]) 
})) 

下面是一个例子：

set.seed(1) 
d1 <- matrix(rpois(36,lambda = 100), 6) 
d2 <- matrix(rpois(36,lambda = 150), 6) 
e1 <- matrix(rpois(36,lambda = 10000), 6) 
e2 <- matrix(rpois(36,lambda = 25000), 6) 

该输出是一个长度为36而不是6x6矩阵的向量。所有四个输入为6x6的矩阵：

(p.vals <- exact.binom.minlike.stripped(d1, d2, e1, e2)) 
[1] 1.935277e-04 9.680425e-08 1.508232e-08 1.227176e-04 1.656111e-02 
[6] 2.310620e-04 2.871150e-05 4.024025e-06 4.804943e-05 1.619866e-02 
[11] 3.610596e-02 1.101247e-04 5.153746e-04 1.350891e-04 8.663191e-06 
[16] 1.384378e-05 2.681715e-06 4.556092e-08 2.270317e-04 2.040001e-04 
[21] 3.330344e-01 4.775055e-05 2.588667e-07 5.647732e-04 1.615861e-03 
[26] 2.438345e-03 2.524692e-04 3.398664e-05 2.001322e-05 4.361194e-03 
[31] 3.909116e-05 1.697943e-03 8.543677e-07 2.992653e-05 2.617216e-04 
[36] 3.106748e-03 

我收集我可以添加dim() s到这一点，并使其恢复到矩阵：

dim(p.vals) <- dim(d1) 

但似乎退而求其次。我可以让Vectorize()返回一个与传递给它的参数尺寸相同的矩阵吗？更好的是，是否有一种方法可以将我在此处所做的操作正确地矢量化并避免隐藏循环（Vectorize()使用mapply()）？

[[编辑]]感谢皮特的伟大建议。这里有一个比较使用的数据在尺寸上更接近什么，我实际上做：

set.seed(1) 
N <-110 
d1 <- matrix(rpois(N^2,lambda = 1000), N) 
d2 <- matrix(rpois(N^2,lambda = 1500), N) 
e1 <- matrix(rpois(N^2,lambda = 10000), N) 
e2 <- matrix(rpois(N^2,lambda = 25000), N) 

system.time(exact.binom.minlike.stripped.2(d1, d2, e1, e2)) 
    user system elapsed 
16.353 1.112 17.635 
system.time(exact.binom.minlike.stripped.3(d1, d2, e1, e2)) 
    user system elapsed 
14.685 0.016 14.715 
system.time({ 
     (p.vals <- exact.binom.minlike.stripped(d1, d2, e1, e2)) 
     (dim(p.vals) <- dim(d1)) 
    }) 
    user system elapsed 
12.541 0.040 12.604 

我在这看着我的系统监控内存使用，只有exact.binom.minlike.stripped.2()是一个记忆体猪。我发现如果我在真实数据上使用这个数据，那么max(n)可能会增大10-20倍，那么我的电脑就会窒息。 （3）不存在这个问题，但由于某种原因，它不如exact.binom.minlike.stripped()那么快。编译（3）并没有让我的系统运行得更快。

[编辑2]：在相同的数据，皮特的新exact.binom.minlike.stripped3()做的工作：

user system elapsed 
    6.468 0.032 6.513 

因此，后来stretegy，预先计算的max(n)日志阶乘，是一个重要的时间-saver。非常感谢皮特！

Answer 1

我可以想出两个想要像这样矢量化的功能的原因：便利性或性能。

为方便起见，以下内容应该有效，但我怀疑如果max(n)非常大，那么所有的内存分配将抵消dbinom调用向量化的任何收益。

exact.binom.minlike.stripped.2 <- function(d1, d2, e1, e2, relErr = 1 + 1e-7) { 

    x <- round(d1) 
    n <- x + round(d2) 
    p <- e1/(e1 + e2) 

    # `binom` is already vectorised. 
    d <- dbinom(x, n, p) 

    # rearrange inputs to `dbinom` so that it works with `outer`. 
    dbinom.rearrange <- function(n, x, p) dbinom(x, n, p) 
    support <- 0:max(n) 
    f <- outer(n, support, dbinom.rearrange, p=p) 

    # repeat `d` enough times to conform with `f`. 
    d <- array(d, dim(f)) 
    f[f > d * relErr] <- 0 

    # extract the required sums. 
    apply(f, c(1,2), sum) 
} 

---

或者，可能是更明智的方式来做到这一点：尽可能将去利用自然矢量化，并限制Vectorize的“不自然”的一部分。这仍然需要在最后修复尺寸。

vector.f <- Vectorize(function(d, n, p, ftable) { 

    x <- 0:n 
    f <- exp(ftable[n+1] - (ftable[x+1] + ftable[n-x+1]) + x*log(p) + (n-x)*log(1-p)) 
    sum(f[f <= d]) 

}, c('d', 'n', 'p')) 

exact.binom.minlike.stripped.3 <- function(d1, d2, e1, e2, relErr = 1 + 1e-7) { 

    x <- round(d1) 
    n <- x + round(d2) 
    p <- e1/(e1 + e2) 

    # `binom` is already vectorised. 
    d <- dbinom(x, n, p) 

    # precompute factorials 
    ftable <- lfactorial(0:max(n)) 

    f <- vector.f(d * relErr, n, p, ftable) 
    dim(f) <- dim(d1) 

    return(f) 
} 

这两种出来关于我的笔记本电脑为您的示例相同的速度，但一个或另一个可能更快取决于你的问题，你的硬件的实际大小。