%PROCESS OPT(TIME), MACRO; /* Marsaglia & Tsang generator for random normals & random exponentials. */ /* Translated from C by Alan Miller (amiller @ bigpond.net.au) */ /* Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating */ /* random variables', J. Statist. Software, v5(8). */ /* This is an electronic journal which can be downloaded from: */ /* http://www.jstatsoft.org/v05/i08 */ /* N.B. It is assumed that all integers are 32-bit. */ /* N.B. The value of M2 has been halved to compensate for the lack of */ /* unsigned integers in Fortran. */ /* Latest version - 1 January 2001 */ /* Translated to PL/I by R. A. Vowels, 28 November 2001. */ /* (SUBSCRIPTRANGE, SIZE): */ test_ziggurat: PROCEDURE OPTIONS (MAIN, REORDER); %DECLARE Digits CHARACTER; %Digits = 15; DECLARE (TRUE VALUE ('1'B), FALSE VALUE ('0'B) ) BIT (1); DECLARE (m1 VALUE (2147483648.0E0), m2 VALUE (2147483648.0E0), half VALUE (0.5E0) ) FLOAT (Digits); DECLARE (dn INITIAL (3.442619855899E0), tn INITIAL (3.442619855899E0), vn INITIAL (0.00991256303526217E0), q, de INITIAL (7.697117470131487E0), te INITIAL (7.697117470131487E0), ve INITIAL (0.3949659822581572E-2) ) FLOAT (Digits); DECLARE (iz, jz, jsr INITIAL (123456789), kn(0:127), ke(0:255), hz) FIXED BINARY (31) STATIC; DECLARE (wn(0:127), fn(0:127), we(0:255), fe(0:255)) FLOAT (Digits) STATIC; DECLARE initialized BIT (1) INITIAL (FALSE) STATIC; zigset: PROCEDURE (jsrseed); DECLARE jsrseed FIXED BINARY (31); DECLARE i FIXED BINARY (31); /* Set the seed */ jsr = jsrseed; /* Tables for RNOR */ q = vn*EXP(half*dn*dn); kn(0) = (dn/q)*m1; kn(1) = 0; wn(0) = q/m1; wn(127) = dn/m1; fn(0) = 1.0; fn(127) = EXP(-half*dn*dn); DO i = 126 TO 1 BY -1; dn = SQRT(-2 * LOG(vn/dn + EXP(-half*dn*dn))); kn(i+1) = (dn/tn)*m1; tn = dn; fn(i) = EXP(-half*dn*dn); wn(i) = dn/m1; END; /* Tables for REXP */ q = ve*EXP(de); ke(0) = (de/q)*m2; ke(1) = 0; we(0) = q/m2; we(255) = de/m2; fe(0) = 1.0; fe(255) = EXP(-de); DO i = 254 TO 1 BY -1; de = -LOG(ve/de + EXP(-de)); ke(i+1) = (de/te)*m2; te = de; fe(i) = EXP(-de); we(i) = de/m2; END; initialized = TRUE; RETURN; END zigset; shr3: PROCEDURE () RETURNS (FIXED BINARY (31)); /* Generate random 32-bit integers */ jz = jsr; jsr = IEOR(jsr, ISLL(jsr, 13)); jsr = IEOR(jsr, ISRL(jsr, 17)); jsr = IEOR(jsr, ISLL(jsr, 5)); (NOFIXEDOVERFLOW): RETURN (jz + jsr); END shr3; uni: PROCEDURE () RETURNS (FLOAT (Digits) ); /* Generate uniformly distributed random numbers. */ RETURN (half + 0.232830600000E-9 * shr3() ); END uni; rnor: PROCEDURE () RETURNS (FLOAT (Digits)); /* Generate random normals. */ DECLARE r VALUE (3.44262000000E0) FLOAT (Digits); DECLARE (x, y) FLOAT (Digits); ON SUBSCRIPTRANGE SNAP BEGIN; PUT SKIP DATA (IZ); END; IF ^initialized THEN CALL zigset(jsr); hz = shr3(); iz = IAND(hz, 127); IF (ABS(hz) < kn(iz)) THEN RETURN (hz * wn(iz)); ELSE DO FOREVER; IF (iz = 0) THEN DO; DO FOREVER; x = -LOG(uni())*0.2904764; y = -LOG(uni()); IF (y+y >= x*x) THEN LEAVE; END; IF (hz > 0) THEN RETURN (r+x); ELSE RETURN (-r-x); END; x = hz * wn(iz); IF (fn(iz) + uni()*(fn(iz-1)-fn(iz)) < EXP(-half*x*x)) THEN RETURN (x); hz = shr3(); iz = IAND(hz, 127); IF (ABS(hz) < kn(iz)) THEN RETURN (hz * wn(iz)); END; END rnor; rexp: PROCEDURE () RETURNS (FLOAT (Digits)); /* generate random exponentials. */ DECLARE x FLOAT (Digits); IF ^initialized THEN CALL zigset(jsr); jz = shr3(); iz = IAND(jz, 255); IF (ABS(jz) < ke(iz)) THEN RETURN ( ABS(jz) * we(iz)); DO FOREVER; IF (iz = 0) THEN RETURN (7.69711000000E0 - LOG(uni())); x = ABS(jz) * we(iz); IF (fe(iz) + uni()*(fe(iz-1) - fe(iz)) < EXP(-x)) THEN RETURN (x); jz = shr3(); iz = IAND(jz, 255); IF (ABS(jz) < ke(iz)) THEN RETURN ( ABS(jz) * we(iz)); END; END rexp; DECLARE (t1, t2) FLOAT (Digits); DECLARE (x, chisq, expctd) FLOAT (Digits); DECLARE (i, ten_million INITIAL (10000000), freq(0:999), j) FIXED BINARY (31); /* Set the random number seed. */ PUT ( ' Enter a non-zero integer: '); GET (i); CALL zigset(i); /* First time the generation of uniform, normal & exponential generators. */ t1 = SECS(); DO i = 1 TO ten_million; x = uni(); END; t2 = SECS (); PUT SKIP EDIT ( ' Time to generate ', ten_million, ' uniform random nos. = ', t2-t1, 'sec.') ( a, F(10), a, F(9,2), a); t1 = SECS; DO i = 1 TO ten_million; x = rnor(); END; t2 = SECS; PUT SKIP EDIT ( ' Time to generate ', ten_million, ' normal random nos. = ', t2-t1, 'sec.') ( a, F(10), a, F(9,2), a); t1 = SECS; DO i = 1 TO ten_million; x = rexp(); END; t2 = SECS; PUT SKIP EDIT ( ' Time to generate ', ten_million, ' exponential random nos. = ', t2-t1, 'sec.') (a, F(10), a, F(9,2), a); /* Now look at the distributions generated. */ PUT SKIP (2) LIST ( 'Std. devn. of chi-squared with 999 d.of.f. = 44.70'); PUT SKIP LIST ( 'Values of chi-squared below should be within about 90. of 999.'); freq = 0; DO i = 1 TO ten_million; j = 1000 * uni(); freq(j) = freq(j) + 1; END; chisq = 0.0; expctd = ten_million / 1000; DO j = 0 TO 999; chisq = chisq + (freq(j) - expctd)**2 / expctd; END; PUT SKIP EDIT ( ' Uniform distribution, chi-squared = ', chisq, ' with 999 deg. of freedom' ) (a, F(9,2), a); freq = 0; DO i = 1 TO ten_million; j = 1000 * alnorm(rnor(), FALSE); freq(j) = freq(j) + 1; END; chisq = 0.0; expctd = ten_million / 1000; DO j = 0 TO 999; chisq = chisq + (freq(j) - expctd)**2 / expctd; END; PUT SKIP EDIT ( ' Normal distribution, chi-squared = ', chisq, ' with 999 deg. of freedom' ) (a, F(9,2), a); freq = 0; DO i = 1 TO ten_million; j = 1000 * EXP(-rexp()); IF (j > 999 | j < 0) THEN PUT SKIP EDIT ( ' i, j = ', i, j) (a, 2 F(10)); ELSE freq(j) = freq(j) + 1; END; chisq = 0.0; expctd = ten_million / 1000; DO j = 0 TO 999; chisq = chisq + (freq(j) - expctd)**2 / expctd; END; PUT SKIP EDIT ( ' Exponential distribution, chi-squared = ', chisq, ' with 999 deg. of freedom' ) (a, F(9,2), a); STOP; alnorm: PROCEDURE (x, upper) RETURNS(FLOAT(Digits)); /* Algorithm AS66 Applied Statistics (1973) vol.22, no.3 */ /* Evaluates the tail area of the standardised normal curve */ /* from x to infinity if upper is .true. or */ /* from minus infinity to x if upper is .false. */ /* ELF90-compatible version by Alan Miller */ /* Latest revision - 29 November 1997 */ DECLARE x FLOAT (Digits); DECLARE upper BIT(1); DECLARE fn_val FLOAT(Digits); /* Local variables */ DECLARE (zero VALUE (0.0E0), one VALUE (1.0E0), half VALUE (0.5E0), con VALUE (1.28000000000E0)) FLOAT (Digits); DECLARE (z, y) FLOAT (Digits); DECLARE up BIT (1); /**** machine dependent constants. */ DECLARE (ltone VALUE (7.0E0), utzero VALUE (18.6600000000E0)) FLOAT (Digits); DECLARE (p VALUE (0.398942280444E0), q VALUE (0.39990348504E0), r VALUE (0.398942280385E0), a1 VALUE (5.75885480458E0), a2 VALUE (2.62433121679E0), a3 VALUE (5.92885724438E0), b1 VALUE (-29.8213557807E0), b2 VALUE (48.6959930692E0), c1 VALUE (-3.805200000000E-8), c2 VALUE (3.98064794000E-4), c3 VALUE (-0.151679116635E0), c4 VALUE (4.8385912808E0), c5 VALUE (0.742380924027E0), c6 VALUE (3.99019417011E0), d1 VALUE (1.00000615302E0), d2 VALUE (1.98615381364E0), d3 VALUE (5.29330324926E0), d4 VALUE (-15.1508972451E0), d5 VALUE (30.789933034E0) ) FLOAT (Digits); up = upper; z = x; IF z < zero THEN DO; up = ^up; z = -z; END; IF (z <= ltone) | (up & z <= utzero) THEN y = half*z*z; ELSE DO; fn_val = zero; GO TO L40; END; IF(z > con) THEN fn_val = r*EXP(-y)/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6)))))); ELSE fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3)))); L40: IF ^up THEN fn_val = one - fn_val; RETURN (fn_val); END alnorm; END test_ziggurat;