Artifact
889654e0db42fa83e2e378cc873eca12d1b8f9c3:
Attachment "norminv.tcl" to
ticket [1838731fff]
added by
erickb
2007-11-26 21:17:25.
##
## Code by Eric Kemp-Benedict 2007
## Arjen and Tcllib team: Do with it as you wish
##
## quicknorminv is from an algorithm in Abramowitz and Stegun,
## otherwise just Newton's method
##
package require math::statistics
proc norminv {p} {
set maxiter 100
set epsilon 1.0e-9
set deltax [expr {100 * $epsilon}]
# Initial value for x
set x [quicknorminv $p]
set niter 0
while {abs([::math::statistics::pnorm $x] - $p) > $epsilon} {
set pstar [::math::statistics::pnorm $x]
set pl [::math::statistics::pnorm [expr {$x - $deltax}]]
set ph [::math::statistics::pnorm [expr {$x + $deltax}]]
set x [expr {$x + 2.0 * $deltax * ($p - $pstar)/($ph - $pl)}]
incr niter
if {$niter == $maxiter} {
error "Inverse normal distribution did not converge after $niter iterations"
return {}
}
}
return $x
}
proc quicknorminv {p} {
set epsilon 1.0e-9
if {$p < $epsilon || $p > 1-$epsilon} {
error "Value out of bounds for inverse normal: $p"
return 0
}
set sign -1
if {$p > 0.5} {
set p [expr {1 - $p}]
set sign 1
}
set t [expr {sqrt(-2.0 * log($p))}]
set a0 2.30753
set a1 0.27061
set b1 0.99229
set b2 0.04481
set num [expr {$a0 + $a1 * $t}]
set denom [expr {1 + $b1 * $t + $b2 * $t * $t}]
return [expr {$sign * ($t - $num/$denom)}]
}
#################################
##
## TEST
##
#################################
console show
puts "Orig x\tProb\tInv Prob"
set t [time {
for {set x -5.0} {$x < 5.0} {set x [expr {$x + 0.1}]} {
lappend xs $x
set p [::math::statistics::pnorm $x]
lappend ps $p
#lappend xprimes $x ;# Use this to time just pnorm, for baseline
lappend xprimes [norminv $p]
#lappend xprimes [quicknorminv $p]
#lappend xprimes [::math::statistics::Inverse-cdf-normal 0.0 1.0 $p]
}
}]
foreach x $xs p $ps xprime $xprimes {
puts "$x\t$p\t$xprime"
}
puts $t
## -------
## Results
## -------
# 1) No inverse calc (baseline): 3332 us/iteration
# 2) quicknorminv: 3817 us/iteration (+485 us/iteration = 0.0005 s)
# 3) norminvv: 26221 us/iteration (+22889 us/iteration = 0.02 s)
# 4) Inverse-cdf-normal: 228324 us/iteration (+224992 us/iteration = 0.2 s)
