You asked a
similar question some time ago. Actually, it's basically the same thing. You had events with multiple outcomes, while now you have just two outcomes for each event. In any case, you want the sum of them.
As BruceZ showed, this is a
convolution problem. Luckily, R supports convolution and your problem can be solved defining a one-liner function:
Code:
genUnProb<-function(ph) {
setNames(Reduce(function(x,y) convolve(x,rev(y),type="open"),Map(c,1-ph,ph)),0:length(ph))
}
where the argument ph is a vector with the heads probabilities. Also a one-liner for a simulation:
Code:
genUnProbSim<-function(ph,nsim=100000) {
setNames(tabulate(colSums(matrix(rbinom(nsim*length(ph),1,ph),ncol=nsim))+1)/nsim,0:length(ph))
}
Now we get the results with your values:
Code:
ph<-c(.54,.59,.66,.50,.49,.48,.45)
convResult<-genUnProb(ph)
simResult<-genUnProbSim(ph)
cbind(convResult,simResult)
# convResult simResult
#0 0.004676563 0.00445
#1 0.038610441 0.03924
#2 0.134457040 0.13646
#3 0.256406970 0.25401
#4 0.289556183 0.29007
#5 0.193854783 0.19301
#6 0.071310213 0.07168
#7 0.011127806 0.01108
As you can see, there is a good agreement between the convolution and the simulation results (if you run the code, you'll have slightly different values for simulation result of course).
