Artifact
f69872047e8725fa3900ba410bdd9da67e603008:
Attachment "chisq.tcl" to
ticket [1861370fff]
added by
erickb
2007-12-31 20:58:21.
# 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