다른 확률로 일련의 베르누이 시행의 가능한 모든 조합 확률을 열거하십시오.

Question

독립 베르누이 시행 성공을위한 일련의 n 확률, 즉 p1! = p2! = ...! = pn과 같은 p1 ~ pn이 있다고 가정합니다. 각 평가판에 고유 한 이름을 지정하십시오.

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

가 나는 포아송 이항 분포 함수를 사용하는 등 CDF, PMF를 찾을 수 스택 교환 (예를 들면, here 및 here)을 검색하지 알고있다.

내가 관심있는 부분은 성공과 실패의 모든 가능한 조합의 정확한 확률입니다. (나는 확률 트리를 그린 일예 경우에, 나는 각 지점의 끝에있는 확률을 알고 싶습니다.)

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입니다.

Answer 1

이 계산을 빠르게하는 핵심은 로그 확률 공간에서 수행하여 트리의 각 분기에 대한 곱이 행렬 곱셈의 내부 합으로 계산 될 수있는 합계가되도록합니다. 이 방식으로 모든 브랜치를 벡터화 된 방식으로 함께 계산할 수 있습니다.

먼저 모든 브랜치의 열거를 구성합니다. 이를 위해, 우리는 R.utils 패키지에서 intToBin 함수를 사용 n는 베르누이 변수의 수는

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

. 각 열은 확률 트리의 브랜치로부터의 결과를 인 매트릭스

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" 

결과하여 예를 들어 n=5. 값이 다르게 log(p) 경우 enum.branches=="1" 및 log(1-p) 여기서 지금

는 enum.branches 같은 크기의 로그 확률 행렬을 구성하는 것을 사용한다. 데이터를 들면, 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.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.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의 맥북 (2009 년경)에 있습니다.계산 자체는 매우 빠르다는 점에 유의해야합니다. 그것은 시간의 대부분을 차지하는 확률 트리의 가지들의 열거입니다 (나는 그 안에서 unlist을 추측 할 것입니다). 지역 사회가 제안한 다른 제안에 대한 제안은 인정 될 것입니다.

Answer 2

그냥 기본 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.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.0.241472 0.004928 0.003528 0.000072 0.014112 0.000288 0.008232 0.000168 0.032928 0.000672