Attachment "pnorm_ekb.tcl" to
ticket [1847125fff]
added by
erickb
2007-12-09 11:33:33.
# 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
##
