PROGRAM t_random ! Test of module RANDOM ! Latest revision - 4 January 1999 ! Programmer - Alan.Miller @ vic.cmis.csiro.au USE random IMPLICIT NONE CHARACTER (LEN = 2) :: option INTEGER, ALLOCATABLE :: ix(:), seed(:) REAL, ALLOCATABLE :: x(:), xmean(:), cov(:), chol_f(:) INTEGER :: i, n, ndf, n_binom, k, ier, pos1, pos2, which REAL :: average, stdvn, shape, a, b, mu, p, pop_mean, pop_sd, & one = 1.0, zero = 0.0, sk, two = 2.0, pcntile(9), & middle, xmax, xmin LOGICAL :: first INTERFACE SUBROUTINE statistics(x, n, mean, median, xmax, xmin, stdev, percentile) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n) REAL, OPTIONAL, INTENT(OUT) :: mean, median, xmax, xmin, stdev, percentile(9) END SUBROUTINE statistics END INTERFACE ! Set the random number seed. CALL RANDOM_SEED(size=k) ALLOCATE (seed(k)) WRITE(*, '(1x, a, i4, a)') 'Enter ', k, ' integers as random number seeds: ' READ(*, *) seed CALL RANDOM_SEED(put=seed) DO WRITE(*, *)' 1. Normal (Gaussian) 2. Gamma' WRITE(*, *)' 3. Exponential 4. Beta' WRITE(*, *)' 5. t 6. Multivariate normal' WRITE(*, *)' 7. Generalized inverse Gaussian 8. Poisson' WRITE(*, *)' 9. Binomial 10. Negative binomial' WRITE(*, *)'11. von Mises 12. Cauchy' WRITE(*, *)'13. Weibull 14. Chi-squared' WRITE(*, *)'15. Random order of integers' WRITE(*, *) WRITE(*, *)' Q Quit' WRITE(*, *) WRITE(*, *)' Enter your choice: ' READ(*, '(a)') option option = ADJUSTL(option) IF (option .EQ. 'q ' .OR. option .EQ. 'Q ') STOP SELECT CASE (option) CASE ('1 ') WRITE(*, *) 'Distribution: Normal (Gaussian)' WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) DO i = 1, n x(i) = random_normal() END DO CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) WRITE(*, '(1x, a, f9.4, 2x, a, f9.4, 2x, a, f9.4)') & 'Mean = ', average, 'Median = ', middle, 'Sample st.dev. = ', stdvn WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile 900 FORMAT(" 10% 20% 30% 40% 50% 60% 70% 80% 90% "/ & 9f8.2) WRITE(*, *) DEALLOCATE (x) CASE ('2 ') WRITE(*, *) 'Distribution: Gamma' WRITE(*, *) 'Enter +ve value for shape parameter: ' READ(*, *) shape WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_gamma(shape, first) first = .false. END DO CALL statistics(x, n, average, middle, stdev=stdvn) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Median = ', middle, ' Std.devn. = ', stdvn pop_mean = shape pop_sd = SQRT(shape) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, *) CASE ('3 ') WRITE(*, *) 'Distribution: Exponential' WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) DO i = 1, n x(i) = random_exponential() END DO CALL statistics(x, n, average, stdev=stdvn, percentile=pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn pop_mean = one pop_sd = one WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, 900) pcntile WRITE(*, *) CASE ('4 ') WRITE(*, *) 'Distribution: Beta' WRITE(*, *) 'Enter values of a, b: ' READ(*, *) a, b WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_beta(a, b, first) first = .false. END DO CALL statistics(x, n, average, stdev=stdvn, percentile=pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn pop_mean = a/(a+b) pop_sd = SQRT(a*b/(a+b+1)) / (a+b) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, 900) pcntile WRITE(*, *) CASE ('5 ') WRITE(*, *) 'Distribution: t' WRITE(*, *) 'How many degrees of freedom?: ' READ(*, *) ndf WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) DO i = 1, n x(i) = random_t(ndf) END DO CALL statistics(x, n, average, stdev=stdvn, percentile=pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn IF (ndf .GT. 2) THEN pop_mean = zero pop_sd = SQRT(ndf / REAL(ndf-2)) WRITE(*, '(a, 2g14.6)') ' Population mean & st.devn. = ', pop_mean, & pop_sd END IF WRITE(*, 900) pcntile WRITE(*, *) CASE ('6 ') WRITE(*, *) 'Distribution: Multivariate normal' WRITE(*, *) 'Enter no. of variables: ' READ(*, *) k n = k*(k+1)/2 ALLOCATE ( x(k), xmean(k), cov(n), chol_f(n) ) WRITE(*, *) 'Enter means of the variables: ' READ(*, *) xmean WRITE(*, *) 'Enter the covariance matrix, row-by-row' pos2 = 0 DO i = 1, k pos1 = pos2 + 1 pos2 = pos2 + i WRITE(*, '(1x, a, i2, a, i2, a, i2, a)') 'Row ', i, ' Columns ', 1, & ' - ', i, ' : ' READ(*, *) cov(pos1:pos2) END DO WRITE(*, *) 'How many random deviates?: ' READ(*, *) n first = .true. DO i = 1, n CALL random_mvnorm(k, xmean, cov, chol_f, first, x, ier) IF (ier .NE. 0) THEN WRITE(*, *) '** Covariance matrix is not +ve definite **' EXIT END IF first = .false. WRITE(*, '(1x, i5, 10f7.3/ (6x, 10f7.3))') i, x END DO DEALLOCATE (x, xmean, cov, chol_f) WRITE(*, *) CASE ('7 ') WRITE(*, *) 'Distribution: Generalized inverse Gaussian' WRITE(*, *) 'Enter values of h, b: ' READ(*, *) a, b WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_inv_gauss(a, b, first) first = .false. END DO CALL statistics(x, n, average, stdev=stdvn, percentile=pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn WRITE(*, 900) pcntile WRITE(*, *) CASE ('8 ') WRITE(*, *) 'Distribution: Poisson' WRITE(*, *) 'Enter +ve value for mean: ' READ(*, *) mu WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_Poisson(mu, first) first = .false. END DO CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn pop_mean = mu pop_sd = SQRT(mu) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) CASE ('9 ') WRITE(*, *) 'Distribution: Binomial' DO WRITE(*, *) 'Enter 1 or 2: ' READ(*, *) which IF (which == 1 .OR. which == 2) EXIT END DO WRITE(*, *) 'Enter parameters N & p: ' READ(*, *) n_binom, p WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. IF (which == 1) THEN DO i = 1, n x(i) = random_binomial1(n_binom, p, first) first = .false. END DO ELSE DO i = 1, n x(i) = random_binomial2(n_binom, p, first) first = .false. END DO END IF CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn pop_mean = n_binom * p pop_sd = SQRT(pop_mean * (one-p)) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) CASE ('10') WRITE(*, *) 'Distribution: Negative binomial' WRITE(*, *) 'Enter parameters k & p: ' READ(*, *) sk, p WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_neg_binomial(sk, p) first = .false. END DO CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Std.devn. = ', stdvn pop_mean = sk * p / (one-p) pop_sd = SQRT(pop_mean / (one-p)) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) CASE ('11') WRITE(*, *) 'Distribution: von Mises' WRITE(*, *) 'Enter +ve value for scale parameter: ' READ(*, *) sk WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) first = .true. DO i = 1, n x(i) = random_von_Mises(sk, first) first = .false. END DO CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) DEALLOCATE (x) WRITE(*, '(1x, a, f9.4, 2x, a, f9.4, 2x, a, f9.4)') & 'Mean = ', average, 'Median = ', middle, 'Sample st.dev. = ', stdvn WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) CASE ('12') WRITE(*, *) 'Distribution: Cauchy' WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE ( x(n) ) DO i = 1, n x(i) = random_Cauchy() END DO CALL statistics(x, n, average, middle, xmax, xmin, percentile=pcntile) WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Mean = ', average, 'Median = ', middle WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) DEALLOCATE ( x ) CASE ('13') WRITE(*, *) 'Distribution: Weibull' WRITE(*, *) 'Enter +ve value for scale parameter: ' READ(*, *) sk WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE (x(n)) DO i = 1, n x(i) = random_Weibull(sk) END DO CALL statistics(x, n, average, middle, stdev=stdvn) DEALLOCATE (x) WRITE(*, *)'Mean = ', average, ' Median = ', middle, ' Std.devn. = ', stdvn pop_mean = EXP( lngamma( DBLE(one/sk + one) ) ) pop_sd = SQRT( lngamma( DBLE(two/sk + one) - pop_mean**2) ) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, *) CASE ('14') WRITE(*, *) 'Distribution: Chi-squared' WRITE(*, *) 'Enter number of degrees of freedom (integer): ' READ(*, *) ndf WRITE(*, *) 'How many random deviates?: ' READ(*, *) n ALLOCATE ( x(n) ) first = .true. DO i = 1, n x(i) = random_chisq(ndf, first) first = .false. END DO CALL statistics(x, n, average, middle, xmax, xmin, stdvn, pcntile) WRITE(*, *)'Mean = ', average, ' Median = ', middle, ' Std.devn. = ', stdvn pop_mean = ndf pop_sd = SQRT(two * ndf) WRITE(*, '(1x, a, 2g14.6)') 'Population mean & st.devn. = ', pop_mean, & pop_sd WRITE(*, '(1x, a, f9.4, 2x, a, f9.4)') & 'Smallest x = ', xmin, 'Largest x = ', xmax WRITE(*, 900) pcntile WRITE(*, *) DEALLOCATE ( x ) CASE ('15') WRITE(*, *) 'Random order of integers 1 .. N' WRITE(*, *) 'Enter N: ' READ(*, *) n ALLOCATE ( ix(n) ) CALL random_order(ix, n) WRITE(*, '(1x, 15I5)') ix WRITE(*, *) DEALLOCATE ( ix ) CASE DEFAULT WRITE(*, *) 'Option ', option, ' not allowed' END SELECT END DO STOP END PROGRAM t_random SUBROUTINE statistics(x, n, mean, median, xmax, xmin, stdev, percentile) ! Calculate the requested statistics for the sample x(1), ..., x(n). ! The percentiles estimated are 10% , 20% , ..., 90%. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n) REAL, OPTIONAL, INTENT(OUT) :: mean, median, xmax, xmin, stdev, percentile(9) INTERFACE SUBROUTINE qsort(a, n, t) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN OUT) :: a(:) INTEGER, INTENT(IN OUT) :: t(:) END SUBROUTINE qsort END INTERFACE ! Local variables INTEGER, ALLOCATABLE :: original_pos(:) REAL, ALLOCATABLE :: xx(:) INTEGER :: i, j REAL :: pos, remainder, half = 0.5, one = 1.0 ! n must be > 0 IF (n <= 0) THEN WRITE(*, *)'N must be > 0, value of n input to subroutine statistics = ', n RETURN END IF ! If requested, calculate the mean & standard deviation IF (PRESENT(mean) .OR. PRESENT(stdev)) THEN mean = SUM(x)/n IF(PRESENT(stdev)) THEN stdev = SUM( (x-mean)**2 ) IF (n > 1) stdev = SQRT(stdev/(n-1)) END IF END IF ! If median or percentiles have been requested, copy x to xx and sort. IF (PRESENT(median) .OR. PRESENT(percentile)) THEN ALLOCATE (xx(n), original_pos(n)) xx = x DO i = 1, n original_pos(i) = i END DO CALL qsort(xx, n, original_pos) IF (PRESENT(median)) THEN IF (n == 2*(n/2)) THEN i = n/2 median = half * (xx(i) + xx(i+1)) ELSE i = (n+1)/2 median = xx(i) END IF END IF IF (PRESENT(percentile)) THEN DO i = 1, 9 pos = i * REAL(n+1) / 10. j = pos remainder = pos - j percentile(i) = (one - remainder)*xx(j) + remainder*xx(j+1) END DO END IF END IF IF (PRESENT(xmin) .OR. PRESENT(xmax)) THEN xmin = MINVAL(x) xmax = MAXVAL(x) END IF RETURN END SUBROUTINE statistics SUBROUTINE qsort(a, n, t) ! NON-RECURSIVE STACK VERSION OF QUICKSORT FROM N. WIRTH'S PASCAL ! BOOK, 'ALGORITHMS + DATA STRUCTURES = PROGRAMS'. ! ALSO CHANGES THE ORDER OF THE ASSOCIATED ARRAY T. ! 30 stacks are enough for about 150,000 < n < 200,000 ! Corrected by Jaoa Martins, Oeiras, Portugal ! Latest revision - 4 January 1999 IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN OUT) :: a(:) INTEGER, INTENT(IN OUT) :: t(:) ! Local Variables INTEGER :: k, low, high, low0, high0, s, stackl(30), stackr(30), ww REAL :: w, x stackl(1) = 1 stackr(1) = n s = 1 low0 = 1 high0 = n 10 low = low0 high = high0 k = (low + high)/2 x = a(k) DO WHILE (low < high) DO WHILE (a(low) < x) low = low + 1 END DO DO WHILE (a(high) > x) high = high - 1 END DO IF (low <= high) THEN w = a(low) ww = t(low) a(low) = a(high) t(low) = t(high) a(high) = w t(high) = ww END IF END DO IF (low0 < high) THEN stackl(s) = high stackr(s) = high0 high0 = stackl(s) - 1 low0 = low0 - 1 IF (low0 <= 0) low0 = 1 s = s + 1 GO TO 10 END IF IF (low < high0) THEN low0 = low + 1 GO TO 10 END IF IF (s > 1) THEN s = s - 1 low0 = stackl(s) high0 = stackr(s) GO TO 10 END IF RETURN END SUBROUTINE qsort