use of cdfnorminv function: how to make use of this function in a fortran code, provide working ex

ll1240 · ‎12-05-2025

! file: mkl_norminv_vml.f90
! ifx mkl_norminv_vml.f90 -qmkl

program test_mkl_cdfnorminv
use, intrinsic :: iso_c_binding
implicit none
interface
subroutine vdCdfNormInv(n, a, y) bind(C, name="cdfnorminv")
!import :: c_int, c_double
integer :: n
real(8) :: a(1)
real(8) :: y(1)
end subroutine vdCdfNormInv
end interface

real(8) :: p(1), x(1)
p(1) = 0.90d0
call vdCdfNormInv(1, p, x) ! x = Φ^{-1}(p)
WRITE(*,'((a10),f18.8)') "out vdCdF", x(1)
!print '(a,f8.4)', '90th percentile of N(0,1): ', x(1)
end program test_mkl_cdfnorminv

*** and provide all the instruction to compile , link and run, please.

i've tried implementing the Acklam algo, feeded with rand() , it provides numerical results but not gaussian distributed.

FUNCTION inv_normal_ack(p, mu, sigma) result(x)

! use of the instruction result

IMPLICIT NONE

REAL(8), INTENT(in) :: p, mu, sigma

REAL(8) :: x, z

REAL(8) :: q, r

! Coefficients for Acklam approximation

REAL(8), PARAMETER :: a1 = -3.969683028665376d0

REAL(8), PARAMETER :: a2 = 2.209460984245205d0

REAL(8), PARAMETER :: a3 = -2.759285104469687d0

REAL(8), PARAMETER :: a4 = 1.383577518672690d0

REAL(8), PARAMETER :: a5 = -3.066479806614716d-1

REAL(8), PARAMETER :: a6 = 2.506628277459239d-2

REAL(8), PARAMETER :: b1 = -5.447609879822406d0

REAL(8), PARAMETER :: b2 = 1.615858368580409d0

REAL(8), PARAMETER :: b3 = -1.556989798598866d-1

REAL(8), PARAMETER :: b4 = 6.680131188771972d-3

REAL(8), PARAMETER :: c1 = -7.784894002430293d-3

REAL(8), PARAMETER :: c2 = -3.223964580411365d-1

REAL(8), PARAMETER :: c3 = -2.400758277161838d0

REAL(8), PARAMETER :: c4 = -2.549732539343734d0

REAL(8), PARAMETER :: c5 = 4.374664141464968d0

REAL(8), PARAMETER :: c6 = 2.938163982698783d0

REAL(8), PARAMETER :: d1 = 7.784695709041462d-3

REAL(8), PARAMETER :: d2 = 3.224671290700398d-1

REAL(8), PARAMETER :: d3 = 2.445134137142996d0

REAL(8), PARAMETER :: d4 = 3.754408661907416d0

! Define breakpoints

REAL(8), PARAMETER :: plow = 0.02425d0

REAL(8), PARAMETER :: phigh = 1.0d0 - plow

IF (p <= 0.0d0 .or. p >= 1.0d0) THEN

WRITE(*,'(a)') "Error: p must be in [0 1]"

stop

END IF

! Compute z for standard normal

IF (p < plow) THEN

q = sqrt(-2.0d0 * log(p))

z = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / &

((((d1*q+d2)*q+d3)*q+d4)*q+1.0d0)

ELSE IF (p > phigh) THEN

q = sqrt(-2.0d0 * log(1.0d0 - p))

z = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / &

((((d1*q+d2)*q+d3)*q+d4)*q+1.0d0)

ELSE

q = p - 0.5d0

r = q*q

z = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / &

(((((b1*r+b2)*r+b3)*r+b4)*r+1.0d0))

END IF

! Scale to given mean and std dev

x = mu + sigma * z

RETURN

END FUNCTION inv_normal_ack

JohnNichols · ‎12-05-2025

! file: mkl_norminv_vml.f90
! ifx mkl_norminv_vml.f90 -qmkl

program test_mkl_cdfnorminv
use, intrinsic :: iso_c_binding
implicit none
interface
subroutine vdCdfNormInv(n, a, y) bind(C, name="cdfnorminv")
!import :: c_int, c_double
integer :: n
real(8) :: a(1)
real(8) :: y(1)
end subroutine vdCdfNormInv
end interface

real(8) :: p(1), x(1)
p(1) = 0.90d0
call vdCdfNormInv(1, p, x) ! x = Φ^{-1}(p)
WRITE(*,'((a10),f18.8)') "out vdCdF", x(1)
!print '(a,f8.4)', '90th percentile of N(0,1): ', x(1)
end program test_mkl_cdfnorminv

  FUNCTION inv_normal_ack(p, mu, sigma) result(x)
            ! use of the instruction result 
        IMPLICIT NONE
        REAL(8), INTENT(in) :: p, mu, sigma
        REAL(8) :: x, z
        REAL(8) :: q, r
 
        ! Coefficients for Acklam approximation
        REAL(8), PARAMETER :: a1 = -3.969683028665376d0
        REAL(8), PARAMETER :: a2 =  2.209460984245205d0
        REAL(8), PARAMETER :: a3 = -2.759285104469687d0
        REAL(8), PARAMETER :: a4 =  1.383577518672690d0
        REAL(8), PARAMETER :: a5 = -3.066479806614716d-1
        REAL(8), PARAMETER :: a6 =  2.506628277459239d-2
        REAL(8), PARAMETER :: b1 = -5.447609879822406d0
        REAL(8), PARAMETER :: b2 =  1.615858368580409d0
        REAL(8), PARAMETER :: b3 = -1.556989798598866d-1
        REAL(8), PARAMETER :: b4 =  6.680131188771972d-3
        REAL(8), PARAMETER :: c1 = -7.784894002430293d-3
        REAL(8), PARAMETER :: c2 = -3.223964580411365d-1
        REAL(8), PARAMETER :: c3 = -2.400758277161838d0
        REAL(8), PARAMETER :: c4 = -2.549732539343734d0
        REAL(8), PARAMETER :: c5 =  4.374664141464968d0
        REAL(8), PARAMETER :: c6 =  2.938163982698783d0
        REAL(8), PARAMETER :: d1 =  7.784695709041462d-3
        REAL(8), PARAMETER :: d2 =  3.224671290700398d-1
        REAL(8), PARAMETER :: d3 =  2.445134137142996d0
        REAL(8), PARAMETER :: d4 =  3.754408661907416d0
 
        ! Define breakpoints
        REAL(8), PARAMETER :: plow = 0.02425d0
        REAL(8), PARAMETER :: phigh = 1.0d0 - plow
 
        IF (p <= 0.0d0 .or. p >= 1.0d0) THEN
            WRITE(*,'(a)') "Error: p must be in [0 1]"
            stop
        END IF
        ! Compute z for standard normal
        IF (p < plow) THEN
            q = sqrt(-2.0d0 * log(p))
            z = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / &
                ((((d1*q+d2)*q+d3)*q+d4)*q+1.0d0)
        ELSE IF (p > phigh) THEN
            q = sqrt(-2.0d0 * log(1.0d0 - p))
            z = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / &
                ((((d1*q+d2)*q+d3)*q+d4)*q+1.0d0)
        ELSE
            q = p - 0.5d0
            r = q*q
            z = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / &
                (((((b1*r+b2)*r+b3)*r+b4)*r+1.0d0))
        END IF
 
        ! Scale to given mean and std dev
        x = mu + sigma * z
      RETURN
      END FUNCTION inv_normal_ack

Question 1: Please explain why I should have to show you how to compile a Fortran program, the usual method is to read the manual.

Question 2: Why do you assume that RAND gives you Random numbers that are Gaussian, we know it is a lousy RND, has been since it was developed for Microsoft Fortran

Question 3: Is this a school assignment, a PhD thesis stuff or work related? All of those mean you should do it on your own to learn?

Also read the manual on posting properly. and *** and provide all the instruction to compile , link and run, please. == why should volunteers who spent years programming and learning do this for you?

These people help, but you have to make the first effort at getting it working.

mecej4O · ‎12-08-2025

You have already got the pieces that you need. Put the following lines into a Fortran source file, say, "tstcdf.f90":

program test_mkl_cdfnorminv
implicit none
integer :: n

real(8) :: p(1), x(1)
p(1) = 0.90d0
call vdCdfNormInv(1, p, x) ! x = F^{-1}(p)
WRITE(*,'((a10),f18.8)') "out vdCdF", x(1)
!print '(a,f8.4)', '90th percentile of N(0,1): ', x(1)
end program test_mkl_cdfnorminv

Compile and link using the command

ifx -qmkl tstcdf.f90

and run using the command

./a.out

and you should see the output

out vdCdF 1.28155157

You do not need Acklam's code -- it does the same thing that the MKL routine vdCdfNormInv does.