[R] Sum of Bernoullis with varying probabilities

[Gabor Grothendieck]

On 10/8/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> On 10/8/06, Ted Harding <Ted.Harding at nessie.mcc.ac.uk> wrote:
> > On 08-Oct-06 Gabor Grothendieck wrote:
> > > Or perhaps its clearer (and saves a bit of space) to use apply...prod
> > > here instead of exp...log:
> > >
> > > fft(apply(mvfft(t(cbind(1-p,p,0,0,0))), 1, prod), inverse = TRUE)/5
> >
> > Yes, that's neat (in either version). With the example I have,
> > where length(p)=16, I applied your suggestion above in the form
> >
> >  v<-fft(apply(mvfft(t(cbind(1-p,p,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))),
> >               1,prod), inverse = TRUE)/17
> >
> > (making the number of columns for cbind equal to 17 = 16+1). Now,
> > comparing this result with the P I got by brute force (subsets
> > selection method), and removing the imaginary parts from v:
> >
> > cbind(r=(0:16),v=Re(v),P)
> >   r            v            P
> >   0 1.353353e-01 1.353353e-01
> >   1 3.007903e-01 3.007903e-01
> >   2 2.976007e-01 2.976007e-01
> >   3 1.747074e-01 1.747074e-01
> >   4 6.826971e-02 6.826971e-02
> >   5 1.884969e-02 1.884969e-02
> >   6 3.804371e-03 3.804371e-03
> >   7 5.720398e-04 5.720398e-04
> >   8 6.463945e-05 6.463945e-05
> >   9 5.489454e-06 5.489454e-06
> >  10 3.473997e-07 3.473997e-07
> >  11 1.607822e-08 1.607822e-08
> >  12 5.262532e-10 5.262532e-10
> >  13 1.148620e-11 1.148618e-11
> >  14 1.494031e-13 1.493761e-13
> >  15 8.887779e-16 8.764298e-16
> >  16 1.434973e-17 1.313261e-19
> >
> > so this calculation gives a better result than convolve().
> > The only "fly in the ointment" (which comes down to how one
> > sets up the arguments to cbind(), so is quite easily handled
> > in general application) is the variable number of columns
> > required according to the length of p.
>
> Try this:
>
> p <- seq(4)/8 # test data
> Re(fft(apply(mvfft(t(cbind(1-p, p, 0*p%o%p[-1]))), 1, prod), inv =
> TRUE)/(length(p)+1))
>

And here is one minor improvement (rbind() rather than t(cbind()):

p <- 1:4/8
Re(fft(apply(mvfft(rbind(1-p, p, 0 * p[-1] %o% p)), 1, prod), inverse
= TRUE) / (length(p) + 1))