##
## 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)
