home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / gnu / g / lib / bug / 739 < prev    next >
Encoding:
Text File  |  1993-01-05  |  2.4 KB  |  91 lines
  1. Newsgroups: gnu.g++.lib.bug
  2. Path: sparky!uunet!zaphod.mps.ohio-state.edu!cis.ohio-state.edu!pic.ucla.edu!jimc
  3. From: jimc@pic.ucla.edu
  4. Subject: Normal.cc may not be normal
  5. Message-ID: <9301050046.AA06293@malibu.pic.ucla.edu>
  6. Sender: gnulists@ai.mit.edu
  7. Organization: GNUs Not Usenet
  8. Distribution: gnu
  9. Date: Mon, 4 Jan 1993 08:46:18 GMT
  10. Approved: bug-lib-g++@prep.ai.mit.edu
  11. Lines: 78
  12.  
  13. Gnu g++ v2.3.2 for Sparc-Sun-SunOS 4.1.1.  libg++ v2.3.  
  14.  
  15. In Normal.cc, I am not sure that the given algorithm gives an exactly
  16. normally distributed pair of results.  I have included a modification
  17. implementing the Box-Muller algorithm, which is exact, at the cost of
  18. a sin/cos evaluation.  
  19.  
  20. James F. Carter        (213) 825-2897
  21. UCLA-Mathnet;  6221 MSA; 405 Hilgard Ave.; Los Angeles, CA, USA  90024-1555
  22. Internet: jimc@math.ucla.edu            BITNET: jimc%math.ucla.edu@INTERBIT
  23. UUCP:...!{ucsd,ames,ncar,gatech,purdue,rutgers,decvax,uunet}!math.ucla.edu!jimc
  24.  
  25. ----------- From libg++.2.3/libg++/src/Normal.cc
  26.  
  27. double Normal::operator()()
  28. {
  29.     
  30.     if (haveCachedNormal == 1) {
  31.     haveCachedNormal = 0;
  32.     return(cachedNormal * pStdDev + pMean );
  33.     } else {
  34.     
  35.     for(;;) {
  36.         double u1 = pGenerator -> asDouble();
  37.         double u2 = pGenerator -> asDouble();
  38.         double v1 = 2 * u1 - 1;
  39.         double v2 = 2 * u2 - 1;
  40.         double w = (v1 * v1) + (v2 * v2);
  41.         
  42. //
  43. //    We actually generate two IID normal distribution variables.
  44. //    We cache the one & return the other.
  45. // 
  46.         if (w <= 1) {
  47.         double y = sqrt( (-2 * log(w)) / w);
  48.         double x1 = v1 * y;
  49.         double x2 = v2 * y;
  50.         
  51.         haveCachedNormal = 1;
  52.         cachedNormal = x2;
  53.         return(x1 * pStdDev + pMean);
  54.         }
  55.     }
  56.     }
  57. }
  58.  
  59. ------------ Suggested change
  60.  
  61. // Reference: Ripley, Brian D.  Stochastic Simulation.  Wiley (NY) 1987.
  62. // p. 54, algorithm 3.1, Box-Muller (not Mu\"ller).
  63. double Normal::operator()()
  64. {
  65.     
  66.     if (haveCachedNormal == 1) {
  67.     haveCachedNormal = 0;
  68.     return(cachedNormal * pStdDev + pMean );
  69.     } else {
  70.  
  71. //    Get 2 flat random variables in 0..1.  Generator gives 0 <= u1 < 1.
  72.     double u1 = pGenerator -> asDouble();
  73. //    This one must never be zero.
  74.     double u2;
  75.     while (0.0 == (u2 = pGenerator -> asDouble())) /*null stmt*/;
  76.         
  77. //
  78. //    We actually generate two IID normal distribution variables.
  79. //    We cache the one & return the other.
  80. // 
  81.     double th = 6.28318530717958648 * u1
  82.     double r = sqrt(-2 * log(u2));
  83.         
  84.     haveCachedNormal = 1;
  85.     cachedNormal = r * sin(th);
  86.     return(pStdDev*r*cos(th) + pMean);
  87.     }
  88. }
  89.  
  90.  
  91.