IMPLEMENTATION MODULE RAND1; (******************************************************************************* Environnement : DOS et TopSpeed Modula-2. ------------- Programmeurs : ------------ Pierre L'Ecuyer, Pierre Poulin, Gaetan Perron, Jean Belanger, Denis Alain *******************************************************************************) FROM Storage IMPORT ALLOCATE; FROM NUM IMPORT Ln, Pow, Sqrt, Sin, Cos, Exp; FROM MYINOUT IMPORT WriteString, WriteLn; FROM RAND IMPORT Gen, Random, Uniform01; FROM SYSTEM IMPORT ADDRESS; (* Ref. : L'Ecuyer, P. et Cote, S. (1987). "A Random Number Package with *) (* Splitting Facilities." Rapport no. DIUL-RR-8705, *) (* Departement d'informatique , Universite Laval. *) CONST Pi = 3.1415926; MaxIndex = 8190; TYPE DiscreteDist = POINTER TO RECORD Nb : CARDINAL; Value, Function : POINTER TO ARRAY [0..MaxIndex] OF LONGREAL END; VAR MemGamma : ARRAY [1..6] OF LONGREAL; K : LONGCARD; PSave, QSave : LONGREAL; Con : ARRAY [1..3] OF LONGREAL; pp, qq, x2, x3, x4, f2, s, r, f4, xl2, xl4, p1, p2, p3 : LONGREAL; (******************************************************************************) PROCEDURE Error (Number : CARDINAL); BEGIN WriteString ( '***** ERROR from module RAND1 : ***** ' ); WriteLn; CASE Number OF 1 : WriteString ("Calling procedure Discrete with a discrete dist."); WriteString (" which has not been created yet.") | 2 : WriteString ("When calling CreateDiscreteDist, the values in the"); WriteString (" array V must be in increasing order."); | 3 : WriteString ("When calling CreateDiscreteDist, the values in the"); WriteString (" array F must be in increasing order."); | 4 : WriteString ("When calling CreateDiscreteDist, the last value"); WriteString (" in the array F must be 1.0."); | 5 : WriteString ("Calling Beta with a negative parameter"); | 6 : WriteString ("Calling Weibull with a negative parameter"); | 7 : WriteString ("Calling Gamma with a negative parameter"); | 8 : WriteString ("Calling Normal with a negative standard deviation"); | 9 : WriteString ("Calling InvStudentDist with u < E-20 or u > 1 - E-20"); | 10 : WriteString ("Calling InvStudentDist with n < 1"); | 11 : WriteString ("Calling Student with DegFreedom < 1") END; WriteLn; WriteLn; HALT END Error; (******************************************************************************) PROCEDURE DiscreteUniform ( N1, N2 : LONGINT ; g : Gen ) : LONGINT ; (* Condition pour la validite de la fonction 0 <= IX < P Valeur aleatoire, IX = Random (g). 0 < L < P Longueur de l'intervalle, L = N2 - N1 + 1. 1 < P < 2**31 Valeur du modulo utilise pour obtenir IX Elle retourne (IX * L) DIV P + N1. Ref : A Guide to Simulation, Second Edition, by Bratley, Fox et Schrage. Edition Springer-Verlag. pp. 221--223. *) CONST M = 2147483647; P = 2147483563; P15 = 32768; P16 = 65536; VAR Nonfini : BOOLEAN; IX, X, Y, U, V, YV1, YV2, R, S, UR1, UR2, VX1, VX2, A, B, AB, K, L : LONGINT; BEGIN IX := Random (g); L := ABS (N2 - N1 + 1); X := L DIV P16; Y := L - P16 * X; U := IX DIV P15; V := IX - P15 * U; K := Y * V; YV1 := K DIV P15; YV2 := K - P15 * YV1; R := Y DIV 2; S := Y - 2 * R; K := U * R; UR1 := K DIV P15; UR2 := K - P15 * UR1; IF X > 0 THEN K := V * X; VX1 := K DIV P15; VX2 := K - P15 * VX1; VX1 := VX1 + U * X ELSE VX1 := 0; VX2 := 0 END; AB := YV1 + 2 * ( UR2 + VX2 ) + S * U; K := AB DIV P16; B := P15 * ( AB - P16 * K ) + YV2; A := K + UR1 + VX1; R := M - P + 1; S := M DIV R; U := A; K := 0; V := B - P; Nonfini := TRUE; WHILE Nonfini DO WHILE V >= 0 DO V := V - P; K := K + 1 END; IF U = 0 THEN Nonfini := FALSE ELSE IF S < U THEN X := S ELSE X := U END; U := U - X; V := V + R * X END END; RETURN A + K + N1 END DiscreteUniform; (******************************************************************************) PROCEDURE CreateDiscreteDist ( VAR D : DiscreteDist; n : CARDINAL; VAR V, F : ARRAY OF LONGREAL ); VAR Addr : ADDRESS; I : CARDINAL; BEGIN NEW (D); WITH D^ DO IF F[n-1] # 1.0 THEN Error (4) END; Nb := n; ALLOCATE (Addr, 8 * n); Value := Addr; Value^[0] := V[0]; FOR I := 1 TO n - 1 DO IF V[I] > V [I-1] THEN Value^[I] := V[I] ELSE Error (2) END END; ALLOCATE (Addr, 8 * n); Function := Addr; Function^[0] := F[0]; FOR I := 1 TO n - 1 DO IF F[I] > F[I-1] THEN Function^[I] := F[I] ELSE Error (3) END END END END CreateDiscreteDist; (******************************************************************************) PROCEDURE Discrete ( D : DiscreteDist; g : Gen ) : LONGREAL; VAR U : LONGREAL; Low, High, Dif, DifDiv2 : CARDINAL; BEGIN U := Uniform01 (g); IF D = NIL THEN Error (1) END; WITH D^ DO IF Function^[0] < U THEN Low := 0; High := Nb - 1; Dif := High - Low; WHILE Dif > 1 DO DifDiv2 := ( Dif DIV 2 ) + Low; IF Function^[DifDiv2] < U THEN Low := DifDiv2 ELSE High := DifDiv2 END; Dif := High - Low END ELSE High := 0 END; RETURN Value^[High] END END Discrete; (******************************************************************************) PROCEDURE Uniform ( A, B : LONGREAL; g : Gen ): LONGREAL; BEGIN RETURN A + ((B - A) * Uniform01 (g)) END Uniform; (******************************************************************************) PROCEDURE Expon ( Mean : LONGREAL; g : Gen ): LONGREAL; BEGIN RETURN -Mean * Ln (1.- Uniform01 (g)) END Expon; (******************************************************************************) PROCEDURE Weibull ( Alpha : LONGREAL; Lambda : LONGREAL; g : Gen ): LONGREAL; BEGIN IF NOT ((Alpha > 0.) AND (Lambda > 0.)) THEN Error (6) END; RETURN Pow ((-Ln (1.0 - Uniform01 (g)) / Lambda ), (1.0 / Alpha)) END Weibull; (******************************************************************************) PROCEDURE Gamma ( Alpha, Lambda : LONGREAL; g : Gen ): LONGREAL; (* Ref. : Bratley, P., Fox, B.L., Schrage, L.E. (1987). *) (* "A Guide to Simulation", Springer-Verlag. *) VAR U1, U2, U : LONGREAL; B, A, P, X, W : LONGREAL; Flag : BOOLEAN; BEGIN IF NOT ((Alpha > 0.) AND (Lambda > 0.)) THEN Error (6) END; Flag := TRUE; IF Alpha <= 1.0 THEN WHILE Flag DO U1 := Uniform01 (g); B := ( 2.718281828 + Alpha) / 2.718281828; P := B * U1; U2 := Uniform01(g); IF P > 1.0 THEN X := -Ln ((B - P) / Alpha); Flag := Ln (U2) > (Alpha - 1.0) * Ln (X) ELSE X := Exp (Ln (P) / Alpha); Flag := U2 > Exp (-X) END END ELSE IF (MemGamma[1] # Alpha) THEN MemGamma[1] := Alpha; K := 1; IF (Alpha > 2.5) THEN K := 2 END; MemGamma[2] := Alpha - 1.0; MemGamma[3] := (Alpha - 1.0 / (6.0 * Alpha)) / MemGamma[2]; MemGamma[4] := 2.0 / MemGamma[2]; MemGamma[5] := MemGamma[4] + 2.0; MemGamma[6] := Sqrt (Alpha) END; WHILE Flag DO IF K = 1 THEN U1 := Uniform01(g); U2 := Uniform01(g) ELSE REPEAT U1 := Uniform01(g); U := Uniform01(g); U2 := U1 + (1.0 - 1.86 * U) / MemGamma[6] UNTIL (U2 > 0.0) AND (U2 < 1.0) END; W := MemGamma[3] * U1 / U2; IF ((MemGamma[4] * U2 - MemGamma[5] + W + 1.0 / W) <= 0.0) THEN Flag := FALSE ELSE Flag := ((MemGamma[4] * Ln(U2) - Ln(W) + W - 1.0) >= 0.0) END END; X := MemGamma[2] * W END; RETURN X / Lambda END Gamma; (******************************************************************************) PROCEDURE Beta ( A, B : LONGREAL; g : Gen ): LONGREAL; (* Ref. : Bratley, P., Fox, B.L., Schrage, L.E. (1987). *) (* "A Guide to Simulation", Springer-Verlag. *) VAR p, q, d, u, u1, u2, v, x, a : LONGREAL; Flag, Flag2 : BOOLEAN; (* PROCEDURE Generate ( g : Gen ); BEGIN u := Uniform01 (g) * p3; v := Uniform01 (g); IF Antithetic[g] THEN v := 1.0 - v END END Generate; *) BEGIN IF NOT ((A > 0.) AND (B > 0.)) THEN Error (5) END; IF (Con[1] <= 0.0) THEN IF (A < B) THEN Con[1] := A ELSE Con[1] := B END; IF (Con[1] <= 1.0) THEN Con[1] := 1.0 / Con[1] ELSE Con[1] := Sqrt ((A + B - 2.0) / (2.0 * A * B - A - B)) END; Con[2] := A + B; Con[3] := A + 1.0 / Con[1] END; REPEAT u1 := Uniform01(g); u2 := Uniform01(g); v := Con[1] * Ln(u1 / (1.0 - u1)); x := A * Exp (v) UNTIL ((Con[2] * Ln(Con[2] / (B + x)) + Con[3] * v - 1.3862944) >= (Ln(u1 * u1 * u2))); x := x / (B + x); (* Apres quelques repetitions de ce programme, dont le nombre depend des valeurs de A, B, g1 et g2, il semble entrer dans une boucle infinie dans le WHILE qui est plus loin. ELSE IF (A # PSave) OR (B # QSave) THEN PSave := A; QSave := B; IF (A <= B) THEN p := A; q := B ELSE p := B; q := A END; pp := p - 1.0; qq := q - 1.0; r := pp + qq; s := r * Log(r); x2 := 0.0; f2 := 0.0; f4 := 0.0; x3 := pp / r; x4 := 1.0; IF (r > 1.0) THEN d := sqrt(pp * qq / (r - 1.0)) / r; IF (d < x3) THEN x2 := x3 - d; xl2 := pp / x2 - qq / (1.0 - x2); f2 := exp(pp * Log(x2 / pp) + qq * Log((1.0 - x2) / qq) + s) END; IF ((x3 + d) < 1.0) THEN x4 := x3 + d; xl4 := qq / (1.0 - x4) - pp / x4; f4 := exp(pp * Log(x4 / pp) + qq * Log((1.0 - x4) / qq) + s) END END; p2 := 0.0; IF (xl2 > 0.0) THEN p2 := f2 / xl2 END; p1 := (x4 - x2) + p2; p3 := p1; IF (xl4 > 0.0) THEN p3 := f4 / xl4 + p1 END END; Flag := TRUE; Flag2 := TRUE; Generate (g); WHILE Flag DO (* Boucle mentionnee plus haut *) IF (u < p2) THEN u := u / p2; x := x2 + Log(u) / xl2; IF (v < ((xl2 * (x - x2) + 1.0) / u)) THEN Flag := FALSE; Flag2 := FALSE ELSE IF (x > 0.0) THEN v := v * f2 * u ELSE Flag2 := FALSE; Generate (g2) END END ELSE IF (u <= p1) THEN x := x2 + (u - p2); IF ((x < x3) AND (v < (f2+(x-x2) * (1.0- f2) / (x3 - x2)))) THEN Flag := FALSE; Flag2 := FALSE ELSE IF ((x >= x3) AND (v < (f4+(x4-x)*(1.0- f4)/(x4 - x3)))) THEN Flag := FALSE; Flag2 := FALSE END END ELSE u := (u - p1) / (p3 - p1); x := x4 - Log(1.0 - u) / xl4; IF (v < ((xl4 * (x4 - x) + 1.0) / (1.0 - u))) THEN Flag := FALSE; Flag2 := FALSE ELSE IF (x >= 1.0) THEN Flag2 := FALSE; Generate (g2) ELSE v := v * f4 * (1.0 - u) END END END END; IF Flag2 THEN a := Log(v); IF (a > (-(x - x3) * (x - x3) * (r + r))) THEN Generate (g2) ELSE IF (a > (pp * Log(x / pp) + qq * Log((1.0 - x) / qq) + s)) THEN Generate (g2) END END END END; IF (A > B) THEN x := 1.0 - x END END; *) RETURN x END Beta; (******************************************************************************) PROCEDURE InvNormalDist ( U : LONGREAL ): LONGREAL; (* Utilise l'inversion et une approximation rationnelle donnant environ 7 *) (* decimales de precision pour 1.0E-20 < U < 1.0 - 1.0E-20. *) (* Ref. : Kennedy and Gentle, "Statistical Computing", Dekker, 1980, p.95. *) CONST P0 = -0.322232431088; Q0 = 0.0993484626060; P1 = -1.000000000000; Q1 = 0.588581570495; P2 = -0.342242088547; Q2 = 0.531103462366; P3 = -0.0204231210245; Q3 = 0.103537752850; P4 = -0.0000453642210148; Q4 = 0.0038560700634; VAR Y, Z : LONGREAL; BEGIN IF U > 0.5 THEN Y := Sqrt ( -2.0 * Ln ( 1.0 - U ) ) ELSE Y := Sqrt ( -2.0 * Ln ( U ) ) END; Z := Y + (((( Y * P4 + P3 ) * Y + P2 ) * Y + P1 ) * Y + P0 ) / (((( Y * Q4 + Q3 ) * Y + Q2 ) * Y + Q1 ) * Y + Q0 ); IF U < 0.5 THEN Z := -Z END; RETURN Z END InvNormalDist; (******************************************************************************) PROCEDURE Normal ( Mean : LONGREAL; StdDeviation : LONGREAL; g : Gen ): LONGREAL; BEGIN IF StdDeviation < 0.0 THEN Error (8) END; RETURN Mean + InvNormalDist (Uniform01 (g)) * StdDeviation; END Normal; (******************************************************************************) PROCEDURE InvStudentDist ( n : LONGINT; u : LONGREAL ): LONGREAL; (* Cette fonction retourne la valeur de t telle que F(t) = u, ou F est *) (* la fonction de repartition d'une loi de Student a n degres de liberte. *) (* Ref. : Hill, G. W., (1970b), Algorithme 396: Student's t-Quantiles, *) (* CACM 13, pp.619-620 *) CONST Limite = 1.E-20; (* On doit avoir Limite < u < 1 - Limite. *) VAR a, b, c, d, e, p, t, x, y : LONGREAL; BEGIN e := FLOAT (n); IF u > 0.5 THEN p := 2. * (1. - u) ELSE p := 2. * u END; IF n < 1 THEN Error (10) END; IF p <= 2.0 * Limite THEN Error (9) END; IF n = 1 THEN t := ABS ( Cos ( Pi * p / 2. ) /Sin ( Pi * p / 2. ) ) ELSIF n = 2 THEN t := Sqrt ( 2. / ( p * ( 2. - p ) ) - 2. ) ELSE (* n > 2 *) a := 1. / ( e - 0.5 ); b := 48. / ( a * a ); c := ( ( 20700. / b * a - 98. ) * a - 16. ) * a + 96.36; d := e * Sqrt (a * Pi / 2.) * ( ( 94.5 / ( b + c ) - 3. ) / b + 1. ); y := Pow (( d * p ), ( 2. / e )); IF y > a + 0.05 THEN IF p = 1. THEN x := 0. (* Permet de sauver un appel a InvNormalDist. *) ELSE x := InvNormalDist ( p * 0.5 ) END; y := x * x; IF n < 5 THEN c := c + 0.3 * ( e - 4.5 ) * ( x + 0.6 ) END; c := ( ( ( 0.05 * d * x - 5. ) * x - 7. ) * x - 2. ) * x + b + c; y := ( ( ( ( ( 0.4 * y + 6.3 ) * y + 36. ) * y + 94.5 ) / c - y - 3. ) / b + 1. ) * x; y := a * ( y * y ); IF y > 0.002 THEN y := Exp ( y ) - 1. ELSE y := ( 0.5 * y * y ) + y END ELSE y := ( ( 1. / ( ( ( e + 6. ) / ( e * y ) - 0.089 * d - 0.822 ) * ( e + 2. ) * 3. ) + 0.5 / ( e + 4. ) ) * y - 1. ) * ( e + 1. ) / ( e + 2. ) + 1. / y END; t := Sqrt ( e * y ) END; IF u < 0.5 THEN RETURN - t ELSE RETURN t END END InvStudentDist; (******************************************************************************) PROCEDURE Student ( DegFreedom : LONGINT; g : Gen ): LONGREAL; VAR U : LONGREAL; BEGIN IF DegFreedom < 1 THEN Error (11) END; U := Uniform01 (g); RETURN InvStudentDist (DegFreedom, U) END Student; (******************************************************************************) BEGIN (* Initialisation *) MemGamma [1] := -1.0; PSave := -1.0; QSave := -1.0; Con [1] := 0.0 END RAND1.