Thank you for your kind offer. I will be as brief as I can and hopefuly
spare you some of the gory details.
My immediate objective is simply gaining a better understanding of the
potential of random numbers in making projections, and particularly in
estimating margins of error in my projections. In studying sample uses of
random numbers by a recognized expert in the field, I will probably improve
my methods.
One of the samples (SSJ examples by Pierre l'Ecuyer) titled "A
discrete-time inventory system' relies on a stream of 1 000 000 Poisson
random numbers(500 tests of 2000 days) where daily purchases are normally
simulated with a Poisson distribution. His examples are written in Java and
as far as I can see the generating of these millions of numbers results form
a hidden call to a DLL probably written in C and streamed into the java
environment. You can google his sample code.
My java test is completed in 0.45 seconds. The following rather crude
adaptation of this code is as follows. It runs in a dismal 65 seconds .
inventory=: 3 : 0
NB. A discrete-time inventory system
res=:500$0
for_t. i.500 do.
res=.( 2000 profits(200,100,2,0.1,10,1,0.95,80,200))t}res
end.
dstat res NB. the orignal also shows variability using a student.
)
profits=: 4 : 0
NB. average profit during m years
'i l c h K k p s S'=.y
NB.x: number of years
NB.i: initial
m=.x[pf=.0
ps=.poissonrand l m
for_t. i.m do.
v=. i<.t{ps NB. sales Poisson with mean l
r=. (S-i-v)*(s>i-v)*.p>?0 NB. bring inventory to level S if it fals
below s with probability p of receiving merchandise
pr=.pr+(c*v)-((i-v)*h)+(r>0)*K+k*r NB. profit for day
i=.0>.i-v+r
end.
pr%m NB. average daily profit
)
inventory=: 3 : 0
NB. A discrete-time inventory system
res=:500$0
for_t. i.500 do.
res=.( 2000 profits(200,100,2,0.1,10,1,0.95,80,200))t}res
end.
dstat res
)
This is over 100 times slower than the java version. I can ease the pain
somewhat by precomputing a Poisson cumulative distribution function, and
cump=.0,+/\poissondist 100 200 NB.cumulative Poisson
Then i remove the request for Poisson random numbers and use instead an
index search and ? :
replace the line : ps=.poissonrand l m
by: ps=.cump I. ?m#0
Thus, the timing is reduced from 65 secs to a more tolerable 16 seconds.
The quality of the result is ok as the target average produced by the java
routine is 84.983 vs 84.982, with an identical standard deviation.
This simple solution is possible because these Poisson random numbers are
integers, and this limits the required size of the above cumulative
distribution without harming the accuracy. Things would probably be
different for other distributions.
At least in this case, it would probably be a good idea to special case
poissonrand when a large number of random numbers is requested.
Cheers.
Robert Cyr
On Thu, Aug 7, 2008 at 5:42 PM, Brian Schott <[EMAIL PROTECTED]> wrote:
> Robert,
>
> If you provide more information about your
> application, you might be surprised at ways this forum could
> help. It may be that the statdist generator can be applied
> more effectively than you are thinking, or a different
> generator may be better for you, for example.
>
> On Thu, 7 Aug 2008, Robert Cyr wrote:
>
> + Performance of ? is fine, but statdist random number generators (normal,
> + poisson, gamma) are just not designed for large scale applications.
> +
> + Using a DLL or R is required.
> +
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
