Raoul said: the polynomial 81 1.25&p. should be more than generous. I agree about the polynomial. My test produced a less generous 20 1.74 _0.002p. after reducing the comparison tolerance. I increased the tolerance in the hope of reducing the extreme random sample results that seem to pop up occasionally.
limit1=:3 : '(+/\(^-y)**/\1,y%}.i.500*y)i.!.(2^_35)1.0' (limit1"0 w)%.(w=.1}.i.150)^/(0 1 2) 20.1244 1.73562 _0.00204 The choice of the "best" polynomial is really not critical. I did not have much luck looking for more information concerning the behavior of I. The problem rarely occurs and apparently only in the extreme ends of the distribution. It would probably not do much harm. On Mon, Aug 11, 2008 at 6:50 PM, Raul Miller <[EMAIL PROTECTED]> wrote: > I have some more observations on this poisson distribution > threshold function: > prthre=: [:~.!.0 [EMAIL PROTECTED]@[ +/\@:* 1 */\@, %&(1+i.1000) > > First off, the plots > plot [EMAIL PROTECTED]"0 i. 713 > plot [EMAIL PROTECTED]"0 i. 800 > are interesting. The first shows a nearly linear trend and the second > shows a catastrophic > decay when we hit this limits of our numeric representation. > > Anyways, if we are concerned about efficiency, maybe we should not > be generating a full 1000 terms when they are not needed. In other words, > since the curve is nearly linear, maybe we should use a first degree > polynomial to estimate the number of terms we need. > > ([EMAIL PROTECTED]"0 i.713) %. (i.713)^/0 1 > 69.4126 1.24226 > >./ ([EMAIL PROTECTED]"0 i.713) - 69.4126 1.24226 p. i. 713 > 11.388 > 11.388 + 69.4126 > 80.8006 > > In other words, the polynomial 81 1.25&p. should be more > than generous. In fact > >./ ([EMAIL PROTECTED]"0 i.713) - 81 1.25 p. i. 713 > _2.25 > > we can use 79 1.25&p. and still have some room to spare. > However, I am going to use 80, because all I am really going > for here is a safe estimate: > > pth=: [: ~.!.0 [EMAIL PROTECTED]@[ +/\@:* 1 */\@, ]% 1 + 80 1.25 i.@>[EMAIL > PROTECTED] ] > > *./ (pth -: prthre)"0 i. 713 > 1 > > This variation produces the same results as before, for integers in > the domain of interest. > > And in some cases this speedup seems significant: > > ts=: 6!:2,7!:[EMAIL PROTECTED] > 10000 ts 'pth 2' > 2.67473e_5 4032 > 10000 ts 'prthre 2' > 0.000248684 25600 > > I could tune this further, using a slightly different polynomial, > for example, but I would expect to get into diminishing returns > rather quickly. > > However, I should get around to looking at normal distributions, > for cases where the floating point exponent of ieee 64 bit floating > point numbers will not go low enough to properly represent the > necessary terms. For example: > > pth 1000 > 0 > > FYI, > > -- > Raul > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
