[R] Sum of Bernoullis with varying probabilities

[Deepayan Sarkar]

On 10/6/06, Ted Harding <Ted.Harding at nessie.mcc.ac.uk> wrote:
> Many thanks for your comments, Deepayan; and I liked your
> recursive solution! Fun indeed.
>
> Just a comment (below) on one of your comments (the rest
> snipped).
>
> On 06-Oct-06 Deepayan Sarkar wrote:
> > On 10/6/06, Ted Harding <Ted.Harding at nessie.mcc.ac.uk> wrote:
> >> Hi Folks,
> >>
> >> Given a series of n independent Bernoulli trials with
> >> outcomes Yi (i=1...n) and Prob[Yi = 1] = Pi, I want
> >>
> >>   P = Prob[sum(Yi) = r]  (r = 0,1,...,n)
> >>
> >> I can certainly find a way to do it:
> >>
> >> Let p be the vector c(P1,P2,...,Pn).
> >> The cases r=0 and r=n are trivial (and also are exceptions
> >> for the following routine).
> >>
> >> For a given value of r in (1:(n-1)),
> >>
> >>   library(combinat)
> >>   Set <- (1:n)
> >>   Subsets <- combn(Set,r)
> >>   P <- 0
> >>   for(i in (1:dim(Subsets)[2])){
> >>     ix <- numeric(n)
> >>     ix[Subsets[,i]] <- 1
> >>     P <- P + prod((p^ix) * ((1-p)^(1-ix)))
> >>   }
> >
> > Don't you want to multiply this with factorial(r) * factorial(n-r)? I
> > assume combn() gives all combinations and not all permutations, in
> > which case you are only counting
> >
> > prod((p^ix) * ((1-p)^(1-ix)))
> >
> > once instead of all the different ways in which ix could be permuted
> > without changing it.
>
> You had me worried for a moment -- the interplay between perms
> and combs is always a head-spinner -- but since I'd verified
> already that I get sum(P)=1 when I do this over all values of r,
> I had to conclude that your comment must be incorrect.
>
> In fact, if you consider the event "r out of the n gave Y=1",
> this can be decomposed into disjoint sub-events
>
>   subset {1,2,3,...,r} gave Y=1, the rest 0
>   subset {1,2,3,...,(r-1),(r=1)} gave Y=1, the rest 0.
>   ....
>
> and so on over all choose(n,r) subsets, and the probability of
> "r out of n" is the sum of the probabilities of the sub-events.
> For any one of these (say the first above), the probability is
>
>   P(Y1=1 & Y2=1 & ... & Yr=1 & Y[r+1]=0 & ... & Y[n]=0)
>
> which is simply the product of their probabilities. No further
> permutation is involved.

You are right. The key point I was missing is that the number of
distinct _permutations_ of r 1's and (n-r) 0's is choose(n, r) (the
places to put the 1's). Sorry for the red herring.

Deepayan