!**************************************************** !* Program to demonstrate inverse normal subroutine * !* ------------------------------------------------ * !* Reference: BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* !* * !* F90 version by J.-P. Moreau, Paris. * !* (www.jpmoreau.fr) * !* ------------------------------------------------ * !* SAMPLE RUN: * !* * !* P(Z>X) X * !* ---------------- * !* 0.50 0.0000 * !* 0.48 0.0500 * !* 0.46 0.1002 * !* 0.44 0.1507 * !* 0.42 0.2015 * !* 0.40 0.2529 * !* 0.38 0.3050 * !* 0.36 0.3580 * !* 0.34 0.4120 * !* 0.32 0.4673 * !* 0.30 0.5240 * !* 0.28 0.5825 * !* 0.26 0.6430 * !* 0.24 0.7060 * !* 0.22 0.7719 * !* 0.20 0.8414 * !* 0.18 0.9152 * !* 0.16 0.9944 * !* 0.14 1.0804 * !* 0.12 1.1751 * !* 0.10 1.2817 * !* 0.08 1.4053 * !* 0.06 1.5551 * !* 0.04 1.7511 * !* 0.02 2.0542 * !* -0.00 INF. * !* * !**************************************************** PROGRAM DEMO_INVNORM integer i real*8 x,y; print *,' ' print *,' P(Z>X) X ' print *,' ----------------' i = 1; y = 0.5d0 do while (y >= -0.01) call Inverse_Normal(x,y) if (x < 0.000001) x = 0.d0 if (x < 1e10) then write(*,50) y,x else write(*,60) y end if if (i.eq.20) then Pause 'Press any key to continue' i = 1 end if i=i+1 y = y - 0.02d0 end do print *, ' ' stop 50 format(' ',F5.2,' ',F6.4) 60 format(' ',F5.2,' INF.'/) end !*********************************************** !* Inverse normal distribution subroutine * !* ------------------------------------------- * !* This program calculates an approximation to * !* the integral of the normal distribution * !* function from x to infinity (the tail). * !* A rational polynomial is used. The input is * !* in y, with the result returned in x. The * !* accuracy is better then 0.0005 in the range * !* 0 < y < 0.5. * !* ------------------------------------------- * !* Reference: Abramowitz and Stegun. * !*********************************************** Subroutine Inverse_Normal(x,y) real*8 c0,c1,c2,d1,d2,d3,x,y,z ! Define coefficients c0 = 2.515517d0 c1 = 0.802853d0 c2 = 0.010328d0 d1 = 1.432788d0 d2 = 0.189269d0 d3 = 0.001308d0 if (y <= 0) then x = 1e13 return end if z = dsqrt(-dlog(y * y)) x = 1.d0 + d1 * z + d2 * z * z + d3 * z * z * z x = (c0 + c1 * z + c2 * z * z) / x x = z - x return end ! End of file invnorm.f90