On 10/8/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote: > On 10/8/06, Ted Harding <[EMAIL PROTECTED]> 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)) ______________________________________________ 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.