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

Re: strange problems with Random

• To: mathgroup at smc.vnet.net
• Subject: [mg47937] Re: [mg47812] strange problems with Random
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Sat, 1 May 2004 02:43:03 -0400 (EDT)
• References: <200404270848.EAA18923@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

After investigating this matter further (I would like to thank several 
persons who sent helpful comments, especially Jens Kuska who also 
performed some very useful tests) I have come to the conclusion that 
the problem is indeed with the subtract-with borrow algorithm that 
Daniel Lichtblau described in 
<http://forums.wolfram.com/mathgroup/archive/2000/May/msg00088.html>. 
Using the workaround he suggests (which involves switching from the 
subtract with borrow algorithm to Wolfram's CA algorithm) indeed cures 
the problem without any noticeable loss of performance. Given the 
importance of this issue (for example, using Monte Carlo methods with 
Mathematica, without the adjustment gives significantly wrong answers 
for prices of financial derivatives ) and the ease with which it can be 
corrected it is really surprising that it has not been done already 
several years after the problem was  first discovered. But this time, 
the problem seems to reveal another curious aspect.
Of course I knew about the problem with Random[] when I observed that 
Random[NormalDistribution[...]] was giving wrong answers, but I 
believed that this could not be the issue here because I could obtain 
correct answers by replacing Random[NormalDistribution[...]] with an 
alternative, which I called RandomNormal in my earlier posting, which 
itself was using Random[]. So if the problem was with Random[] it would 
seem strange that a function defined in terms of Random[] and 
mathematically equivalent to  Random[NormalDistribution[...]] should be 
giving correct answers. Yet it does  seem to do so consistently. I have 
since considered this issue more carefully and now I seem to have a 
plausible explanation, but not being a statistician I would like to 
hear if people who know more about this field would agree with it. The 
point is this. To generate a normally distributed random variable from 
the uniformly distributed random variable Random[] built-in into 
Mathematica, the package Statistics`NormalDistribution` uses the well 
known Box-Muller method. This method begins by choosing a random point 
in the unit square using the given uniformly distributed random 
variable (Random[]). The generated normally distributed random variable 
(actually a pair of random variables) is expressed in terms of 
trigonometric functions. In order to avoid this, and thus improve 
performance (computing trigonometric functions is a relatively slow 
process) Marsaglia introduced another method, in which one samples 
repeatedly points in the unit square until one obtains a point in the 
unit circle. My function RandomNormal was an implementation of the 
Marsaglia method.  In the case of this method usually the gain in 
performance due to the fact that one no longer needs to compute trig 
functions outweighs the loss due to the necessity to re-sample in 
unlucky cases. However, it appears that the method has another 
advantage: it appears to eliminate the "non-randomness" in the 
algorithm in Random[]. It would seem as if the "additional layer" of 
randomness, due to the occasional repeated sampling, needed in the 
Marsaglia method, was just enough to offset the insufficient randomness 
in Random[].At least, that is my current hypothesis. I anybody has any 
thoughts on this I would like to hear them.
(What is still very puzzling, is the problem reported my Mark Coleman. 
In his case Random[UniformDistribution[]] was giving incorrect answers 
while Random[] worked fine! )

Andrzej Kozlowksi


On 27 Apr 2004, at 17:48, Andrzej Kozlowski wrote:

> Recently Mark Coleman wrote a message with the subject "Differences in
> Random Numbers" in which he described his experience of getting
> statistically incorrect answers when using
> Random[UniformDistribution[0,1]] instead of simple Random[]. Jens Kuska
> and Bill Rowe pointed out that the package  ContinuousDistributions
> contains the code
>
> UniformDistribution/: Random[UniformDistribution[min_:0, max_:1]] :=
> Random[Real, {min, max}]
>
> and that seemed to be it. However, I have now encountered a problem
> that makes me feel strongly that there is more to this problem than it
> seemed to me at first.
>
> I have been doing some Monte Carlo computations of values of of
> financial options for the Internet based course on this subject I am
> teaching for Warsaw University.  In the case of the simplest options,
> the so called vanilla european options, an exact formula (due to black
> and Scholes) is known. Neverheless, computations with Mathematica using
> the function
>
> Random[NormalDistribution[mu,sigma]] (for some fixed mu and sigma) have
> been giving answers which are consistently higher than the
> Black-Scholes value. On the other hand, when
> Random[NormalDistribution[mu,sigma]]  is replaced by the following code
> using a well known algorithm due to Marsaglia:
>
> RandomNormal=Compile[{mu, sigma}, Module[{va = 1., vb, rad =
>        2.0, den = 1.}, While[rad ³ 1.00, (va = 2.0*Random[] - 1.0;
>              vb = 2.0*Random[] - 1.0;
>              rad = va*va + vb*vb)];
>          den = Sqrt[-2.0*Log[rad]/rad];
>          mu + sigma*va*den]];
>
> the answers agree very closely with the Black-Scholes ones.
> I am not sure if this problem is related to the one mark mentioned but
> it seems quite likely.  One could in principle test it by replacing
> Random[] by Random[Uniformdistribution[0,1]] (after loading the
> Statistics`ContinuousDistributions`) package, but if one does that the
> code will no longer compile. That makes the function to slow for tests
> that are sufficiently precise, so I have not tested it.
>
> In any case, I am sure that using Random[Normal  ] consistently gives
> wrong answers. I looked at the Statistics`NormalDistributionPackage`
> but everything looks O.K. to me. I now tend to believe that the problem
> lies with Random itself, that is that  Random called with any
> parameters or perhaps some particular parameters may not be behaving
> correctly. Random[] called with no parameters does not seem to  suffer
> from any such problems.
>
>
> Andrzej Kozlowski
> Chiba, Japan
> http://www.mimuw.edu.pl/~akoz/
>
>