home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / modula2 / library / queuem2 / normrng.mod < prev    next >
Text File  |  1989-08-02  |  6KB  |  175 lines
  1. (* source: h:\modula\code\random\NormRNG.MOD    v1.0b       revised: 88.06.09
  2.    author: G.Greene, AGCS D/429 (NS/TS), 312/681-7783       created: 88.06.09
  3.  
  4.    function:
  5.     This is the implementation code for a package which provides generators
  6.     of pseudo-random variates of the normal and related distributions in the
  7.     form of double-precision reals.
  8.  
  9.    history:
  10.     88.06.09  1.0a  initial release.
  11. *)
  12.  
  13. IMPLEMENTATION MODULE  NormRNG;
  14.  
  15.  
  16. FROM  MATHLIB   IMPORT (*PROC*) Log, Exp, Sqrt;
  17.  
  18. FROM  UnifRNG   IMPORT (*PROC*) UniformVariate;
  19.  
  20.  
  21. VAR
  22.   HaveSpare:  BOOLEAN;
  23.   SpareValue: LONGREAL;
  24.  
  25.  
  26. (*  Return as the function value a normally-distributed LongReal value, with
  27.     mean=0 and standard deviation=1.  If this function has previously been
  28.     invoked an even number of times (including zero), two values will be
  29.     generated.  One value is returned;  the other is saved for future use in
  30.     a static variable, local to this package (not exported).  On alternate
  31.     invocations, the saved value will be returned.  The local static boolean
  32.     variable HaveSpare is toggled to indicate the source of the next value.
  33.  
  34.     The algorithm used is "algorithm P" in section 3.4.1 of "The Art of
  35.     Computer Programming: vol II (Knuth, 1969).  Effectively, we start by
  36.     generating points (x, y) within a square with coordinates in (-1, +1).
  37.     Acceptable points for normal variate generation are those within the unit
  38.     circle [center at (0, 0), radius=1).  Thus the condition S<1.0 will hold
  39.     pi/4 (>78.5%) of the time.  More than 2 REPEAT loop repetitions will be
  40.     required less than 5% of the time.
  41. *)
  42.  
  43. PROCEDURE  StdNormalVariate ( ): LONGREAL;
  44.  
  45.   VAR
  46.     X, Y,                (* coordinates of point *)
  47.     S,                   (* distance of point (x, y) from (0, 0) *)
  48.     SFactor: LONGREAL;   (* common factor for the two random variates *)
  49.  
  50. (*                                                                         [2]
  51.  source: h:\modula\code\random\NormRNG.MOD    v1.0b       revised: 88.06.09 *)
  52.  
  53.  
  54.   BEGIN   (* StdNormalVariate *)
  55.     IF  HaveSpare  THEN
  56.       HaveSpare := FALSE;
  57.       RETURN ( SpareValue );
  58.  
  59.     ELSE
  60.       REPEAT
  61.         X := 2.0 * UniformVariate ( )  -  1.0;
  62.         Y := 2.0 * UniformVariate ( )  -  1.0;
  63.         S := X * X  +  Y * Y;
  64.       UNTIL  S < 1.0;
  65.  
  66.       (* NB:  S>0, since X and Y cannot both be zero;  ln(S)<0, since S<1 *)
  67.       SFactor    := Sqrt ( -2.0 * Log ( S ) / S );
  68.       SpareValue := X * SFactor;
  69.       HaveSpare  := TRUE;
  70.  
  71.       RETURN  Y * SFactor;
  72.     END;
  73.   END  StdNormalVariate;
  74.  
  75.  
  76.  
  77. (*  Return as the function value a normally-distributed real value, with
  78.     mean given by the first parameter, and standard deviation given by the
  79.     second parameter.
  80. *)
  81.  
  82. PROCEDURE  NormalVariate (
  83.                   (*in*)  Mean:   REAL;
  84.                   (*in*)  StdDev: REAL ): LONGREAL;
  85.  
  86.   BEGIN
  87.     RETURN  LONGREAL ( Mean )  +  LONGREAL ( StdDev ) * StdNormalVariate ( )
  88.   END  NormalVariate;
  89.  
  90.  
  91.  
  92. (*  Return as the function value a lognormally-distributed real value, with
  93.     location and shape parameters given by the function parameters.  More
  94.     commonly, one would want to specify the mean and standard deviation of
  95.     the values.  Since the computations involved in providing that form
  96.     for the parameters would substantially slow the process of generating
  97.     variates, a separate procedure, LognormalParameters, is provided below
  98.     to do the conversion once for a given distribution.  The output values
  99.     from that procedure may then be used as the parameters for this one.
  100. *)
  101.  
  102. PROCEDURE  LognormalVariate (
  103.                      (*in*)  Location: REAL;
  104.                      (*in*)  Shape:    REAL ): LONGREAL;
  105.  
  106.   BEGIN
  107.     RETURN  Exp ( LONGREAL ( Location )  +
  108.                   LONGREAL ( Shape ) * StdNormalVariate ( ) );
  109.   END  LognormalVariate;
  110.  
  111. (*                                                                         [3]
  112.  source: h:\modula\code\random\NormRNG.MOD    v1.0b       revised: 88.06.09 *)
  113.  
  114.  
  115. (*  The expected value (variate mean) and standard deviation of values for a
  116.     lognormal distribution are given as the first and second parameters,
  117.     respectively.  Set the third and fourth parameters to the corresponding
  118.     location and shape parameters for use with procedure LognormalVariate,
  119.     above.
  120. *)
  121.  
  122. PROCEDURE  LognormalParameters (
  123.                    (*in*)       Mean:     REAL;
  124.                    (*in*)       StdDev:   REAL;
  125.                   (*out*)  VAR  Location: REAL;
  126.                   (*out*)  VAR  Shape:    REAL );
  127.  
  128.   VAR
  129.     Variance: REAL;
  130.  
  131.   BEGIN
  132.     Variance := REAL ( Log ( 1.0  +
  133.                              LONGREAL ( StdDev * StdDev ) /
  134.                                                 LONGREAL ( Mean * Mean ) ) );
  135.     Location := REAL ( Log ( LONGREAL ( Mean ) ) )  -  Variance * 0.5;
  136.     Shape    := REAL ( Sqrt ( LONGREAL ( Variance ) ) );
  137.   END  LognormalParameters;
  138.  
  139.  
  140.  
  141. (*  Return as the function value a Cauchy-distributed real value, with
  142.     location (median) and scale parameters given by the function parameters.
  143.     The implementation assumes that the StdNormalVariate function will not
  144.     generate two successive zero values.  For the current implementation of
  145.     StdNormalVariate, this is equivalent to saying that UniformVariate will
  146.     not generate consecutive zero values.
  147. *)
  148.  
  149. PROCEDURE  CauchyVariate (
  150.                   (*in*)  Median: REAL;
  151.                   (*in*)  Scale:  REAL ): LONGREAL;
  152.  
  153.   VAR
  154.     Normal1,
  155.     Normal2,
  156.     StdCauchy: LONGREAL;
  157.  
  158.   BEGIN
  159.     Normal1 := StdNormalVariate ( );
  160.     Normal2 := StdNormalVariate ( );
  161.  
  162.     IF  Normal1 = 0.0
  163.       THEN  StdCauchy := Normal1 / Normal2
  164.       ELSE  StdCauchy := Normal2 / Normal1;
  165.     END;
  166.  
  167.     RETURN  LONGREAL ( Median )  +  StdCauchy * LONGREAL ( Scale );
  168.   END  CauchyVariate;
  169.  
  170.  
  171. BEGIN  (* initialization *)
  172.   HaveSpare  := FALSE;
  173.   SpareValue := 0.0;
  174. END  NormRNG.
  175.