hfnran Subroutine

public subroutine hfnran(N, Iseed, X)

NAME

hfnran(3f) - [M_datapac:RANDOM] generate half-normal random numbers

SYNOPSIS

   SUBROUTINE HFNRAN(N,Iseed,X)

    INTEGER,intent(in)         :: N
    INTEGER,intent(inout)      :: Iseed
    REAL(kind=wp),intent(out)  :: X(:)

DESCRIPTION

HFNRAN(3f) generates a random sample of size n from the halfnormal
distribution.

The prototype halfnormal distribution used herein has mean = sqrt(2/pi)
= 0.79788456 and standard deviation = 1. this distribution is defined
for all non-negative X and has the probability density function

    f(X) = (2/sqrt(2*pi)) * exp(-X*X/2).

The prototype halfnormal distribution used herein is the distribution
of the variate X = abs(Z) where the variate Z is normally distributed
with mean = 0 and standard deviation = 1.

INPUT ARGUMENTS

N      = The desired integer number of random numbers to be generated.

ISEED An integer iseed value. Should be set to a non-negative value to start a new sequence of values. Will be set to -1 on return to indicate the next call should continue the current random sequence walk.

OUTPUT ARGUMENTS

X      = A  vector (of dimension at least N) into which the generated
         random sample from the halfnormal distribution will be placed.

EXAMPLES

Sample program:

program demo_hfnran
use M_datapac, only : hfnran
implicit none
! call hfnran(x,y)
end program demo_hfnran

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCES

  • TOCHER, THE ART OF SIMULATION, 1963, pages 14-15.
  • HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, 1964, page 36.
  • JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE DISTRIBUTIONS–1, 1970, pages 53, 59, 81, 83.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N
integer, intent(inout) :: Iseed
real(kind=wp), intent(out) :: X(:)

Source Code

subroutine hfnran(N,Iseed,X)
integer,intent(in)         :: N
integer,intent(inout)      :: Iseed
real(kind=wp),intent(out)  :: X(:)
real(kind=wp)              :: arg1 , arg2 , sqrt1 , u1 , u2 , y(2) , z1 , z2
integer                    :: i , ip1
!
!  CHECK THE INPUT ARGUMENTS FOR ERRORS
!
   if ( N<1 ) then
      write (G_io,99001)
      99001 format (' ***** FATAL ERROR--The first input argument to HFNRAN(3f) is non-positive *****')
      WRITE (G_IO,99002) N
      99002 FORMAT (' ***** The value of the argument is ',I0,' *****')
      RETURN
   else
!
!  GENERATE N UNIFORM (0,1) RANDOM NUMBERS;
!  THEN GENERATE 2 ADDITIONAL UNIFORM (0,1) RANDOM NUMBERS
!  (TO BE USED BELOW IN FORMING THE N-TH NORMAL RANDOM NUMBER WHEN THE DESIRED SAMPLE SIZE N HAPPENS TO BE ODD).
!
      call uniran(N,Iseed,X)
      call uniran(2,Iseed,y)
!
!  GENERATE N NORMAL RANDOM NUMBERS USING THE BOX-MULLER METHOD.
!
      do i = 1 , N , 2
         ip1 = i + 1
         u1 = X(i)
         if ( i==N ) then
            u2 = y(2)
         else
            u2 = X(ip1)
         endif
         arg1 = -2.0_wp*LOG(u1)
         arg2 = 2.0_wp*G_pi*u2
         sqrt1 = SQRT(arg1)
         z1 = sqrt1*COS(arg2)
         z2 = sqrt1*SIN(arg2)
         X(i) = z1
         IF ( i/=N ) X(ip1) = z2
      enddo
!
!  GENERATE N HALFNORMAL RANDOM NUMBERS USING THE DEFINITION THAT A HALFNORMAL VARIATE
!  EQUALS THE ABSOLUTE VALUE OF A NORMAL VARIATE.
!
      do i = 1 , N
         if ( X(i)<0.0_wp ) X(i) = -X(i)
      enddo
   endif
!
end subroutine hfnran