Brian Schott wrote:
>
> 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
>
>
normalrand5 =:({.,)(2 1&o."1 0@(+:@[EMAIL PROTECTED]) * [:%:[EMAIL
PROTECTED]:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
One variation on normalrand5:
normalrand5a=:({.,)(2 1&o."1 0@:+:@[EMAIL PROTECTED] * [:%:[EMAIL
PROTECTED]:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
The fastest and shortest tacit code I could find is:
normalrand6 =:({.,)(%:@[EMAIL PROTECTED]:@[EMAIL PROTECTED] * 2 1 o."1
0(o.2)*rand01)@>[EMAIL PROTECTED]:
On my computer I get:
ts'normalrand3 1000000'
0.257033785 67112512
ts'normalrand4 1000000'
1.934153642 20972672
ts'normalrand5 1000000'
1.296739949 50332928
ts'normalrand5a 1000000'
1.884908836 20972672
ts'normalrand6 1000000'
1.135214671 46138624
As you already noticed, some tacit routines crash the J session.
One way to do it is to give a large argument:
ts'normalrand3 100000000' NB. non-tacit
|out of memory: normalrand3
| y{.,((2&o.,.1&o.)+:o.rand01 n)*"1 0%:-+: ^.rand01 n
ts'normalrand4 100000000' NB. crashes the J session
The most extreme case I've encountered is:
normalrandX=:({.,)(%:@[EMAIL PROTECTED]:@[EMAIL PROTECTED] * 2 1&o."0 [EMAIL
PROTECTED]:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
normalrandY=:({.,)(%:@[EMAIL PROTECTED]:@[EMAIL PROTECTED] * 2 1 o."0 1
+:@[EMAIL PROTECTED])@>[EMAIL PROTECTED]:
The former works fine, the later crashes the J session for any value of the
argument.
--
View this message in context:
http://www.nabble.com/generating-standard-normal-variates-tp18898408s24193p18989334.html
Sent from the J Programming mailing list archive at Nabble.com.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
