Tk Library Source Code

Artifact 3751eae7786cebd3b99aee91a1a3fffef57fe717:

Attachment "invchisqr.tcl" to ticket [1861370fff] added by erickb 2008-01-01 05:13:43.

---

source chisq.tcl

namespace eval chisq {
    variable maxiter 100
    
    # Values generated using R ver. 2.4.1
    st::addproc invchisq {
        df  ;# Degrees of freedom
        p   ;# probability
    } where {
        {{2 0.5} ~ 1.38629436111989}
        {{3 0.0} ~ 0.0}
        {{15 0.7} ~ 17.3216944984992}
        {{3 0.1} ~ 0.584374374155183}
        {{14 0.01} ~ 4.66042506265777}
        {{70 0.10} ~ 55.3289395719096}
        {{2 1.0} -> Inf}
        {{40 1.0} -> Inf}
    }
}

proc chisq::quicknorminv {p} {
    if {$p < 0.0 || $p > 1.0} {
        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)}]

}
 
proc chisq::invchisq {df p {tol 1.0e-10}} {
    variable maxiter
    
    if {$p < 0.0 || $p > 1.0} {
        error "Value out of bounds for inverse chi-square: $p"
        return 0.0
    }
    
    if {$p > 1.0 - $tol} {
        return [expr Inf]
    }
    
    set deltax [expr {100 * $tol}]

    # Initial value for x -- use mean = dof, unless large #dof
    set x $df
    if {$df > 30} {
        # Use approximate value from Abramowitz & Stegun for df > 30
        set p_norm [expr {1.0 - $p}]
        set xp [quicknorminv $p_norm]
        set y [expr {2.0/(9.0 * double($df))}]
        set x [expr {$df * (1.0 - $y + $xp * sqrt($y))}]
    }
    
    set niter 0
    while {1} {
        set pstar [cdf-chisquare $df $x]
        set dpdx [pdf-chisquare $df $x]
        if {abs($dpdx) < $tol} {
            set x [expr {$x + $deltax}]
            set dpdx [pdf-chisquare $df $x]
        }
        
        # Sometimes the slope can be quite large -- damp down size of steps if necessary
        set newx [expr {$x + ($p - $pstar)/$dpdx}]
        set newx [expr {$newx > $x + 0.5 * $df ? $x + 0.5 * $df : $newx}]
        set newx [expr {$newx < $x - 0.5 * $df ? $x - 0.5 * $df : $newx}]
        
        if {$newx < 0.0} {
            set newx 0.0
        }
        
        if {abs($newx - $x) < $tol} {
            set x $newx
            break
        } else {
            set x $newx
        }
        
        incr niter
        if {$niter == $maxiter} {
            error "Inverse chi-squared distribution did not converge after $niter iterations"
            return {}
        }
    }

    return $x

}

console show
st::tolerance 1.0e-9
namespace eval chisq {
    st::print
}