home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #16 / NN_1992_16.iso / spool / sci / math / stat / 1529 < prev    next >
Encoding:
Text File  |  1992-07-27  |  4.5 KB  |  141 lines
  1. Path: sparky!uunet!cs.utexas.edu!sdd.hp.com!usc!news
  2. From: ajayshah@almaak.usc.edu (Ajay Shah)
  3. Newsgroups: sci.math.stat
  4. Subject: Re: random number generator
  5. Date: 27 Jul 1992 20:59:51 -0700
  6. Organization: University of Southern California, Los Angeles, CA
  7. Lines: 128
  8. Sender: ajayshah@almaak.usc.edu (Ajay Shah)
  9. Message-ID: <l79hhnINNf0k@almaak.usc.edu>
  10. References: <1992Jul27.132508.217@ginger.hnrc.tufts.edu> <1992Jul27.201230.984@organpipe.uug.arizona.edu>
  11. NNTP-Posting-Host: almaak.usc.edu
  12.  
  13. paul@music.sie.arizona.edu (Paul J. Sanchez) writes:
  14.  
  15. >There are a couple of morals to this story:
  16. >a) If you want to use Box-Muller, use something other than a LCG; and
  17.  
  18. Or, if you want to use a LCG then don't use the Box-Muller?
  19.  
  20. Could someone comment on the algorithm in Algorithm 488, CACM, Vol.
  21. 17, #12, p704, which draws from N(0,1).  
  22.  
  23. C source for it is here:
  24.  
  25. #define REAL float
  26. /* Algorithm 488, CACM, Vol. 17, #12, p704.
  27.  
  28. Except on the first call ranstdn2 returns a pseudo-random number
  29. having a gaussian (i.e.  normal) distribution with zero mean and unit
  30. standard deviation.
  31. Thus, the density is f(x) = exp(-0.5*x**2)/sqrt(2.0*pi).
  32. The first call initializes ranstdn2 and returns zero
  33. (no hassles -- it's a legit draw from N(0,1)).
  34.  
  35. It calls a function ranpm(), and it is assumed that successive calls
  36. to ranpm(0) give independent pseudo- random numbers distributed
  37. uniformly on (0, 1), possibly including 0 (but not 1).  the method
  38. used was suggested by Von Neumann, and improved by Forsythe, Ahrens,
  39. Dieter and Brent.  on the average there are 1.37746 calls of ranpm() for
  40. each call of ranstdn().
  41.  
  42. Warning - dimension and data statements below are machine-dependent.
  43. The dimension of d must be at least the number of bits in the fraction
  44. of a REAL number.  Thus, on most machines the data statement
  45. below can be truncated.
  46.  
  47. How is d calculated:
  48. If the integral of sqrt(2.0/pi)*exp(-0.5*x**2) from a(i) to infinity
  49. is 2**(-i), then d(i) = a(i) - a(i-1).   Thus it seems that we could
  50. do better than the constants in this code by calculating d[] in
  51. double precision using a modern pnorm1.
  52.  
  53. On a ELC using gcc 1.4 -O, it gives 2400 draws a CPU second.
  54. */
  55.  
  56. REAL ranstdn2(void)
  57. {
  58.     /* Initialized data */
  59.     static REAL d[60] = { .67448975f,.47585963f,.383771164f,.328611323f,
  60.         .291142827f,.263684322f,.242508452f,.225567444f,.211634166f,
  61.         .199924267f,.189910758f,.181225181f,.1736014f,.166841909f,
  62.         .160796729f,.155349717f,.150409384f,.145902577f,.141770033f,
  63.         .137963174f,.134441762f,.13117215f,.128125965f,.12527909f,
  64.         .122610883f,.12010356f,.117741707f,.115511892f,.113402349f,
  65.         .11140272f,.109503852f,.107697617f,.105976772f,.104334841f,
  66.         .102766012f,.101265052f,.099827234f,.098448282f,.097124309f,
  67.         .095851778f,.094627461f,.093448407f,.092311909f,.091215482f,
  68.         .090156838f,.089133867f,.088144619f,.087187293f,.086260215f,
  69.         .085361834f,.084490706f,.083645487f,.082824924f,.082027847f,
  70.         .081253162f,.080499844f,.079766932f,.079053527f,.078358781f,
  71.         .077681899f };
  72.     static REAL u = 0;
  73.  
  74.     /* Local variables */
  75.     static REAL a;
  76.     static int i;
  77.     static REAL v, w;
  78.  
  79. /* U MUST BE PRESERVED BETWEEN CALLS. */
  80. /* INITIALIZE DISPLACEMENT A AND COUNTER I. */
  81.     a = 0;
  82.     i = 0;
  83. /* INCREMENT COUNTER AND DISPLACEMENT IF LEADING BIT */
  84. /* OF U IS ONE. */
  85. L10:
  86.     u += u;
  87.     if (u < 1) goto L20;
  88.     u -= 1;
  89.     ++i;
  90.     a -= d[i - 1];
  91.     goto L10;
  92. /* FORM W UNIFORM ON 0 .LE. W .LT. D(I+1) FROM U. */
  93. L20:
  94.     w = d[i] * u;
  95. /* FORM V = 0.5*((W-A)**2 - A**2). NOTE THAT 0 .LE. V .LT. LOG(2). */
  96.     v = w * (w * .5f - a);
  97. /* GENERATE NEW UNIFORM U. */
  98. L30:
  99.     u = ranpm();
  100. /* ACCEPT W AS A RANDOM SAMPLE IF V .LE. U. */
  101.     if (v <= u) goto L40;
  102. /* GENERATE RANDOM V. */
  103.     v = ranpm();
  104. /* LOOP IF U .GT. V. */
  105.     if (u > v) goto L30;
  106. /* REJECT W AND FORM A NEW UNIFORM U FROM V AND U. */
  107.     u = (v - u) / (1 - u);
  108.     goto L20;
  109. /* FORM NEW U (TO BE USED ON NEXT CALL) FROM U AND V. */
  110. L40:
  111.     u = (u - v) / (1 - v);
  112. /* USE FIRST BIT OF U FOR SIGN, RETURN NORMAL VARIATE. */
  113.     u += u;
  114.     if (u < 1) return (a-w);
  115.     else {
  116.         u -= 1; return (w-a);
  117.     }
  118. } /* ranstdn2() */
  119.  
  120. #ifdef TESTING
  121. #include <stdlib.h>
  122. int main(int argc, char **argv)
  123. {
  124.     int N, i; 
  125.     if (argc != 3) {
  126.         fprintf(stderr, "Usage: %s [12] N\n", argv[0]);
  127.         return 1;
  128.     }
  129.     N = atoi(argv[2]);
  130.     if (atoi(argv[1]) == 2) {
  131.         (void) ranstdn2(); /* swallow the 0 */
  132.         for (i=0; i<N; i++) printf("%.6g\n", ranstdn2());
  133.     } else
  134.       for (i=0; i<N; i++) printf("%.6g\n", ranstdn());
  135.     return 0;
  136. }
  137. #endif /* TESTING */
  138.  
  139. -- 
  140. Ajay Shah, (213)749-8133, ajayshah@usc.edu
  141.