2

p1 != p2 != ... != pn となる一連の n 個の独立したベルヌーイ試行の成功確率 (p1 から pn) があるとします。各トライアルに一意の名前を付けます。

    p <- c(0.5, 0.12, 0.7, 0.8, .02)
    a <- c("A","B","C","D","E")

ポアソン二項分布関数を使用して cdf、pmf などを見つけることができることは、スタック交換 (たとえば、こことここ) を検索することでわかっています。

私が興味を持っているのは、成功と失敗のあらゆる可能な組み合わせの正確な確率です。(例: 確率木を描いた場合、各枝の最後にある確率を知りたい。)

    all <- prod(p)
    all
    [1] 0.000672
    o1 <- (0.5 * (1-0.12) * 0.7 * 0.8 * .02)
    o1
    [1] 0.004928
    o2 <- (0.5 * 0.12 * (1-0.7) * 0.8 * .02)
    o2
    [1] 0.000288

...成功/失敗のすべての 2^5 の可能な組み合わせ。

Rでこれを行う効率的な方法は何ですか?

私の実際のデータ セットの場合、試行回数は 19 です。したがって、確率ツリーの合計 2^19 パスについて話していることになります。

4

2 に答える 2

1

この計算を高速化するための鍵は、対数確率空間で実行して、ツリーの各分岐の積が行列乗算の内和として計算できる合計になるようにすることです。このようにして、すべての分岐をベクトル化された方法でまとめて計算できます。

まず、すべてのブランチの列挙を作成します。このために、パッケージのintToBin関数を使用します。R.utils

library(R.utils)
enum.branches <- unlist(strsplit(intToBin(seq_len(2^n)-1),split=""))

ここnで、 はベルヌーイ変数の数です。あなたの例ではn=5：

matrix(enum.branches, nrow=n)
##     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
##[1,] "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"   "0"   "0"   "0"   "0"   "0"   "0"   "1"  
##[2,] "0"  "0"  "0"  "0"  "0"  "0"  "0"  "0"  "1"  "1"   "1"   "1"   "1"   "1"   "1"   "1"   "0"  
##[3,] "0"  "0"  "0"  "0"  "1"  "1"  "1"  "1"  "0"  "0"   "0"   "0"   "1"   "1"   "1"   "1"   "0"  
##[4,] "0"  "0"  "1"  "1"  "0"  "0"  "1"  "1"  "0"  "0"   "1"   "1"   "0"   "0"   "1"   "1"   "0"  
##[5,] "0"  "1"  "0"  "1"  "0"  "1"  "0"  "1"  "0"  "1"   "0"   "1"   "0"   "1"   "0"   "1"   "0"  
##     [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
##[1,] "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"  
##[2,] "0"   "0"   "0"   "0"   "0"   "0"   "0"   "1"   "1"   "1"   "1"   "1"   "1"   "1"   "1"  
##[3,] "0"   "0"   "0"   "1"   "1"   "1"   "1"   "0"   "0"   "0"   "0"   "1"   "1"   "1"   "1"  
##[4,] "0"   "1"   "1"   "0"   "0"   "1"   "1"   "0"   "0"   "1"   "1"   "0"   "0"   "1"   "1"  
##[5,] "1"   "0"   "1"   "0"   "1"   "0"   "1"   "0"   "1"   "0"   "1"   "0"   "1"   "0"   "1"

各列が確率ツリーの分岐からの結果である行列になります。

enum.branches次に、それを使用して、値がlog(p)ifenum.branches=="1"およびother である場所と同じサイズの対数確率の行列を作成しlog(1-p)ます。を使用したデータの場合p <- c(0.5, 0.12, 0.7, 0.8, .02)、これは次のとおりです。

logp <- matrix(ifelse(enum.branches == "1", rep(log(p), 2^n), rep(log(1-p), 2^n)), nrow=n)

次に、対数確率を合計し、指数をとって確率の積を取得します。

result <- exp(rep(1,n) %*% logp)
##         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]   [,10]
##[1,] 0.025872 0.000528 0.103488 0.002112 0.060368 0.001232 0.241472 0.004928 0.003528 7.2e-05
        [,11]    [,12]    [,13]    [,14]    [,15]    [,16]    [,17]    [,18]    [,19]    [,20]
##[1,] 0.014112 0.000288 0.008232 0.000168 0.032928 0.000672 0.025872 0.000528 0.103488 0.002112
        [,21]    [,22]    [,23]    [,24]    [,25]   [,26]    [,27]    [,28]    [,29]    [,30]
##[1,] 0.060368 0.001232 0.241472 0.004928 0.003528 7.2e-05 0.014112 0.000288 0.008232 0.000168
        [,31]    [,32]
##[1,] 0.032928 0.000672

はresult、 の枝の番号付けと同じ順序になりenum.branchesます。

計算を関数にカプセル化できます。

enum.prob.product <- function(n, p) {
  enum.branches <- unlist(strsplit(intToBin(seq_len(2^n)-1),split=""))
  exp(rep(1,n) %*% matrix(ifelse(enum.branches == "1", rep(log(p), 2^n), rep(log(1-p), 2^n)), nrow=n))
}

19独立したベルヌーイ変数を使用してこれをタイミングします。

n <- 19
p <- runif(n)
system.time(enum.prob.product(n,p))
##   user  system elapsed 
## 24.064   1.470  26.082

これは私の 2 GHz MacBook (2009 年頃) にあります。計算自体は非常に高速であることに注意してください。unlist時間の大半を占めているのは、確率ツリーの枝の列挙です (その中にあると思います)。それを行うための別のアプローチに関するコミュニティからの提案は大歓迎です。

于 2016-10-06T21:05:10.387 に答える
1

ベースRでこれを試してください：

p <- c(0.5, 0.12, 0.7, 0.8, .02)
a <- c("A","B","C","D","E")
n <- length(p)
apply(expand.grid(replicate(n,list(0:1)))[n:1], 1, 
                  function(x) prod(p[which(x==1)])*prod(1-p[which(x==0)]))

#[1] 0.025872 0.000528 0.103488 0.002112 0.060368 0.001232 0.241472 0.004928 0.003528 0.000072 0.014112 0.000288 0.008232 0.000168 0.032928 0.000672 0.025872
#[18] 0.000528 0.103488 0.002112 0.060368 0.001232 0.241472 0.004928 0.003528 0.000072 0.014112 0.000288 0.008232 0.000168 0.032928 0.000672
于 2016-10-07T06:54:24.310 に答える