Re: [R] Sum of Bernoullis with varying probabilities

From: Deepayan Sarkar <deepayan.sarkar_at_gmail.com>
Date: Fri 06 Oct 2006 - 20:46:38 GMT

On 10/6/06, Ted Harding <Ted.Harding@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.

> 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.

I don't see how. Someone has to evaluate the n-choose-r distinct terms to be added, and your code seems to be doing that as directly as possible. Just for fun, here's a cute but very inefficient implementation using recursion (even n=8 is very slow):

u <-

    function(n, r, p)
{

    if (r == 0) return(prod(1-p))
    if (r == n) return(prod(p))
    ans <- 0
    for (i in 1:n)
    {

        ans <-
            ans +
                p[i] * u(n-1, r-1, p[-i]) +
                    (1 - p[i]) * u(n-1, r, p[-i])
    }
    ans / n
}

-Deepayan

>
> And also wondering if someone has implemented it!
>
> Statutory Declaration: I have been on to R Site Search.
>
> Best wishes to all,
> Ted.



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat Oct 07 07:01:14 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 09 Oct 2006 - 15:30:10 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.