2

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?

whuber
  • 281,159
  • 54
  • 637
  • 1,101
ryan
  • 21
  • 1

1 Answers1

1

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.

Greg Snow
  • 46,563
  • 2
  • 90
  • 159