[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/
>
>
Prev by Date:
**RE: on ColorFunction and combined images?**
Next by Date:
**Re: Re: bug in IntegerPart ?**
Previous by thread:
**Re: on ColorFunction and combined images?**
Next by thread:
**Re: Re: bug in IntegerPart ?**
| |