NORPPF – M

public subroutine NORPPF(P, Ppf)

NAME

norppf(3f) - [M_datapac:PERCENT_POINT] compute the normal percent point function

SYNOPSIS

   SUBROUTINE NORPPF(P,Ppf)

DESCRIPTION

norppf(3f) computes the percent point function value for the normal
(gaussian) distribution with mean = 0 and standard deviation = 1.

this distribution is defined for all x and has the probability
density function

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

note that the percent point function of a distribution is identically
the same as the inverse cumulative distribution function of the
distribution.

OPTIONS

 X   description of parameter
 Y   description of parameter

EXAMPLES

Sample program:

program demo_norppf
use M_datapac, only : norppf
implicit none
! call norppf(x,y)
end program demo_norppf

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

ODEH AND EVANS, THE PERCENTAGE POINTS OF THE NORMAL DISTRIBUTION, ALGORITHM 70, APPLIED STATISTICS, 1974, pages 96-97.
EVANS, ALGORITHMS FOR MINIMAL DEGREE POLYNOMIAL AND RATIONAL APPROXIMATION, M. SC. THESIS, 1972, UNIVERSITY OF VICTORIA, B. C., CANADA.
HASTINGS, APPROXIMATIONS FOR DIGITAL COMPUTERS, 1955, pages 113, 191, 192.
NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS SERIES 55, 1964, page 933, FORMULA 26.2.23.
FILLIBEN, SIMPLE AND ROBUST LINEAR ESTIMATION OF THE LOCATION PARAMETER OF A SYMMETRIC DISTRIBUTION (UNPUBLISHED PH.D. DISSERTATION, PRINCETON UNIVERSITY), 1969, pages 21-44, 229-231.
FILLIBEN, ‘THE PERCENT POINT FUNCTION’, (UNPUBLISHED MANUSCRIPT), 1970, pages 28-31.
JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE DISTRIBUTIONS–1, 1970, pages 40-111.
THE KELLEY STATISTICAL TABLES, 1948.
OWEN, HANDBOOK OF STATISTICAL TABLES, 1962, pages 3-16.
PEARSON AND HARTLEY, BIOMETRIKA TABLES FOR STATISTICIANS, VOLUME 1, 1954, pages 104-113.

NAME

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

SYNOPSIS

   SUBROUTINE NORRAN(N,Iseed,X)

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

DESCRIPTION

NORRAN(3f) generates a random sample of size N from the normal
(Gaussian) distribution with mean = 0 and standard deviation = 1.

Internally, it uses the Box-Muller algorithm.

This distribution is defined for all X and has the probability
density function

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

INPUT ARGUMENTS

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

ISEED An integer seed 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 of size N from the normal distribution with mean =
     0 and standard deviation = 1 will be placed.

EXAMPLES

Sample program:

program demo_norran
use M_datapac, only : norran, label, plotxt, sort, norplt, plott
implicit none
integer,parameter :: N=300
real              :: x(N), y(N)
real              :: mu, sigma
integer           :: Iseed
integer           :: i
   Iseed=1234
   sigma=1.00000
   mu=0.0
   call label('norran')
   call norran(N,Iseed,x)
   x = sigma*x
   x = x + mu
   call plotxt(x,n)
   call sort(x,n,y) ! sort and replot to better discern distribution
   call plott([(real(i),i=1,n)],y,n)
end program demo_norran

Results:

 THE FOLLOWING IS A PLOT OF X(I) (VERTICALLY) VERSUS I (HORIZONTALLY
                   I-----------I-----------I-----------I-----------I
  0.3016713E+01 -                                               X
  0.2787551E+01 I
  0.2558388E+01 I
  0.2329226E+01 I     X
  0.2100063E+01 I
  0.1870901E+01 I     X   X XX      X XX      XX X
  0.1641738E+01 -         X    X         X
  0.1412575E+01 I    X        X X X  XX       X  X X         X
  0.1183413E+01 I                     X X        XX    X XXX   XX
  0.9542503E+00 I    X   XX          X             X    X  XX X X
  0.7250879E+00 I   X  XX X      X  X        XXX      XX     X X X
  0.4959254E+00 I     XX X  XXX   XXXXX   X  XX    X    X XX  XX  X
  0.2667627E+00 -    X XX  XXX X   XXX X X XX   X XXXX X  X     XX
  0.3760028E-01 I   X X    X   XX XXX  X   XXX X  X XXXX XX XX X  XX
 -0.1915622E+00 I  XX  X   X  X   X   X X X  X XXXX XX  XX X  X   X
 -0.4207249E+00 I  XX  XX   XX XXXX X   XX XX   X XXXX X X XXX XXX
 -0.6498873E+00 I        X XXX  XX  XX    XXXXXX    X XX    X     XX
 -0.8790498E+00 I   XX  X X   X  X X    XXX      X   X  XX       XX
 -0.1108212E+01 -        X     XXX     XXX  X         X        X
 -0.1337375E+01 I  X X    X                 X       X X   X XX X X
 -0.1566537E+01 I    X X          X       X               XX
 -0.1795700E+01 I  X    X   X   XX        X                        X
 -0.2024862E+01 I        X  X      X
 -0.2254025E+01 I                                     X         XX
 -0.2483188E+01 -            X
                   I-----------I-----------I-----------I-----------I
            0.1000E+01  0.7575E+02  0.1505E+03  0.2252E+03  0.3000E+03

 The following is a plot of Y(I) (vertically) versus X(I) (horizontally)
                   I-----------I-----------I-----------I-----------I
  0.3000000E+03 -                                        XX  X     X
  0.2875417E+03 I                                    XXXXX
  0.2750833E+03 I                                  XXX
  0.2626250E+03 I                                XXX
  0.2501667E+03 I                              XXX
  0.2377083E+03 I                             XX
  0.2252500E+03 -                           XXX
  0.2127917E+03 I                           X
  0.2003333E+03 I                          XX
  0.1878750E+03 I                         XX
  0.1754167E+03 I                        XX
  0.1629583E+03 I                        X
  0.1505000E+03 -                       XX
  0.1380417E+03 I                      XX
  0.1255833E+03 I                     XX
  0.1131250E+03 I                    XX
  0.1006667E+03 I                    X
  0.8820834E+02 I                   XX
  0.7575000E+02 -                  X
  0.6329167E+02 I                 XX
  0.5083334E+02 I               XX
  0.3837500E+02 I              XX
  0.2591669E+02 I           XXX
  0.1345834E+02 I       XXXXX
  0.1000000E+01 -  X X XX
                   I-----------I-----------I-----------I-----------I
           -0.2483E+01 -0.1108E+01  0.2668E+00  0.1642E+01  0.3017E+01

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

Box and Muller, ‘A Note on the Generation of Random Normal Deviates’, Journal of the Association for Computing Machinery, 1958, pages 610-611.
Tocher, The Art of Simulation, 1963, pages 33-34.
Hammersley and Handscomb, Monte Carlo Methods, 1964, page 39.
Johnson and Kotz, Continuous Univariate Distributions–1, 1970, pages 40-111.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp)				::	P
real(kind=wp)				::	Ppf

Source Code

SUBROUTINE NORPPF(P,Ppf)
REAL(kind=wp) :: aden , anum , P , p0 , p1 , p2 , p3 , p4 , Ppf , q0 , q1 , q2 , q3 , q4 , r , t
!
!     INPUT ARGUMENTS--P      = THE  VALUE
!                                (BETWEEN 0.0 AND 1.0)
!                                AT WHICH THE PERCENT POINT
!                                FUNCTION IS TO BE EVALUATED.
!     OUTPUT ARGUMENTS--PPF    = THE  PERCENT
!                                POINT FUNCTION VALUE.
!     OUTPUT--THE  PERCENT POINT
!             FUNCTION VALUE PPF.
!     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
!     RESTRICTIONS--P SHOULD BE BETWEEN 0.0 AND 1.0, EXCLUSIVELY.

!     COMMENTS--THE CODING AS PRESENTED BELOW
!               IS ESSENTIALLY IDENTICAL TO THAT
!               PRESENTED BY ODEH AND EVANS
!               AS ALGORITHM 70 OF APPLIED STATISTICS.
!               THE PRESENT AUTHOR HAS MODIFIED THE
!               ORIGINAL ODEH AND EVANS CODE WITH ONLY
!               MINOR STYLISTIC CHANGES.
!             --AS POINTED OUT BY ODEH AND EVANS
!               IN APPLIED STATISTICS,
!               THEIR ALGORITHM REPRESENTES A
!               SUBSTANTIAL IMPROVEMENT OVER THE
!               PREVIOUSLY EMPLOYED
!               HASTINGS APPROXIMATION FOR THE
!               NORMAL PERCENT POINT FUNCTION--
!               THE ACCURACY OF APPROXIMATION
!               BEING IMPROVED FROM 4.5*(10**-4)
!               TO 1.5*(10**-8).
!
!---------------------------------------------------------------------
!
      DATA p0 , p1 , p2 , p3 , p4/ - .322232431088_wp , -1.0_wp ,             &
     &     -.342242088547_wp , -.204231210245E-1_wp , -.453642210148E-4_wp/
      DATA q0 , q1 , q2 , q3 , q4/.993484626060E-1_wp , .588581570495_wp ,    &
     &     .531103462366_wp , .103537752850_wp , .38560700634E-2_wp/
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( P<=0.0_wp .OR. P>=1.0_wp ) THEN
         WRITE (G_IO,99001)
99001    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE NORPPF SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
         WRITE (G_IO,99002) P
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,       &
     &           ' *****')
         RETURN
!
!-----START POINT-----------------------------------------------------
!
      ELSEIF ( P/=0.5_wp ) THEN
!
         r = P
         IF ( P>0.5_wp ) r = 1.0_wp - r
         t = SQRT(-2.0_wp*LOG(r))
         anum = ((((t*p4+p3)*t+p2)*t+p1)*t+p0)
         aden = ((((t*q4+q3)*t+q2)*t+q1)*t+q0)
         Ppf = t + (anum/aden)
         IF ( P<0.5_wp ) Ppf = -Ppf
         GOTO 99999
      ENDIF
      Ppf = 0.0_wp
      RETURN
!
99999 END SUBROUTINE NORPPF

NORPPF Subroutine

Contents

Source Code

public subroutine NORPPF(P, Ppf)

NAME

SYNOPSIS

DESCRIPTION

OPTIONS

EXAMPLES

AUTHOR

MAINTAINER

LICENSE

REFERENCES

NAME

SYNOPSIS

DESCRIPTION

INPUT ARGUMENTS

OUTPUT ARGUMENTS

EXAMPLES

AUTHOR

MAINTAINER

LICENSE

REFERENCES

Arguments

Source Code