normal distribution random number generation

[Chris Kunkel]

Hello,

Recently I have encountered a problem in Mathematica's normal
distribution random number generator. The problem arises when I look
at the distribution of averages of a list of theses numbers. That is,
I generate 1000 numbers and take their average. I do this a number of
times and the plot a frequency distribution. Consistently it seems to
be skewed positive. Specifically, the number of occurrences less than
3 standard deviations is consistent with the expected number, but the
number greater than 3 is always larger (in some cases about twice the
expected number).

I was wondering if anyone else has noticed this and knows of a fix.
Also, does anyone have code for a good quality normal distribution
random number generator?

Chris

[Andrzej Kozlowski]

The problem is not actually with the way Mathematica's
NormalDistribution but with the uniformly distributed Random[]
function.
NormalDistribution itself is generated (in the package
Statistics`NormalDistribution`) by means of the very classical so
called Boox-Muller method (although actually the Marsaglia variant
below works better). You can try downolading my little
RandomReplacement package form one of my web sites:

http://www.akikoz.net/~andrzej//Mathematica/ in Japan

or

http://www.mimuw.edu.pl/~akoz/Mathematica/ in Poland

The package is based on the post of daniel Lichtblau, which explains
the nature of the problem:

<http://library.wolfram.com/mathgroup/archive/2000/May/msg00088.html>

The package is named init.m and can be loaded automatically at the
start of each Mathematica session.
However, in fact if you are only concerned with normal distributions it
may be enough to the following. First load the (unmodified)
NormalDistribution package using

<<Statistics`NormalDistribution`

and then evaluate the code below. It will replace the Box-Muller
generator by the Marsaglia one. I have found that this is usually
enough. But if your problem still remains try the RandomReplacement
package. (Or you can use both of course). Here is the Marsaglia code
for normal distribution:

Statistics`NormalDistribution`Private`normal=
Compile[{{mu,_Real},{sigma,_Real},{q1,_Real},{q2,_Real}},
Module[{va=1.,vb,rad=2.,den=1.},
While[rad>=1.,va=2.*Random[]-1.;vb=2.*Random[]-1.;
rad=va*va+vb*vb];
den=Sqrt[-2.*(Log[rad]/rad)];
mu+sigma*va*den]];

On 5 Oct 2004, at 17:36, Chris Kunkel wrote:

> *This message was transferred with a trial version of CommuniGate(tm)
> Pro*

[Bill Rowe]

On 10/5/04 at 4:36 AM, cjku...@marauder.millersville.edu (Chris
Kunkel) wrote:

>Recently I have encountered a problem in Mathematica's normal
>distribution random number generator. The problem arises when I
>look at the distribution of averages of a list of theses numbers.
>That is, I generate 1000 numbers and take their average. I do this
>a number of times and the plot a frequency distribution.
>Consistently it seems to be skewed positive. Specifically, the
>number of occurrences less than 3 standard deviations is consistent
>with the expected number, but the number greater than 3 is always
>larger (in some cases about twice the expected number).

>I was wondering if anyone else has noticed this and knows of a fix.
>Also, does anyone have code for a good quality normal distribution
>random number generator?

I don't think this is a problem with the code for generating normal deviates. Instead, I think it is a problem with the algorithm to generate uniform deviates used by the code to generate normal deviates.

For reals, Mathematica uses a subtract with borrow algorithm. Unforunately, this algorithm is flawed. This has been discussed previously in this forum. The work around is to generate random integers between 0 and 2^30-1. For these integers, Mathematica uses the Wolfram rule 30 cellular automaton which does not seem to have any problems.

If you do the following:

Unprotect[Random];

Random[]:=(Random[Integer,2^30-1]/2^30.+Random[Integer,2^30-1])/2^30.

Random[Real,{a_Real,b_Real}]:=a+Random[]*(b-a)

Random[Real,b_Real]:=Random[Real,{0,b}]

Random[Real]:=Random[Real,{0,1}]

Protect[Random];

you will have modified Random to use the Wolfram rule 30 cellular automaton and avoid the subtract with borrow algorithm. The main consequence of this is Random will now be considerably slower.

The other alternative would be to create you own pseudo-random number generator and patch Random as above with it. But be aware, creating pseudo-random number generators is not an easy task. Unless you know what you are doing, you are quite likely to wind up with something even worse than the subtract with borrow algorithm Mathematica uses.
--
To reply via email subtract one hundred and four

[Ray Koopman]

Bill Rowe <readn...@earthlink.net> wrote in message
news:<ck0ccp$o1u$1...@smc.vnet.net>...
> [...]

> you will have modified Random to use the Wolfram rule 30 cellular
> automaton and avoid the subtract with borrow algorithm. The main
> consequence of this is Random will now be considerably slower.

> [...]

If time is an issue, you might want to consider generating integers
on 0...2^n-2 instead of 0...2^n-1. It's always much faster. And if
you're willing to spend a little of the time you've saved, you can
add a half and avoid ever having to worry about getting a zero.

In[1]:= ToString[TableForm[Table[With[{m1 = 2^n - 1, m2 = 2^n - 2},
{n, First[Timing[Do[Random[Integer,m1],{1*^6}]]]/.Second->1.,
First[Timing[Do[Random[Integer,m2],{1*^6}]]]/.Second->1.}],
{n,2,30}],TableSpacing->{0,2}]]

Out[1]= 2 1.96 1.42
3 2.12 1.5
4 2.38 1.61
5 2.66 1.73
6 2.91 1.86
7 3.16 2.
8 3.41 2.1
9 3.68 2.19
10 3.92 2.35
11 4.21 2.56
12 4.5 2.68
13 4.79 2.82
14 5.07 3.02
15 5.34 3.08
16 5.56 3.26
17 5.84 3.38
18 6.09 3.53
19 6.33 3.64
20 6.57 3.77
21 6.84 3.87
22 7.1 4.03
23 7.33 4.2
24 7.63 4.25
25 7.89 4.37
26 8.15 4.56
27 8.4 4.61
28 8.56 4.79
29 8.95 4.95
30 9.16 5.07

In[2]:= ran1 = With[{m = 2.^-30, m1 = 2^30 - 1},
Compile[{},(Random[Integer,m1]*m + Random[Integer,m1])*m]];

In[3]:= ran2 = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2},
Compile[{},(Random[Integer,m2]*m1 + Random[Integer,m2])*m1]];

In[4]:= ran2h = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2},
Compile[{},((Random[Integer,m2]+.5)*m1+Random[Integer,m2])*m1]];

In[5]:= First/@{Timing@Do[ran1[],{1*^5}],Timing@Do[ran2[],{1*^5}],
Timing@Do[ran2h[],{1*^5}]}
Out[5]= {2.03 Second, 1.05 Second, 1.08 Second}

[DrBob]

I'm trying to combine that idea with Andrzej Kozlowski's recent fix for Random, and here's what I came up with:

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[]
]
Random[Real, {a_Real, b_Real}] := a + Random[]*(b - a)
Random[Real, b_Real] := Random[Real, {0, b}]
Random[Real] := Random[Real, {0, 1}]
Random[Complex, {a_Complex | a_Real | a_Integer, b_Complex | b_Real | \
b_Integer}] := a + Random[]*Re[(b - a)] + Random[]*Im[(b - a)]*I
Random[Complex] := Random[Complex, {0, 1 + I}]
Protect[Random];

I wanted NOT to use a Global (randomSubstitutionFunction) for the Compiled function, but I haven't stumbled on a way to manage it.

Bobby

On Sat, 9 Oct 2004 04:18:30 -0400 (EDT), Ray Koopman <koo...@sfu.ca> wrote:

> Bill Rowe <readn...@earthlink.net> wrote in message
> news:<ck0ccp$o1u$1...@smc.vnet.net>...
>> [...]

>> you will have modified Random to use the Wolfram rule 30 cellular
>> automaton and avoid the subtract with borrow algorithm. The main
>> consequence of this is Random will now be considerably slower.

--
Dr...@bigfoot.com
www.eclecticdreams.net

[Andrzej Kozlowski]

It's probably best to make a proper package out of this. The function
randomSubstitutionFunction should be in the Private` context, invisible
to the user while Random[] will be exported, thus replacing the built
in Random[]. One can still have the package load at startup, if one
choses to. I think I will put this version on my site instead of the
earlier fix (with an acknowledgement of all the contributors of course
;-) )

Also, I am going to include a similar short package that will replace
the Box-Muller normal random generator in the
Statistics`NormalDistribution` package by the Marsaglia version. As I
have mentioned before, I have found that when generating normally
distributed data this replacement seems to make up for the weakness of
the Mathematica Random[] function and is actually a little faster than
the Box-Muller normal distribution generator.

Andrzej Kozlowski

On 10 Oct 2004, at 14:21, DrBob wrote:

>
> I'm trying to combine that idea with Andrzej Kozlowski's recent fix
> for Random, and here's what I came up with:
>
> 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[]
> ]
> Random[Real, {a_Real, b_Real}] := a + Random[]*(b - a)
> Random[Real, b_Real] := Random[Real, {0, b}]
> Random[Real] := Random[Real, {0, 1}]
> Random[Complex, {a_Complex | a_Real | a_Integer, b_Complex | b_Real | \
> b_Integer}] := a + Random[]*Re[(b - a)] + Random[]*Im[(b - a)]*I
> Random[Complex] := Random[Complex, {0, 1 + I}]
> Protect[Random];
>
> I wanted NOT to use a Global (randomSubstitutionFunction) for the
> Compiled function, but I haven't stumbled on a way to manage it.
>
> Bobby
>
> On Sat, 9 Oct 2004 04:18:30 -0400 (EDT), Ray Koopman <koo...@sfu.ca>
> wrote:
>
>> Bill Rowe <readn...@earthlink.net> wrote in message
>> news:<ck0ccp$o1u$1...@smc.vnet.net>...
>>> [...]

>>> you will have modified Random to use the Wolfram rule 30 cellular
>>> automaton and avoid the subtract with borrow algorithm. The main
>>> consequence of this is Random will now be considerably slower.

[Mark Fisher]

FYI: I've just a little testing and I find that Mathematica ignors the user
defined rules for Random in Table[Random[],{n}] when n >= 250.

[Andrzej Kozlowski]

I have put on both my web sites a new version of the Random[] fix
incorporating the ideas of Ray and Bob. The package is now called
RadnomReplacement.m and is a proper Mathematica package, which defines
its own context: RandomReplacement`. When loaded, it will replace the
built in Random[] function with one based on the Wolfram CA algorithm.

I have also put on the same page another "package' called
MarsagliaNormal.m which changes the definition of the normal
distribution generator in the Statistics`NormalDistribution` package to
that of Marsaglia's. This is not a "true package", in that it does not
define any new contexts, and has to be read in using Get and its full
name including the extension *.m. It can be read in after or before the
Statistics`NormalDistribution` package. In the latter case it will
automatically load the package and seamlessly redifine the normal
generator. Doing this seems to have two advantages. It seems to be
marginally faster than the Box-Muller generator, though this may be
hardware dependent. It also seems to "eliminate' the problem caused by
"insufficient randomness' of Random[] so if one is dealing only with
normally distributed quantities it may not be necessary to read in the
RandomReplacement package.

Andrzej

On 10 Oct 2004, at 14:21, DrBob wrote:

>
> I'm trying to combine that idea with Andrzej Kozlowski's recent fix
> for Random, and here's what I came up with:
>
> 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[]
> ]
> Random[Real, {a_Real, b_Real}] := a + Random[]*(b - a)
> Random[Real, b_Real] := Random[Real, {0, b}]
> Random[Real] := Random[Real, {0, 1}]
> Random[Complex, {a_Complex | a_Real | a_Integer, b_Complex | b_Real | \
> b_Integer}] := a + Random[]*Re[(b - a)] + Random[]*Im[(b - a)]*I
> Random[Complex] := Random[Complex, {0, 1 + I}]
> Protect[Random];
>
> I wanted NOT to use a Global (randomSubstitutionFunction) for the
> Compiled function, but I haven't stumbled on a way to manage it.
>
> Bobby
>
> On Sat, 9 Oct 2004 04:18:30 -0400 (EDT), Ray Koopman <koo...@sfu.ca>
> wrote:
>
>> Bill Rowe <readn...@earthlink.net> wrote in message
>> news:<ck0ccp$o1u$1...@smc.vnet.net>...
>>> [...]

>>> you will have modified Random to use the Wolfram rule 30 cellular
>>> automaton and avoid the subtract with borrow algorithm. The main
>>> consequence of this is Random will now be considerably slower.

[DrBob]

I really didn't contribute enough to earn a mention, Andrzej.

I'm still scratching my head trying to eliminate the named helper function.

Bobby

>>>> you will have modified Random to use the Wolfram rule 30 cellular
>>>> automaton and avoid the subtract with borrow algorithm. The main
>>>> consequence of this is Random will now be considerably slower.

--
Dr...@bigfoot.com
www.eclecticdreams.net

[Andrzej Kozlowski]

It's clearly because of packed arrays:

<<Developer`

PackedArrayQ[Table[Random[],{249}]]

False

PackedArrayQ[Table[Random[],{250}]]

True

Obviously Mathematica uses a different method of generating random
packed arrays. There are various ways of getting round this (e.g.
constructing lists as joins of lists of length less than 250 etc) but
one will loose the benefits of packed arrays and with that, presumably,
a lot of speed. This seems to be something that only WRI can change.
Note however that not all situations in which one uses Random[] require
generating this type of lists.
Moreover, for problems involving the normal distribution it seems to me
that using the Marsaglia generator gives good results even with the
built in Random[]. But certainly we should hope that something will
finally be done do deal with this issue in version 6.

Andrzej

On 11 Oct 2004, at 17:24, DrBob wrote:

> He's right! How the devil does that happen?
>
> Here's a test, with Andrzej's package loaded in my Init file.
>
> First, with n=250:
>
> Quit
> SeedRandom[5]
> Table[Random[],{250}];
> Last@%
>
> 0.107874
>
> Unprotect[Random];
> Clear@Random
> SeedRandom[5]
> Table[Random[],{250}];
> Last@%
>
> 0.107874
>
> Now, with n=249:
>
> Quit
> SeedRandom[5]
> Table[Random[],{249}];
> Last@%
>
> 0.656099
>
> Unprotect[Random];
> Clear@Random
> SeedRandom[5]
> Table[Random[],{249}];
> Last@%
>
> 0.0708373
>
> Bobby

>
> On Mon, 11 Oct 2004 01:25:24 -0400 (EDT), Mark Fisher
> <ma...@markfisher.net> wrote:
>
>> FYI: I've just a little testing and I find that Mathematica ignors
>> the user
>> defined rules for Random in Table[Random[],{n}] when n >= 250.
>>
>> Andrzej Kozlowski wrote:

> --
> Dr...@bigfoot.com
> www.eclecticdreams.net
>

[Maxim]

Mark Fisher <ma...@markfisher.net> wrote in message news:<ckd5pr$50l$1...@smc.vnet.net>...

> FYI: I've just a little testing and I find that Mathematica ignors the user
> defined rules for Random in Table[Random[],{n}] when n >= 250.

This is indeed weird. It is easy to reproduce this behaviour:

In[1]:=
Unprotect[Random]
Random[]:=0
Table[Random[], {249}] // Short[#, 1]&
Table[Random[], {250}] // Short[#, 1]&

Out[3]=
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, <<227>>, 0, 0, 0, 0, 0,
0, 0}

Out[4]=
{0.6238474115604462, 0.9251918536300756, <<247>>, 0.7177852283781632}

In the last case the user definitions for Random are simply ignored.
This seems to affect only Random, not other internal functions.
However, I'd say that this example is not unique. Consider:

In[1]:=
(*restart the kernel*)
Unprotect[Power];
ClearAttributes[Power, Listable];
(A_?MatrixQ) ^ p_ := MatrixPower[A, p]
E ^ (A_?MatrixQ) := MatrixExp[A];
SetAttributes[Power, Locked];
LaplaceTransform[1, t, p];
Attributes[Power]

Out[7]=
{Listable, Locked, NumericFunction, OneIdentity, Protected}

The definitions for A^p and E^A are given only to show why removing
Listable may be useful. What this example demonstrates is that calling
LaplaceTransform (only the first time, not the subsequent calls)
restores the Listable attribute, and what's more, it ignores the
attribute Locked! Since the evaluator checks the attributes before
searching for user-defined rules, it means that the definitions for
A^p and E^A break down. Overall, redefining the built-in functions
seems to be extremely unreliable.

Maxim Rytin
m...@inbox.ru

[DrBob]

First, with n=250:

0.107874

0.107874

Now, with n=249:

0.656099

0.0708373

Bobby

> FYI: I've just a little testing and I find that Mathematica ignors the user
> defined rules for Random in Table[Random[],{n}] when n >= 250.
>

--
Dr...@bigfoot.com
www.eclecticdreams.net

[Andrzej Kozlowski]

Are you using the latest version of the RandomReplacement package? (The
one I posted last night).
The earlier one had a bug and ddi not work with calls to
Random[Real,{a,b}] where a and b are integers. You had to use
Random[Real,N[{a,b}]].

I think (I hope) I fixed that in the latest version.

Andrzej

On 13 Oct 2004, at 10:06, sean kim wrote:

> *This message was transferred with a trial version of CommuniGate(tm)
> Pro*

> I know... It's werid for me... which i wanted to ask
> you all.
>
> when i used the code below I still got the low 24 bin
> value.
>
> In[1]:=<<RandomReplacement`
>
> In[2]:=
> SeedRandom[1111111 ];
> vec1=Table[Random[Real,{0,1000}],{10^6}];
> vec2=Map[If[#<44,1,0]&,vec1];
> ones=Flatten[Position[vec2,1]];
> plotvec=Map[Length,Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];
> ListPlot[Log[plotvec],PlotRange->All]
> ListPlot[plotvec,PlotRange->All]
>
> used it exactly as it is shown.
>
>
> but the below works, it seems...
>
> SeedRandom[1111111];
> vec1 = Table[Random[Integer,1000],{10^6}];
> vec2 = Map[If[#<=44,1,0]&, vec1];
> ones = Flatten[Position[vec2,1]];
> Length[ones]/10.^6
> plotvec =
> Map[Length,
> Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];
> ListPlot[
> Log[plotvec], PlotRange->All]
> ListPlot[plotvec, PlotRange->All]
>
> maybe other people can try it and see if it's just me?
> (I took out the init.m btw)
>
> sean

>
>
>
> --- Andrzej Kozlowski <ak...@mimuw.edu.pl> wrote:
>
>> *This message was transferred with a trial version
>> of CommuniGate(tm) Pro*

>> 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)
>>>
>>> http://groups.google.com/groups?
>>>
>>
> hl=en&lr=&newwindow=1&safe=off&threadm=a5t1fq%245vt%241%40smc.vnet.net&
>>
>>> rnum=1&prev=/
>>>
>>
> groups%3Fhl%3Den%26lr%3D%26newwindow%3D1%26safe%3Doff%26q%3Dchange%2Bth
>>
>>>
>>
> is%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253
>>
>>> Dcomp.soft-sys.math.mathematica
>>>
>>> shows a thread dealing with random reals and
>> getting a
>>> much smaller frequency at position 24 than others
>>> positions.
>>>
>>> Daniel posted a work around showing the use of
>> Random
>>> Integers to fix the problem, which as far I
>>> understand, the Random replacement also
>> implements.
>>>
>>> It seems the RandomReplacement doesn't take care
>> of
>>> the appearance of the lower than expected
>> frequency at
>>> position 24(or interval 24 as the OP had noted).
>>>
>>> I must be missing something.
>>>
>>> thanks all for any insights.
>>>
>>> sean

>>>
>>>
>>> --- Andrzej Kozlowski <ak...@mimuw.edu.pl> wrote:
>>>
>>>> *This message was transferred with a trial
>> version
>>>> of CommuniGate(tm) Pro*

>>>> I think I have now fixed the problem reported by
>>>> Mark Fisher. The new
>>>> version of the RandomReplacement now works like
>>>> this:
>>>>
>>>>
>>>> SeedRandom[5]
>>>>
>>>>
>>>> Timing[First[Table[Random[],{10000}]]]
>>>>
>>>>
>>>> {0.01 Second,0.786599}
>>>>
>>>>
>>>> <<RandomReplacement`
>>>>
>>>> Timing[First[Table[Random[],{10000}]]]
>>>>
>>>>
>>>> {0.02 Second,0.409222}
>>>>
>>>>
>>>> As you can see there is a loss of speed involved.
>>>> There may also be
>>>> still some as yet undiscovered problems. However,
>>>> the new package is
>>>> available form my sites, though I still should
>> write
>>>> some documentation
>>>> of the latest changes (a few bugs in the previous
>>>> version have also
>>>> been fixed).
>>>>
>>>> Andrzej Kozlowski
>>>>
>>>> Andrzej Kozlowski
>>>> Chiba, Japan
>>>> http://www.akikoz.net/~andrzej/
>>>> http://www.mimuw.edu.pl/~akoz/
>>>>
>>>>
>>>>
>>>> On 12 Oct 2004, at 14:57, Andrzej Kozlowski

>> wrote:
>>>>
>>>>> *This message was transferred with a trial
>> version
>>>> of CommuniGate(tm)
>>>>> Pro*

>>> __________________________________
>>> Do you Yahoo!?
>>> Yahoo! Mail Address AutoComplete - You start. We
>> finish.
>>> http://promotions.yahoo.com/new_mail
>>>
>>
>>
>
>
> __________________________________________________
> Do You Yahoo!?
> Tired of spam? Yahoo! Mail has the best spam protection around
> http://mail.yahoo.com
>

[sean kim]

i just took a look at Part[plotvec, 24]

the respective values, for me, are,

375 and 686.

so i did a test with the following.

in the fresh kernel, I get exactly the same thing as
you have.

but I decided to evaluate the code

In[5]:=

SeedRandom[1111111 ];
vec1=Table[Random[Real,{0,1000}],{10^6}];
vec2=Map[If[#<44,1,0]&,vec1];
ones=Flatten[Position[vec2,1]];

plotvec=Map[Length,

Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];

Part[plotvec, 24]

ListPlot[Log[plotvec],PlotRange->All]
ListPlot[plotvec,PlotRange->All]

Out[10]=
375

and then

reevaluate the codes

In[13]:= SeedRandom[5]
In[14]:= Timing[First[Table[Random[],{10000}]]]
Out[14]= {0.01 Second,0.786599}
In[15]:= <<RandomReplacement`

Timing[First[Table[Random[],{10000}]]]

Out[16]= {0.01 Second,0.611425}

and I get the different value than before(one with the
fresh kernel)

i wonder what's going on here?

sean

__________________________________
Do you Yahoo!?
Read only the mail you want - Yahoo! Mail SpamGuard.
http://promotions.yahoo.com/new_mail

[Andrzej Kozlowski]

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*

[sean kim]

I know... It's werid for me... which i wanted to ask
you all.

when i used the code below I still got the low 24 bin
value.

In[1]:=<<RandomReplacement`

In[2]:=

SeedRandom[1111111 ];
vec1=Table[Random[Real,{0,1000}],{10^6}];
vec2=Map[If[#<44,1,0]&,vec1];
ones=Flatten[Position[vec2,1]];

plotvec=Map[Length,Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];

ListPlot[Log[plotvec],PlotRange->All]
ListPlot[plotvec,PlotRange->All]

used it exactly as it is shown.

but the below works, it seems...

SeedRandom[1111111];
vec1 = Table[Random[Integer,1000],{10^6}];
vec2 = Map[If[#<=44,1,0]&, vec1];
ones = Flatten[Position[vec2,1]];
Length[ones]/10.^6
plotvec =

Map[Length,
Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];

ListPlot[
Log[plotvec], PlotRange->All]
ListPlot[plotvec, PlotRange->All]

maybe other people can try it and see if it's just me?
(I took out the init.m btw)

sean

> > http://groups.google.com/groups?
> >
>
hl=en&lr=&newwindow=1&safe=off&threadm=a5t1fq%245vt%241%40smc.vnet.net&
>
> > rnum=1&prev=/
> >
>
groups%3Fhl%3Den%26lr%3D%26newwindow%3D1%26safe%3Doff%26q%3Dchange%2Bth
>
> >
>
is%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253
>
> > Dcomp.soft-sys.math.mathematica
> >
> > shows a thread dealing with random reals and
> getting a
> > much smaller frequency at position 24 than others
> > positions.
> >
> > Daniel posted a work around showing the use of
> Random
> > Integers to fix the problem, which as far I
> > understand, the Random replacement also
> implements.
> >
> > It seems the RandomReplacement doesn't take care
> of
> > the appearance of the lower than expected
> frequency at
> > position 24(or interval 24 as the OP had noted).
> >
> > I must be missing something.
> >
> > thanks all for any insights.
> >

> > sean

> >
> >
> > --- Andrzej Kozlowski <ak...@mimuw.edu.pl> wrote:
> >
> >> *This message was transferred with a trial
> version
> >> of CommuniGate(tm) Pro*

> >> I think I have now fixed the problem reported by
> >> Mark Fisher. The new

> >> version of the RandomReplacement now works like
> >> this:
> >>
> >>
> >> SeedRandom[5]
> >>
> >>

> >> Timing[First[Table[Random[],{10000}]]]
> >>
> >>
> >> {0.01 Second,0.786599}
> >>
> >>
> >> <<RandomReplacement`
> >>
> >> Timing[First[Table[Random[],{10000}]]]
> >>
> >>
> >> {0.02 Second,0.409222}
> >>
> >>
> >> As you can see there is a loss of speed involved.
> >> There may also be
> >> still some as yet undiscovered problems. However,
> >> the new package is
> >> available form my sites, though I still should
> write
> >> some documentation
> >> of the latest changes (a few bugs in the previous
> >> version have also
> >> been fixed).
> >>

> >> Andrzej Kozlowski
> >>
> >> Andrzej Kozlowski
> >> Chiba, Japan
> >> http://www.akikoz.net/~andrzej/
> >> http://www.mimuw.edu.pl/~akoz/
> >>
> >>
> >>

> >> On 12 Oct 2004, at 14:57, Andrzej Kozlowski

> wrote:
> >>
> >>> *This message was transferred with a trial
> version
> >> of CommuniGate(tm)
> >>> Pro*

> >>>>>> On 5 Oct 2004, at 17:36, Chris Kunkel wrote:
> >>>>>>
> >>>>>>
> >>>>>>> *This message was transferred with a trial
> >> version of
> >>>>>>> CommuniGate(tm)
> >>>>>>> Pro*

[Andrzej Kozlowski]

[Andrzej Kozlowski]

[sean kim]

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)

http://groups.google.com/groups?hl=en&lr=&newwindow=1&safe=off&threadm=a5t1fq%245vt%241%40smc.vnet.net&rnum=1&prev=/groups%3Fhl%3Den%26lr%3D%26newwindow%3D1%26safe%3Doff%26q%3Dchange%2Bthis%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253Dcomp.soft-sys.math.mathematica

shows a thread dealing with random reals and getting a
much smaller frequency at position 24 than others
positions.

Daniel posted a work around showing the use of Random
Integers to fix the problem, which as far I
understand, the Random replacement also implements.

It seems the RandomReplacement doesn't take care of
the appearance of the lower than expected frequency at
position 24(or interval 24 as the OP had noted).

I must be missing something.

thanks all for any insights.

sean

--- Andrzej Kozlowski <ak...@mimuw.edu.pl> wrote:

[Andrzej Kozlowski]

I know. I now have some ideas about how to make it a bit faster, but
need to test them first. But there is no way it can ever be as fast as
the built-in Random. I do hope the decision makers at WRI take this
issue more seriously and do something about it in Mathematica 6. As
computers grow faster the problem is becoming more visible and the
number of users affected by it seems to be growing.
Andrzej

On 15 Oct 2004, at 14:57, DrBob wrote:

> The speed difference can be pretty dramatic:
>
> Quit
>
> a=First@Timing[Table[Random[Real,{0,1000}],{10^6}];]
> <<RandomReplacement`
> b=First@Timing[Table[Random[Real,{0,1000}],{10^6}];]
> b/a
>
> 0.375 Second
>
> 10.078 Second
>
> 26.8747
>
> I can probably live with it, however, in cases where I'm concerned
> about built-in Random's non-random behavior.
>
> Bobby
>
> On Thu, 14 Oct 2004 06:35:37 -0400 (EDT), Andrzej Kozlowski

> --
> Dr...@bigfoot.com
> www.eclecticdreams.net
>

[DrBob]