Fortran erfinv and erfinv functions

1

Dear all,

I was checking the methods of using uniform random number to generate gaussian random number, since I try to replace the uniform random number with the quasi-uniform random number (such as Sobol sequence) to generate quasi gaussian random number. It may generate some more ‘gaussian’ distributed gaussian random numbers.

Now, a simple approach to do so, can be using the erfinv and erfcinv functions as simply described in MARLAB page as below,

image
image3285×1477 143 KB

So they use the inverse of error function erfinc, and the inverse of the complimentary error function erfcinv.

I know Fortran has erf and erfc implemented already, but are there some Fortran implementations of erfinv and erfcinv functions?

Many thanks in advance!

1 Like
2

The functions sought are in the JPL Math77 library. See https://netlib.org/math/docpdf/ch02-13.pdf .

1 Like
3

A reasonable recent paper, with references to prior work, is

J. S. Dyer and S. A. Dyer, “Approximations to inverse error functions,” in IEEE Instrumentation & Measurement Magazine, vol. 11, no. 5, pp. 32-36, October 2008, doi: 10.1109/MIM.2008.4630740.

Abstract: We have presented four approximations to the inverse cumulative distribution function of a standard normal random variable. The first two were meant purely as a tool to motivate the other approximations, although either could be useful when extreme accuracy is not required and one doesn’t venture too far into the tails of the distribution. There are many other such “simple” approximations to F X -1 , suitable for use on a hand calculator.

1 Like
4

there is also stats_distribution_normal – Fortran-lang/stdlib

1 Like
5

I ported the erfinv function from Julia to Fortran, but lacked the motivation to push it across the finish line for stdlib. You can find it here:

The method in Julia uses the tables from Blair et al. (1976):

Blair, J. M., Edwards, C. A., & Johnson, J. H. (1976). Rational Chebyshev approximations for the inverse of the error function. Mathematics of Computation, 30(136), 827-830. AMS :: Mathematics of Computation

I’ve also seen an implementation for GPU before:

Giles, M. (2012). Approximating the erfinv function. In GPU Computing Gems Jade Edition (pp. 109-116). Morgan Kaufmann. Redirecting

1 Like
6

Alan Miller’s as241.f90 does this. Here it is with a short driver.

module normdev_mod
implicit none
private
public :: dp, nprob, normdev, normal
integer, parameter :: dp = selected_real_kind(15, 307)
contains
elemental function normal(z) result(prob)
real (dp), intent(in) :: z
real (dp)             :: prob
call nprob(z, p=prob)
end function normal
!
elemental subroutine nprob(z, p, q, pdf)
! For a given number (z) of standard deviations from the mean, the
! probabilities to the left (p) and right (q) of z are calculated,
! and the probability density (pdf).

! Programmer - Alan J. Miller
! Latest revision - 7 August 1997
! This version is compatible with the ELF90 subset of Fortran 90.

! REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE,
! ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969.

real(dp), intent(in)            :: z
real(dp), intent(out), optional :: p, q, pdf

! Local variables
real(dp), parameter :: a0 = 0.5d0, a1 = 0.398942280444d0, a2 = 0.399903438504d0, &
                    a3 = 5.75885480458d0, a4 = 29.8213557808d0,  &
                    a5 = 2.62433121679d0, a6 = 48.6959930692d0,  &
                    a7 = 5.92885724438d0, b0 = 0.398942280385d0, &
                    b1 = 3.8052d-8, b2 = 1.00000615302d0, b3 = 3.98064794d-4, &
                    b4 = 1.98615381364d0, b5 = 0.151679116635d0, &
                    b6 = 5.29330324926d0, b7 = 4.8385912808d0,   &
                    b8 = 15.1508972451d0, b9 = 0.742380924027d0, &
                    b10 = 30.789933034d0, b11 = 3.99019417011d0, &
                    zero = 0.d0, one = 1.d0
real(kind=dp) :: zabs, y, pp, qq, ppdf
zabs = abs(z)
if (zabs < 12.7d0) then
  y = a0*z*z
  ppdf = exp(-y)*b0
  if (present(pdf)) pdf = ppdf

!     Z BETWEEN -1.28 AND +1.28

  if (zabs < 1.28) then
    qq = a0 - zabs*(a1-a2*y/(y+a3-a4/(y+a5+a6/(y+a7))))
    if(z < zero) then
      pp = qq
      qq = one - pp
    else
      pp = one - qq
    end if
    if (present(p)) p = pp
    if (present(q)) q = qq
    return
  end if

!     ZABS BETWEEN 1.28 AND 12.7

  qq = ppdf/(zabs-b1+b2/(zabs+b3+b4/(zabs-b5+b6/(zabs+b7-b8/ &
            (zabs+b9+b10/(zabs+b11))))))
  IF(z < zero) THEN
    pp = qq
    qq = one - pp
  ELSE
    pp = one - qq
  END IF
  IF (PRESENT(p)) p = pp
  IF (PRESENT(q)) q = qq
  RETURN

!     Z FAR OUT IN TAIL

ELSE
  ppdf = zero
  IF(z < zero) THEN
    pp = zero
    qq = one
  ELSE
    pp = one
    qq = zero
  END IF
  IF (PRESENT(p)) p = pp
  IF (PRESENT(q)) q = qq
  IF (PRESENT(pdf)) pdf = ppdf
  RETURN
END IF
END SUBROUTINE nprob
!
elemental FUNCTION normdev(p) RESULT(fn_val)
REAL (dp), INTENT(IN) :: p
REAL (dp)             :: fn_val
! Local variable
INTEGER :: ifault
CALL ppnd16(p, fn_val, ifault)
if (ifault /= 0) fn_val = huge(fn_val)
end function normdev
!
elemental SUBROUTINE ppnd16 (p, normal_dev, ifault)

! ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3

! Produces the normal deviate Z corresponding to a given lower
! tail area of P; Z is accurate to about 1 part in 10**16.

! The hash sums below are the sums of the mantissas of the
! coefficients.   They are included for use in checking
! transcription.

! This ELF90-compatible version by Alan Miller - 20 August 1996
! N.B. The original algorithm is as a function; this is a subroutine

REAL (dp), INTENT(IN)  :: p
INTEGER, INTENT(OUT)   :: ifault
REAL (dp), INTENT(OUT) :: normal_dev

! Local variables

REAL (dp), parameter :: zero = 0.d0, one = 1.d0, half = 0.5d0,  &
                    split1 = 0.425d0, split2 = 5.d0, const1 = 0.180625d0, &
                    const2 = 1.6d0

! Coefficients for P close to 0.5

REAL (dp), parameter :: a0 = 3.3871328727963666080D0, &
                    a1 = 1.3314166789178437745D+2, &
                    a2 = 1.9715909503065514427D+3, &
                    a3 = 1.3731693765509461125D+4, &
                    a4 = 4.5921953931549871457D+4, &
                    a5 = 6.7265770927008700853D+4, &
                    a6 = 3.3430575583588128105D+4, &
                    a7 = 2.5090809287301226727D+3, &
                    b1 = 4.2313330701600911252D+1, &
                    b2 = 6.8718700749205790830D+2, &
                    b3 = 5.3941960214247511077D+3, &
                    b4 = 2.1213794301586595867D+4, &
                    b5 = 3.9307895800092710610D+4, &
                    b6 = 2.8729085735721942674D+4, &
                    b7 = 5.2264952788528545610D+3
! HASH SUM AB           55.8831928806149014439

! Coefficients for P not close to 0, 0.5 or 1.

REAL (dp), parameter :: c0 = 1.42343711074968357734D0, &
                    c1 = 4.63033784615654529590D0, &
                    c2 = 5.76949722146069140550D0, &
                    c3 = 3.64784832476320460504D0, &
                    c4 = 1.27045825245236838258D0, &
                    c5 = 2.41780725177450611770D-1, &
                    c6 = 2.27238449892691845833D-2, &
                    c7 = 7.74545014278341407640D-4, &
                    d1 = 2.05319162663775882187D0, &
                    d2 = 1.67638483018380384940D0, &
                    d3 = 6.89767334985100004550D-1, &
                    d4 = 1.48103976427480074590D-1, &
                    d5 = 1.51986665636164571966D-2, &
                    d6 = 5.47593808499534494600D-4, &
                    d7 = 1.05075007164441684324D-9
! HASH SUM CD           49.33206503301610289036

! Coefficients for P near 0 or 1.

REAL (dp), parameter :: e0 = 6.65790464350110377720D0, &
                    e1 = 5.46378491116411436990D0, &
                    e2 = 1.78482653991729133580D0, &
                    e3 = 2.96560571828504891230D-1, &
                    e4 = 2.65321895265761230930D-2, &
                    e5 = 1.24266094738807843860D-3, &
                    e6 = 2.71155556874348757815D-5, &
                    e7 = 2.01033439929228813265D-7, &
                    f1 = 5.99832206555887937690D-1, &
                    f2 = 1.36929880922735805310D-1, &
                    f3 = 1.48753612908506148525D-2, &
                    f4 = 7.86869131145613259100D-4, &
                    f5 = 1.84631831751005468180D-5, &
                    f6 = 1.42151175831644588870D-7, &
                    f7 = 2.04426310338993978564D-15
real(kind=dp) :: q, r
! HASH SUM EF           47.52583317549289671629

ifault = 0
q = p - half
if (abs(q) <= split1) then
  r = const1 - q * q
  normal_dev = q * (((((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2)*r + a1)*r + a0) / &
           (((((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2)*r + b1)*r + one)
  return
else
  if (q < zero) then
    r = p
  else
    r = one - p
  end if
  if (r <= zero) then
    ifault = 1
    normal_dev = zero
    return
  end if
  r = sqrt(-log(r))
  if (r <= split2) then
    r = r - const2
    normal_dev = (((((((c7*r + c6)*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + c0) / &
             (((((((d7*r + d6)*r + d5)*r + d4)*r + d3)*r + d2)*r + d1)*r + one)
  else
    r = r - split2
    normal_dev = (((((((e7*r + e6)*r + e5)*r + e4)*r + e3)*r + e2)*r + e1)*r + e0) / &
             (((((((f7*r + f6)*r + f5)*r + f4)*r + f3)*r + f2)*r + f1)*r + one)
  end if
  if (q < zero) normal_dev = - normal_dev
  return
end if
end subroutine ppnd16
end module normdev_mod
program xnormdev
! 12/27/2021 07:56 PM driver for normdev
use normdev_mod, only: dp, normdev
implicit none
real(kind=dp), parameter :: p(*) = [0.4_dp,0.5_dp,0.6_dp,0.7_dp,0.99_dp, &
   0.999_dp,0.9999_dp,0.99999_dp,0.999999_dp]
print "(/,a5,*(1x,f9.6))","p",p
print "(a5,*(1x,f9.6))","z",normdev(p)
end program xnormdev

output with gfortran or ifort:

    p  0.400000  0.500000  0.600000  0.700000  0.990000  0.999000  0.999900  0.999990  0.999999
    z -0.253347  0.000000  0.253347  0.524401  2.326348  3.090232  3.719016  4.264891  4.753424
1 Like
7

According to Generating low-discrepancy sequences from the normal distribution: Box–Muller or inverse transform? (2011) by Ökten and Göncü you can use Box-Muller on quasi-random uniform variates to generate quasi-random normal variates.

1 Like