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であるかによって異なりますので、私は一度に1つずつ比較することに固執しています。私ができる最善のことは、すべてを内部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)

この出力は、6x6行列ではなく、長さ36のベクトルです。4つの入力はすべて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

これにsを追加dim()して、マトリックスに戻すことができます。

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

しかし、それは2番目に良いようです。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）はこの問題を回避しませんが、何らかの理由で。ほど高速ではありませんexact.binom.minlike.stripped()。（3）をコンパイルしても、私のシステムではそれ以上速く実行されませんでした。

[[編集2]]：同じデータで、ピートの新しいexact.binom.minlike.stripped3()ものは次の場所で仕事をします：

   user  system elapsed 
  6.468   0.032   6.513

したがって、の対数階乗を事前に計算する後の戦略max(n)は、大幅な時間の節約になります。ピートに感謝します！

4

1 に答える 1

1

このようなベクトル化された関数が必要な理由は2つ考えられます。利便性とパフォーマンスです。

以下は便宜上機能するはずですが、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 に答える