# pnorm_ekb --
#     Calculate the cumulative distribution function (cdf)
#     for the standard normal distribution like in the statistical
#     software 'R' (mean=0 and sd=1)
#
# Arguments:
#     x               Value for which the cdf should be calculated
#
# Result:
#     Value of the statistic D
#
## The "_ekb" version includes edits by Eric Kemp-Benedict, who
## had to evaluate values at large x. This now switches between
## the power series (there already) and a continued fraction expansion.
##
## Continued fraction is from Abramowitz & Stegun, eqn 26.2.14
##
## Also some minor changes for clarity
proc pnorm_ekb {x {tol 1.0e-20}} {
    set overflow 1.0e37
    #
    # cumulative distribution function (cdf)
    # for the standard normal distribution like in the statistical software 'R'
    # (mean=0 and sd=1)
    #
    # x -> value for which the cdf should be calculated
    #
    
    # Calculate Z(x), the pdf at x. The constant is ln(sqrt(2 * pi))
    set Z [expr {exp(-0.5 * $x*$x - 0.91893853320467274178)}]
    
    # Somewhat arbitrary cutoff: value here is 1 within 1e-5
    if {abs($x) < 4.0} {
        # Power series
        set sum [expr {double($x)}]
        set oldSum 0.0
        set i 1
        set denom 1.0
        while {abs($sum - $oldSum) > $tol} {
                set oldSum $sum
                incr i 2
                set denom [expr {$denom*$i}]
                set sum [expr {$oldSum + pow($x,$i)/$denom}]
        }
        return [expr {0.5 + $sum *$Z}]
    } else {
        # Continued fraction
        set t [expr {abs($x)}]
        set a 1
        set first true
        set A0 1.0
        set B0 0.0
        set A1 0.0
        set B1 1.0
        set exp [expr {1.0 * $A1/$B1}]
        while {1} {
            set A2 [expr {$t * $A1 + $a * $A0}]
            set B2 [expr {$t * $B1 + $a * $B0}]
            if {$first} {
                set first false
            } else {
                incr a
            }
            set newexp [expr {1.0 * $A2/$B2}]
            if {abs($newexp - $exp) < $tol} {
                break
            }
            set exp $newexp
            set A0 $A1
            set B0 $B1
            set A1 $A2
            set B1 $B2
            # Too big? Rescale
            if {abs($A2) >= $overflow} {
                set A0 [expr {$A0 / $overflow}]
                set A1 [expr {$A1 / $overflow}]
                set B0 [expr {$B0 / $overflow}]
                set B1 [expr {$B1 / $overflow}]
            }
        }
        if {$x > 0} {
            return [expr {1.0 - $Z * $exp}]
        } else {
            return [expr {$Z * $exp}]
        }
    }
}

####################################################################################
##
## TEST
##
## Run for a large range of values and print result
##
###################################################################################
console show
for {set x -10} {$x <= 10} {set x [expr {$x + 0.2}]} {puts [pnorm_ekb $x]}

## The output from the test compared to Excel's NORMSDIST() function, with
##    selected values from R:
 # x    pnorm_ekb       Excel   Selected R values
 # -10.0 7.6198530E-24	0.0000000E+00  From R: 7.619853e-24
 # -9.8	5.6292823E-23	0.0000000E+00
 # -9.6	3.9972212E-22	0.0000000E+00
 # -9.4	2.7281536E-21	0.0000000E+00
 # -9.2	1.7897488E-20	0.0000000E+00
 # -9.0	1.1285884E-19	0.0000000E+00
 # -8.8	6.8408077E-19	0.0000000E+00
 # -8.6	3.9858050E-18	0.0000000E+00
 # -8.4	2.2323932E-17	0.0000000E+00
 # -8.2	1.2019352E-16	1.1102230E-16
 # -8.0	6.2209606E-16	6.6613381E-16
 # -7.8	3.0953588E-15	3.1086245E-15
 # -7.6	1.4806537E-14	1.4988011E-14
 # -7.4	6.8092249E-14	6.8611783E-14
 # -7.2	3.0106280E-13	3.0320191E-13
 # -7.0	1.2798125E-12	1.2880808E-12  From R: 1.279813e-12
 # -6.8	5.2309575E-12	5.2615690E-12
 # -6.6	2.0557889E-11	2.0665247E-11
 # -6.4	7.7688476E-11	7.8049012E-11
 # -6.2	2.8231580E-10	2.8347136E-10
 # -6.0	9.8658765E-10	9.9012187E-10
 # -5.8	3.3157460E-09	3.3260517E-09
 # -5.6	1.0717590E-08	1.0746217E-08
 # -5.4	3.3320448E-08	3.3396123E-08
 # -5.2	9.9644263E-08	9.9834400E-08
 # -5.0	2.8665157E-07	2.8710500E-07   From R: 2.866516e-07
 # -4.8	7.9332815E-07	7.9435267E-07
 # -4.6	2.1124547E-06	2.1146434E-06
 # -4.4	5.4125439E-06	5.4169531E-06
 # -4.2	1.3345749E-05	1.3354097E-05
 # -4.0	3.1671242E-05	3.1686035E-05
 # -3.8	7.2348044E-05	7.2372434E-05
 # -3.6	1.5910859E-04	1.5914571E-04
 # -3.4	3.3692927E-04	3.3698082E-04
 # -3.2	6.8713794E-04	6.8720208E-04
 # -3.0	1.3498980E-03	1.3499672E-03
 # -2.8	2.5551303E-03	2.5551906E-03
 # -2.6	4.6611880E-03	4.6612218E-03
 # -2.4	8.1975359E-03	8.1975289E-03
 # -2.2	1.3903448E-02	1.3903399E-02
 # -2.0	2.2750132E-02	2.2750062E-02
 # -1.8	3.5930319E-02	3.5930266E-02
 # -1.6	5.4799292E-02	5.4799289E-02
 # -1.4	8.0756659E-02	8.0756711E-02
 # -1.2	1.1506967E-01	1.1506973E-01
 # -1.0	1.5865525E-01	1.5865526E-01
 # -0.8	2.1185540E-01	2.1185533E-01
 # -0.6	2.7425312E-01	2.7425306E-01
 # -0.4	3.4457826E-01	3.4457830E-01
 # -0.2	4.2074029E-01	4.2074031E-01
 # 0.0	5.0000000E-01	5.0000000E-01
 # 0.2	5.7925971E-01	5.7925969E-01
 # 0.4	6.5542174E-01	6.5542170E-01
 # 0.6	7.2574688E-01	7.2574694E-01
 # 0.8	7.8814460E-01	7.8814467E-01
 # 1.0	8.4134475E-01	8.4134474E-01
 # 1.2	8.8493033E-01	8.8493027E-01
 # 1.4	9.1924334E-01	9.1924329E-01
 # 1.6	9.4520071E-01	9.4520071E-01
 # 1.8	9.6406968E-01	9.6406973E-01
 # 2.0	9.7724987E-01	9.7724994E-01
 # 2.2	9.8609655E-01	9.8609660E-01
 # 2.4	9.9180246E-01	9.9180247E-01
 # 2.6	9.9533881E-01	9.9533878E-01
 # 2.8	9.9744487E-01	9.9744481E-01
 # 3.0	9.9865010E-01	9.9865003E-01
 # 3.2	9.9931286E-01	9.9931280E-01
 # 3.4	9.9966307E-01	9.9966302E-01
 # 3.6	9.9984089E-01	9.9984085E-01
 # 3.8	9.9992765E-01	9.9992763E-01
 # 4.0	9.9996833E-01	9.9996831E-01
 # 4.2	9.9998665E-01	9.9998665E-01
 # 4.4	9.9999459E-01	9.9999458E-01
 # 4.6	9.9999789E-01	9.9999789E-01
 # 4.8	9.9999921E-01	9.9999921E-01
 # 5.0	9.9999971E-01	9.9999971E-01
 # 5.2	9.9999990E-01	9.9999990E-01
 # 5.4	9.9999997E-01	9.9999997E-01
 # 5.6	9.9999999E-01	9.9999999E-01
 # 5.8	1.0000000E+00	1.0000000E+00
 # 6.0	1.0000000E+00	1.0000000E+00
 # 6.2	1.0000000E+00	1.0000000E+00
 # 6.4	1.0000000E+00	1.0000000E+00
 # 6.6	1.0000000E+00	1.0000000E+00
 # 6.8	1.0000000E+00	1.0000000E+00
 # 7.0	1.0000000E+00	1.0000000E+00
 # 7.2	1.0000000E+00	1.0000000E+00
 # 7.4	1.0000000E+00	1.0000000E+00
 # 7.6	1.0000000E+00	1.0000000E+00
 # 7.8	1.0000000E+00	1.0000000E+00
 # 8.0	1.0000000E+00	1.0000000E+00
 # 8.2	1.0000000E+00	1.0000000E+00
 # 8.4	1.0000000E+00	1.0000000E+00
 # 8.6	1.0000000E+00	1.0000000E+00
 # 8.8	1.0000000E+00	1.0000000E+00
 # 9.0	1.0000000E+00	1.0000000E+00
 # 9.2	1.0000000E+00	1.0000000E+00
 # 9.4	1.0000000E+00	1.0000000E+00
 # 9.6	1.0000000E+00	1.0000000E+00
 # 9.8	1.0000000E+00	1.0000000E+00
 # 10.0	1.0000000E+00	1.0000000E+00
 ##
