Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Re: Factor 2 error in Inverse Laplace Transform
  • Next by Date: Re: Re: Re: normal distribution random number generation
  • Previous by thread: Re: Re: Re: Re: normal distribution random number generation
  • Next by thread: Re: Re: Re: normal distribution random number generation