R：機能を選択する、効率的

Question

私は自分のRコードをスピードアップする必要があります。私のボトルネックはchoose関数を使用する必要がある関数です。

P_ni <- function(Pn,Pi,eta1,eta2,p,d=NA) 
{ 
if(is.na(d)) d <- 1-p 
if(Pn==Pi) output <- p^Pn 
else 
{ 
    if(Pi==1)seq1 <- seq_len(Pn-1) 
    if(Pi>1)seq1 <- seq_len(Pn-1)[-seq_len(Pi-1)] 
    output <- sum(choose((Pn-Pi-1),c(seq1-Pi))*choose(Pn,seq1)* 
    (eta1/(eta1+eta2))^c(seq1-Pi)* 
    (eta2/(eta1+eta2))^c(Pn-seq1)*p^seq1*d^c(Pn-seq1) 
) 
} 
return(output) 
} 

この関数は、異なるPnとPiで数回呼び出す必要があります。ここでの問題は、PnとPiだけが単一の数値をとり、ベクトルでは動作できないということです。これはchoose（） - 関数によって発生します。

私は現時点でforループでこれを行いますが、完全に動作しますが、遅いです。

for(i in 1:nrow(n_k_matrix_p)) 
{ 
    n_k_matrix_p[i,4] <- P_ni(n_k_matrix_p[i,1],n_k_matrix_p[i,2],eta1,eta2,p) 
} 

それが再現可能にするために：

eta1 <- 10 
eta2 <- 5 
p <- 0.4 
n_k_matrix <- expand.grid(c(1:20),c(1:20)) 
n_k_matrix <- n_k_matrix[n_k_matrix[,1] >=n_k_matrix[,2],] 
n_k_matrix <- n_k_matrix[order(n_k_matrix[,1]),] 

n_k_matrixがPnとパイのための私の番号が含まれているため、ループ ザ・は次のようになります。 残念ながら、ループは適用を使用するよりもまだ高速です。 誰かが物事をスピードアップする方法を知っていますか？

Answer 1

計算を再集計または事前計算することができます。

P_ni2 <- function(n, eta1, eta2, p, d = 1 - p) { 

    res <- matrix(0, n, n) 
    diag(res) <- p^seq_len(n) 

    C1 <- eta1/eta2 * p/d 
    C2 <- eta2/(eta1 + eta2) * d 
    C3 <- eta1/(eta1 + eta2) 
    C2_n <- C2^seq_len(n) 
    C3_n <- C3^seq_len(n) 
    precomputed <- outer(0:n, 0:n, choose) 

    for (j in seq_len(n)) { 
    for (i in seq_len(j - 1)) { 
     seq1 <- seq(i, j - 1) 
     res[i, j] <- sum(
     precomputed[j-i, seq1-i+1] * precomputed[j+1, seq1+1] * C1^seq1 
    ) * C2_n[j]/C3_n[i] 
    } 
    } 

    res 
} 

VERIF：私は最初の二乗行列の上三角部分に結果を格納

> system.time({ 
+ n_k_matrix[[3]] <- sapply(1:nrow(n_k_matrix), function(i) { 
+  P_ni(n_k_matrix[i,1], n_k_matrix[i,2], eta1, eta2, p) 
+ }) 
+ }) 
utilisateur  système  écoulé 
    11.799  0.000  11.797 

> system.time({ 
+ test <- P_ni2(400, eta1, eta2, p) 
+ n_k_matrix[[4]] <- test[as.matrix(n_k_matrix[, 2:1])] 
+ }) 
utilisateur  système  écoulé 
     2.328  0.003  2.341 


> all.equal(n_k_matrix[[3]], n_k_matrix[[4]]) 
[1] TRUE 

注意。次に、あなたのデータフレームフォーマットで変換します（あなたは行列を呼び出します）。

このソリューションはn = 400の5倍高速です。私はあなたがRcppのダブルループ（のみ）を記録することでそれを改善できると思います。