I'd like to plot the probability distribution for a set of binomial trials in R, the catch is that each trial has an independent probability of success (which I have in vector form).

So, in R if I were doing this with coin tosses: plot(x,dbinom(x,10,.5)) would work, where x is 0:10. This shows the probability distribution by plotting number of successful trials on the x-axis with percent of the time that specific results is achieved as the y (4 successes is 20.5%).

That said, how do I plot the same graph for 10 discrete trials with varying probabilities. For instance if odds<-c(.1,.8,.2,.2,.3,.7,.9,.99,.05,.5), was the probability of success for each independent toss, and I wanted to see a probability distribution?


Here is an approach that uses combn to calculate the probabilities:

mybinom <- function(x, probs) {
    n <- length(probs)
    if(x == 0) return( prod(probs) )
    sum(combn( n, x, FUN = function(x) prod( probs[x] ) * prod( 1-probs[-x] )  
        ) )

# test
out1 <- sapply(0:10, mybinom, probs=rep(.5,10))
out2 <- dbinom(0:10, 10, .5)
round( out1 - out2, 10 )

# plot

plot( 0:10, sapply( 0:10, mybinom, probs=c(.1,.8,.2,.2,.3,.7,.9,.99,.05,.5) ) )

You could also use Vectorize if you want it to work more like the regular vectorized functions.


Not the answer you're looking for? Browse other questions tagged or ask your own question.