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