This is an implementation of the probability density function and a random number generator for Student's t-distribution. (See the Wikipedia article [L1 ].)
The function cdf-students-t is already implemented in tcllib, as is a procedure for calculating the t-statistic by calculating the inverse CDF.
This code implements pdf-students-t and random-students-t. Included is the code for the functions as well as some test code. The test code uses Simple Test.
package require math
namespace import math::ln_Gamma
source simpletest.tcl
namespace eval tdist {
variable pi [expr {atan2(0, -1)}]
st::addproc pdf-students-t {
degrees ;# Degrees of freedom (allow non-integer)
value ;# svalue where dist is to be evaluated
} where {
{{1 0.1} ~ 0.315158303152268}
{{0.5 0.1} ~ 0.265700672177405}
{{4 3.2} ~ 0.0156821741652879}
{{3 2.0} ~ 0.0675096606638929}
{{3 7.5} ~ 0.000942291548015668}
}
st::addproc random-students-t {
degrees ;# Degrees of freedom: may be non-integer, but >= 1
number ;# Number of random deviates to produce
} where {} ;# No tests: compare visually as histogram
}
proc tdist::pdf-students-t {degrees value} {
variable pi
set nplus1over2 [expr {0.5 * ($degrees + 1)}]
set f1 [expr {exp([ln_Gamma $nplus1over2] - \
[ln_Gamma [expr {$nplus1over2 - 0.5}]])}]
set f2 [expr {1.0/sqrt($degrees * $pi)}]
expr {$f1 * $f2 * pow(1.0 + $value * $value/double($degrees), -$nplus1over2)}
}
# Use method from Appendix 4.3 in Dagpunar, J.S.,
# "Simulation and Monte Carlo: With Applications in
# Finance and MCMC"
proc tdist::random-students-t {degrees number} {
variable pi
if {$degrees < 1} {
error "random-students-t: Degrees of freedom must be >= 1"
}
set dd [expr {double($degrees)}]
set k [expr {2.0/($dd - 1.0)}]
for {set i 0} {$i < $number} {incr i} {
set r1 [expr {rand()}]
if {$degrees > 1} {
set r2 [expr {rand()}]
set c [expr {cos(2.0 * $pi * $r2)}]
lappend retval [expr {sqrt($dd/ \
(1.0/(1.0 - pow($r1, $k)) \
- $c * $c)) * $c}]
} else {
lappend retval [expr {tan(0.5 * $pi * ($r1 + $r1 - 1))}]
}
}
set retval
} package require math::statistics
source tdist.tcl
st::tolerance 1.0e-9
console show
namespace eval tdist {
st::print
}
##########################################################################################
##
## TESTING
##
## For random numbers, generate histograms:
##
##########################################################################################
canvas .c
pack .c -side top
frame .f
pack .f -side bottom
label .f.dofl -text "dof"
entry .f.dofe -textvariable dof
pack .f.dofl -side left
pack .f.dofe -side left
button .f.run -text "Run" -command runtest
pack .f.run -side left
proc runtest {} {
set step 0.5
set vals [tdist::random-students-t $::dof 5000]
set remainder 5000
for {set x -10.0} {$x < 10.0} {set x [expr {$x + $step}]} {
lappend bins $x
set distval [tdist::pdf-students-t $::dof [expr {$x - 0.5 * $step}]]
# Total trials (5000) * step * value
set distval [expr {int(5000 * $step * $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 -7 7 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 dof 3
runtest