> 5*y should be statistically valid > for far larger samples than our psuedo-random number > generator would be good for.
Can you please explain that? The default RNG used in ?y has a period of _1+2^19937 (4.32e6001). http://www.jsoftware.com/help/dictionary/d640.htm ----- Original Message ----- From: Raul Miller <[EMAIL PROTECTED]> Date: Friday, August 8, 2008 20:54 Subject: Re: [Jprogramming] general Gamma distribution To: Programming forum <[email protected]> > On Fri, Aug 8, 2008 at 4:53 PM, Brian Schott > <[EMAIL PROTECTED]> wrote: > > In Fishman's book > the first algorithm is like yours, > > but does not seem to guess at 'n=. 5*x' > > That does not seem to be a serious issue. > > I implemented an approach that would cache the threshold list, and > extend it if necessary: > > coclass 'poissondist' > > create=: 3 :0 > n=: 2*y > possible=: (^-y) +/\@:* 1 */\@, y&% > thresh=: possible 1+i.y > ) > > rand=: 3 :0 > r=. thresh I. t=. y [EMAIL PROTECTED] 0 > while. #i=. I. r=#thresh do. > thresh=: possible 1+i.n=: 2*n > r=. (thresh I. i{t) i} r > end. > ) > > cocurrent 'base' > And then I exercised it rather extensively: > d100=: 100 conew 'poissondist' > rand__d100 10 > 105 91 107 96 105 86 110 104 112 125 > (several repeats) > dstat rand__d100 1e4 > sample size: 10000 > minimum: 67 > maximum: 141 > median: 100 > mean: 99.7288 > std devn: 9.9441 > skewness: 0.0625085 > kurtosis: 3.01035 > (several repeats) > ts=:6!:2,7!:[EMAIL PROTECTED] > 10 ts 'rand__d100 1e5' > 0.00891996 2.23373e6 > $thresh__d100 > 201 > 10 ts 'rand__d100 1e6' > 0.065828 1.78313e7 > $thresh__d100 > 201 > 10 ts 'rand__d100 1e7' > 0.643403 2.85218e8 > $thresh__d100 > 201 > > That's after generating over a hundred million random > numbers... And to properly test that code, I had to reduce > my initial threshold list seed from 5*y to y, and I decided > that 2*y was plenty. 5*y should be statistically valid > for far larger samples than our psuedo-random number > generator would be good for. ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
