I have been studying the variations on normalrand
some more and have a few improvements and some puzzlement
about J's behavior, too. A little below I list the various
verb definitions, which all produce good results. They are
all variations on the Box-Muller (BM) definition that is so
widely known and used. The one defined as normalrand2 is
closest to the original BM definition, but see the details
below.
The original normalrand does NOT produce (a,b) pairs
of variates, but all the others do produce the variates in
pairs. The others produce (a,b) pairs of random varates in
different orders, they all produce valid pairs of variates.
The two possible pair patterns are as follows where in both
cases there are >.-:y a's and >.-:y b's, produced
altogether. Most produce the first pattern, btw.
abababababab....
aaaa..... bbbb......
The one defined as normalrand2 is most like the
original BM algorithm except that the output of normalrand2
follows the second pair pattern which is not exactly how BM
is most often described. The alternating ababab... pattern
could be produced with the commented out line in
normalrand2, but is not needed to get good variates. From my
comparisons using ts, normalrand2 is the overall best
performer, requiring about 2/3 of the time and 1/2 of the
space as does the current normalrand in statdist.ijs . Even
the current one, normalrand, does a quick job with modest
use of space, btw.
It is surprising to me that normalrand3 and
normalrand3a are not almost identical in performance and
that 3a requires more than twice the processing time and
almost twice the space that 3 takes, when the main
difference between the two is the fork 2&o. ,. 1&o. versus 2
1&o. .
It is also surprising to me that normalrand4 is 7
times slower than normalrand, (while it did gain on space
considerably). I was able to speed it up a little with the
added parentheses in normalrand5, but lost much more
relatively in space. I cannot be sure whether the slowness
is due more to the tacit aspect of 4 and 5 or in the
currying of the trig functions in conjunction with the
rand01 arguments. I would have been very pleased if these
two tacit definitions were better performers because they
are so concise in J.
Btw, if I leave out the rank "1 0 in either
normalrand4 or normalrand5, I crash J when I do ts in the
way I did these tests.
NB. ***************************************
rand01 =: [EMAIL PROTECTED] 0:
normalrand=: 3 : 0
(2 o. +: o. rand01 y) * %: - +: ^. rand01 y
)
normalrand2=: 3 : 0
n =. >. -: y
a=. %: - +: ^. rand01 n
b=. +: o. rand01 n
r1 =. a * 2 o. b
r2 =. a * 1 o. b
NB. y{.,r1,.r2 NB. not required
y{.r1,r2
)
normalrand3=: 3 : 0
n =. >. -: y
y{.,((2&o.,.1&o.) +: o. rand01 n) *"1 0 %: - +: ^. rand01 n
)
normalrand3a=: 3 : 0
n =. >. -: y
y{.,((2 1&o.)"1 0 +: o. rand01 n) *"1 0 %: - +: ^. rand01 n
)
normalrand4=: ({.,)(2 1&o."1 0@ +:@[EMAIL PROTECTED] * [:%:[EMAIL
PROTECTED]:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
normalrand5=: ({.,)(2 1&o."1 0@(+:@[EMAIL PROTECTED]) * [:%:[EMAIL
PROTECTED]:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
NB. ***************************************
ts 'normalrand 1000'
0.00134699 84480
ts 'normalrand2 1000'
0.00091691 35072
ts 'normalrand3 1000'
0.00100492 69184
ts 'normalrand3a 1000'
0.00270991 102080
ts 'normalrand4 1000'
0.00959293 21632
ts 'normalrand5 1000'
0.00741104 50432
NB. ***************************************
There are variations on the BM method which are
allegedly faster than the basic one characterized here, in
some programming languages. These faster ones finesse the
sine and cosine functions but rely on possible repetitions
of approximations using natural logs. The possible
repetitions required, might not be very consistent with
array processing, I fear, so I have not pursued these.
NB. ***************************************
Oops, the citation [1] I gave in my previous message
(below) is incorrect because I read the wrong second
author's name. The correct authors are as follows.
Box, GEB and ME Muller
On Fri, 8 Aug 2008, Brian Schott wrote:
[1] Box, GEP and GM Jenkins, "A Note on the Generation of
Random Normal Deviates," Annals of Mathematical Statistics,
Vol 29, 1958, pp 610-611.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
