Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51332] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Thu, 14 Oct 2004 06:37:26 -0400 (EDT)
- References: <20041012184801.53136.qmail@web60509.mail.yahoo.com> <4AC453EE-1C99-11D9-BF82-000A95B4967A@mimuw.edu.pl>
- Sender: owner-wri-mathgroup at wolfram.com
Thanks to Sean Kim I now realize that creating this fix is a lot harder
that I ever imagined. The main culprit is somethig otherwise very
useful: Mathematica's packed array technology. The essential point, I
think, is this:
When Mathematica 'sees" that it is asked to generate "sufficiently
large' array, where sufficiently large means large enough for the
PackedArray technology to be used, it "ignores" any custom rules for
Random[] and reverts to using the built in Random[] generator.
Actually, I suspect this is not exactly what happens; rather than
"reverting" probably Mathematica uses special code for generating
packed arrays containing random numbers, but since this code uses the
same random number generator as random[] the effect is the same as if
the user defined rules were being ignored. What makes the problem
harder to overcome is that the same principle applies for nested lists,
but it is harder to work out precisely when the PackedArray technology
is going to be used.
What follows is a rather long explanation of what is going on including
the history of my recent attempts to solve this problem. At the end
there current best "solution'. It results in a substantial loss of
speed compared with simple usage of the unmodified Random[] function.
While it seems clear to me that no satisfactory solution can be found
until WRI provides a built in one, I hope someone will think of some
improvements to this temporary fix, as i think the matter is of some
importance.
First let us choose a particular random seed that we shall use
throughout to compare "random" outcomes with different definitions of
Random[].
SeedRandom[5]
standardValue=Random[]
0.786599
Next we shall redifine Random[], according to the original idea of
Daniel Lichblau, in the form proposed by Bobby Treat and including a
suggestion of Ray Koopman.
Unprotect[Random];
With[{m1=1/(2.^30-1.),m2=2^30-2},randomSubstitutionFunction=
Compile[{},((Random[Integer,m2]+.5)*m1+Random[Integer,m2])*m1];
Random[]:=randomSubstitutionFunction[]]
With this definition we now get:
SeedRandom[5];Random[]
0.66205
At this point almost all of us believed we had a fix, until Mark Fisher
pointed out that;
In[6]:=
SeedRandom[5];First[a=Table[Random[],{249}]]
0.66205
SeedRandom[5];First[b=Table[Random[],{250}]]
0.786599
we got the standardValue again. The explanation seems to lie here:
Developer`PackedArrayQ[a]
False
Developer`PackedArrayQ[b]
True
As I already stated above whenever Mathematica generates a packed array
it reverts to using the built in Random[] no matter what rules we
define. Now for lists the point at which the packed array technology
enters seems to be at length 250. So I thought I could solve the
problem by adding a rule for random:
Random /: Table[Random[], {n_}] /; n ³ 250 := Developer`ToPackedArray[
Flatten[{Table[Random[], {i, 1, Quotient[n, 249]}, {j, 1, 249}], \
Table[Random[], {i, 1, Mod[n, 249]}]}]]
I thought I had it solved:
SeedRandom[5];First[b=Table[Random[],{300}]]
0.66205
No problem. Then Sean Kim sent me a message pointing out that things
were still not working. At first I thought he must be using the old
package (which in any case was full of bugs) but eventually I cam to
realize that the problem was still there:
SeedRandom[5];First[b=Table[Random[],{1000}]]
0.786599
Developer`PackedArrayQ[b]
True
The point is that I had thought that I could avoid the problem by
constructing long lists as lists of lists of 249 elements (plus a
shorter list), but Mathematica anticipates this and at some point the
PackedArray technology again kicks in and the problem returns. To see
it more clearly let's look at the following way of generating a table:
We shall first clear Random and redefine it form the beginning;
Clear[Random]
With[{m1=1/(2.^30-1.),m2=2^30-2},randomSubstitutionFunction=
Compile[{},((Random[Integer,m2]+.5)*m1+Random[Integer,m2])*m1];
Random[]:=randomSubstitutionFunction[]]
Now compare:
n=10^4;SeedRandom[5];
First[a=Flatten[{Table[Table[Random[],{j,1,249}],{i,
1,Quotient[n,249]}],Table[Random[],{i,1,Mod[n,249]}]}]]
0.66205
n=10^5;SeedRandom[5];
First[b=Flatten[{Table[Table[Random[],{j,1,249}],{i,
1,Quotient[n,249]}],Table[Random[],{i,1,Mod[n,249]}]}]]
0.786599
So we clearly see that for lists of lists the problem again occurs,
although for much larger lists: of size around 10000.
What about the solution? Well, at the moment the only thing that comes
to my mind is to generate tables using something like this:
SeedRandom[5];Timing[First[Flatten[ReleaseHold[{Table[
Hold[Table[Random[],{j,1,249}]],{i,1,Quotient[n,249]}],Table[Random[],{i
,\
1,Mod[n,249]}]}]]]]
=
{1.57 Second,0.147012}
if we clear Random we get a different answer, a lot faster:
Clear[Random]
SeedRandom[5];
Timing[First[a=Flatten[ReleaseHold[{Table[Hold[Table[Random[],{
j,1,249}]],{i,1,Quotient[n,249]}],Table[Random[],{i,1,Mod[n,249]}]}]]]]
{0.28 Second,0.160058}
Note that this is no longer what we called standardvalue, since this
time the first element of the list is not the first random number that
was generated. But at least for now it seems to me that I have got a
rather slow fix. I am no longer convinced that this is necessarily the
best way to approach the problem: there should be a faster way to
"randomize' the uniform random number generator than this. But we
obviously need WRI to do something about this matter soon.
Andrzej Kozlowski
Andrzej Kozlowski
Chiba, Japan
http://www.akikoz.net/~andrzej/
http://www.mimuw.edu.pl/~akoz/
On 13 Oct 2004, at 06:54, Andrzej Kozlowski wrote:
> This is exactly what it is mean to do. The problem was due to the
> usage of the SWB algortithm by the built in Random[] for generating
> random reals. RandomReplace makes Mathematica use the Wolfran CA
> random number generator, which Mathematica uses for generating random
> integers instead of the SWB. This "cures" the problem reported int
> hat message. Have you got any reason to think it does not?
>
> Andrzej
>
>
> On 13 Oct 2004, at 03:48, sean kim wrote:
>
>> *This message was transferred with a trial version of CommuniGate(tm)
>> Pro*
>> Hi andrej
>>
>> thank you so much for making your package available.
>>
>> I think I may have found another problem... ("may have
>> found" being operative phrase, or non-operative for
>> that matter)
>>
>>