包中有一个函数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)
,是一个主要的节省时间的方法。非常感谢皮特!