# Implements the (central) chi-squared distribution
# as a special case of the gamma distribution

source simpletest.tcl
source gammadist.tcl

namespace eval chisq {
    
    # Values were produced from R ver. 2.4.1
    
    st::addproc cdf-chisquare {
        df  ;# Degrees of freedom
        x   ;# chisquared value; find prob that < this
    } where {
        {{2 3.5} ~ 0.826226056549555}
        {{5 2.2} ~ 0.179164030785504}
        {{5 100} ~ 1.0}
        {{3.9 4.2} ~ 0.634682741547709}
        {{1 2.0} ~ 0.842700792949715}
        {{3 -2.0} -> 0.0}
    }
    
    st::addproc pdf-chisquare {
        df  ;# Degrees of freedom
        x   ;# chisquared value; eval prob dens at this value
    } where {
        {{3 1.75} ~ 0.219999360547348}
        {{10 2.9} ~ 0.0216024880121444}
        {{4 17.45} ~ 0.000708787557977144}
        {{2.5 1.8} ~ 0.218446210041615}
        {{1 0.0} -> Inf}
        {{5 -1.5} -> 0.0}
    }
}

proc chisq::cdf-chisquare {df x} {
    set alpha [expr {0.5 * $df}]
    set beta 0.5
    
    if {$x < 0.0} {return 0.0}
    
    cdf-gamma $alpha $beta $x
}

proc chisq::pdf-chisquare {df x} {
    set alpha [expr {0.5 * $df}]
    set beta 0.5
    
    if {$x < 0.0} {return 0.0}
    
    pdf-gamma $alpha $beta $x
}

# Prodce "number" random deviates distributed according to a chisq dist
proc chisq::random-chisquare {df number} {
    set alpha [expr {0.5 * $df}]
    set beta 0.5
    
    random-gamma $alpha $beta $number    
}

console show
namespace eval chisq {
    st::print
}

### Test random deviates

##########################################################################################
##
## TESTING
##
## Can test pdf & cdf by running in a console. For random numbers, generate histograms:
##
##########################################################################################
canvas .c
pack .c -side top
frame .f
pack .f -side bottom
label .f.dfl -text "Degrees of freedom: "
entry .f.dfe -textvariable df
pack .f.dfl -side left
pack .f.dfe -side left
button .f.run -text "Run" -command runtest
pack .f.run -side left

proc runtest {} {
    set vals [chisq::random-chisquare $::df 5000]
    set remainder 5000
    for {set x 0.0} {$x < 20.0} {set x [expr {$x + 0.25}]} {
        lappend bins $x
        set distval [chisq::pdf-chisquare $::df $x]
        # Total trials (5000) * step (0.25) * value
        set distval [expr {int(5000 * 0.25 * $distval)}]
        lappend distcounts $distval
    }
    # Assume none are left
    lappend distcounts 0.0
    set bincounts [::math::statistics::histogram $bins $vals]
    .c delete hist
    .c delete dist
    math::statistics::plot-scale .c 0 20 0 [math::statistics::max $bincounts]
    math::statistics::plot-histogram .c $bincounts $bins hist
    math::statistics::plot-histogram .c $distcounts $bins dist
    .c itemconfigure dist -fill {} -outline green
}

set df 3
runtest
