24 views

Skip to first unread message

Oct 5, 2004, 4:55:01â€¯AM10/5/04

to

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

Oct 6, 2004, 4:47:00â€¯AM10/6/04

to

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:

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*

Oct 6, 2004, 5:03:21â€¯AM10/6/04

to

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

Kunkel) wrote:

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

Oct 9, 2004, 4:20:43â€¯AM10/9/04

to

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.

> [...]

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}

Oct 10, 2004, 2:12:25â€¯AM10/10/04

to

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.

Oct 10, 2004, 8:55:40â€¯AM10/10/04

to

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

;-) )

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.

Oct 11, 2004, 1:30:35â€¯AM10/11/04

to

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.

defined rules for Random in Table[Random[],{n}] when n >= 250.

Oct 11, 2004, 1:31:36â€¯AM10/11/04

to

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.

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.

Oct 11, 2004, 1:32:36â€¯AM10/11/04

to

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.

Oct 12, 2004, 2:07:19â€¯AM10/12/04

to

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

>

Oct 12, 2004, 2:16:30â€¯AM10/12/04

to

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.

> 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

Oct 12, 2004, 2:02:14â€¯AM10/12/04

to

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.

>

Oct 14, 2004, 7:11:50â€¯AM10/14/04

to

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}]].

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

>

Oct 14, 2004, 7:10:49â€¯AM10/14/04

to

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

Oct 14, 2004, 7:20:00â€¯AM10/14/04

to

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:

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*

Oct 14, 2004, 7:09:47â€¯AM10/14/04

to

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

you all.

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*

Oct 14, 2004, 7:06:44â€¯AM10/14/04

to

Oct 14, 2004, 6:56:35â€¯AM10/14/04

to

Oct 14, 2004, 7:03:41â€¯AM10/14/04

to

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)

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:

Oct 15, 2004, 3:45:21â€¯AM10/15/04

to

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

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

>

Oct 15, 2004, 3:44:20â€¯AM10/15/04

to

Reply all

Reply to author

Forward

0 new messages