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.
--
Raul
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
