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 回答