Here's an implementaiton of poisson distributed pseudo-random
numbers that seems to be generating the same kind of distribution
as poissonrnd from stats.
pssrnd0=:4 :0"0
m=. >.0>.x
n=. [EMAIL PROTECTED]@[EMAIL PROTECTED]&0 1.3* >:&.(p:inv) m (%:@<. * >.) y
'r t'=. ((,: {:"[EMAIL PROTECTED])~ >:&x i."1 1:) +/\"1 (y,m)$n
i=. I. x > t
t=. i { t
while. # i do.
m=. >.0>. 1.7 * (+/%#) v=. (i{x) - t
n=. [EMAIL PROTECTED]@[EMAIL PROTECTED]&0 >:&.(p:inv) m (%:@<. * >.) y
'r2 t2'=. ((,: {:"1)~ >:&v i."1 1:) +/\"1 (y,m)$n
b=. (i{x) > t2+t
r=. (r2+i{r) i} r
i=. b#i
t=. b#t+t2
end.
r
)
pssrnd=:4 :0"0
blksiz=. >. y <. 1e5 % x
r=. $0
while. y > #r do.
r=. r, x pssrnd0 blksiz <. y-#r
end.
)
ts 'poissonrand 100 1e4'
1.60453 690880
10 ts '100 pssrnd 1e4'
0.123605 2.89389e6
Since this approach can use lots of memory, I have arranged for it to
break the problem into blocks in some cases. Note, however,
that this routine still has significant room for improvement over
some parts of its domain.
Finally, note that I have sacrificed some independence between
differing parts of the result, but that is characteristic of all
psuedo-random number routines. And I believe that I have
arranged so that there are no obvious impacts from this
implementation choice.
--
Raul
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
