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) >> >>