[R] Sum of Bernoullis with varying probabilities

[(Ted Harding)]

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)))
  }

OK, some abbreviation of the above can be achieved with
the 'apply' function (and I've spelt out the details for
clarity). But I feel the basis of the approach (i.e. 'combn')
is crude (also tending to cause problems if n is large);
and I'm wondering if there is a nicer way.

And also wondering if someone has implemented it!

Statutory Declaration: I have been on to R Site Search.

Best wishes to all,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 06-Oct-06                                       Time: 20:29:02
------------------------------ XFMail ------------------------------