!********************************************************************** !* NORMAL AND INVERSE NORMAL PROBABILITY FUNCTIONS * !* ------------------------------------------------------------------ * !* SAMPLE RUN: * !* * !* Value of variable: 0.2 * !* * !* Probability = 0.391042693975456 * !* * !* Verify: * !* X = 0.200000286102295 * !* * !* F90 Version By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !********************************************************************** Program Test_Normal_Law real*8 P, X, Y, PHI, NORMAL print *,' ' write(*,10,advance='no'); read *, X P = PHI(X) print *,' ' print *,' Probability = ', P print *,' ' print *,' Verify:' Y = NORMAL(P) print *,' X = ', Y print *,' ' stop 10 format(' Value of variable: ') END real*8 Function PHI (u) real*8 u, PI ! Standard Normal Probability Function PI = 4.d0*datan(1.d0) PHI = (1.d0/dsqrt(2.d0 * PI))*dexp(-0.5d0*u*u) return end real*8 Function NORMAL (P) real*8 P, PHI ! Inverse Standard Normal Function real*8 chimax,chisqval,epsilon,minchisq,maxchisq epsilon = 1.d-6 chimax = 1.d6 minchisq = 0.d0 maxchisq = chimax; if (P <= 0.d0) then NORMAL = maxchisq return else if (P >= 1.0) then NORMAL = 0.d0 return end if end if chisqval = 0.5d0 do while (dabs(maxchisq - minchisq) > epsilon) if (PHI(chisqval) < P) then maxchisq = chisqval else minchisq = chisqval end if chisqval = (maxchisq + minchisq) * 0.5d0 end do NORMAL = chisqval return end ! end of file normal.f90