3

包中有一个函数exactci，我想将参数作为矩阵传递给它并取回一个矩阵。事实上，所有参数只能是长度为 1 的向量。我挖掘了源代码并找到了这个片段，即我实际使用的函数（这里修改和减少了参数）：

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])
}

minlike（这将返回 ap 值，用于使用该方法对泊松率进行相等性的双边检验）

我看到我无法传入矩阵并取回矩阵的原因是因为support在内部创建了向量。我将部分剥离dbinom()为以下内容：

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

这会返回相同的向量，f很好，花花公子，甚至更快，但它似乎并没有解决我的问题 - 至少我看不到将其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()返回一个与传递给它的参数相同维度的矩阵吗？更好的是，有没有办法正确矢量化我在这里所做的事情并完全避免隐藏 for 循环（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) 没有 avthis 问题，但由于某种原因它不如exact.binom.minlike.stripped(). 编译 (3) 并没有使它在我的系统上运行得更快。

[[Edit 2]]：在相同的数据上，Pete's newexact.binom.minlike.stripped3()的工作如下：

   user  system elapsed 
  6.468   0.032   6.513

因此，后期的策略，预先计算 的对数阶乘max(n)，是一个主要的节省时间的方法。非常感谢皮特！

4

1 回答 1

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)
}

对于您的示例，这两种方法在我的笔记本电脑上的速度大致相同，尽管根据您的问题和硬件的实际大小，其中一种可能更快。

于 2012-11-18T00:50:01.240 回答