'=======================================================================
' program probit
'
' purpose:  estimates the probit function
' source:   Abramowitz M, Stegun IA, Handbook of mathematical functions,
'           Eq. 26.2.23 (p.933).  Washington: U. S. Government Printing
'           Office, 1972
'
' input:    user supplies proportion or probability, p, at prompt
' output:   z value such that the integral under the standard normal
'           curve from -inf to z equals p.
' 
' note:     |error| < 4.5E-4
'
' written:     02/13/93
' revised:     change(s):
' --------     ---------------------------------------------------------
'=======================================================================
      defdbl a-z

      cls
      print "(to quit, enter 0)"

  10  print ""
      input "Probability: ", p

      if (p <= 0 or p >= 1 ) then 100

      gosub 1000
      print using "z-score:     ##.###"; xp
      goto 10
 100  end

1000 ' probit function approximation with Abramowitz and Stegan, Eq. 26.2.23

       c0 =  2.515517
       c1 =   .802853
       c2 =   .010328

       d1 =  1.432788
       d2 =   .189269
       d3 =   .001308

       iflag  =   1
       if (p  >  .5) then p = 1. - p : iflag = 0

       t = (log(1./p^2))^.5

      xnum =       c0 + c1 * t + c2 * t^2
       xden =  1.000000 + d1 * t + d2 * t^2 + d3 * t^3
       xp   =  t - xnum/xden

       if (iflag) then xp = -xp

       return
       end