EKB This is a small package that implements some functions for gamma distributions. (It's also been submitted to tcllib.)
It implements cdf-gamma, pdf-gamma, and random-gamma (which generates gamma-distributed random deviates). It has a little test code at the bottom to generate histogram plots of the (pseudo)randomly-generated values.
It uses the code for the Incomplete gamma, read in as "incgamma.tcl".
EKB Already found an error in the random-number generation routine. Now fixed.
package require math
package require math::statistics
source "incgamma.tcl"
# Set of functions for the gamma distribution
# Implemented by Eric Kemp-Benedict, 2007
#
# This uses the following parameterization for the gamma:
# GammaDist(x) = beta * (beta * x)^(alpha-1) e^(-beta * x) / GammaFunc(alpha)
# Here, alpha is the shape parameter, and beta is the rate parameter
# Alternatively, a "scale parameter" theta = 1/beta is sometimes used
proc pdf-gamma {alpha beta x} {
if {$beta < 0} {
error "Rate parameter 'beta' must be be positive"
return
}
set prod [expr {1.0 * $x * $beta}]
set Galpha [expr {exp([::math::ln_Gamma $alpha])}]
expr {(1.0 * $beta/$Galpha) * pow($prod, ($alpha - 1.0)) * exp(-$prod)}
}
proc cdf-gamma {alpha beta x} {
incompleteGamma [expr {$beta * $x}] $alpha
}
# Generate a list of gamma-distributed random deviates
# Use Cheng's envelope rejection method, as documented
# in Dagpunar, J.S. 2007. "Simulation and Monte Carlo:
# With Applications in Finance and MCMC"
proc random-gamma {alpha beta number} {
if {$alpha <= 1} {
set lambda $alpha
} else {
set lambda [expr {sqrt(2.0 * $alpha - 1.0)}]
}
set retval {}
for {set i 0} {$i < $number} {incr i} {
while {1} {
# Two rands: one for deviate, one for acceptance/rejection
set r1 [expr {rand()}]
set r2 [expr {rand()}]
# Calculate deviate from enveloping proposal distribution (a Lorenz distribution)
set lnxovera [expr {(1.0/$lambda) * (log(1.0 - $r1) - log($r1))}]
if {![catch {expr {$alpha * exp($lnxovera)}} x]} {
# Apply acceptance criterion
if {log(4.0*$r1*$r1*$r2) < ($alpha - $lambda) * $lnxovera + $alpha - $x} {
break
}
}
}
lappend retval [expr {1.0 * $x/$beta}]
}
return $retval
}
##########################################################################################
##
## 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.alphal -text "alpha"
entry .f.alphae -textvariable alpha
pack .f.alphal -side left
pack .f.alphae -side left
label .f.betal -text "beta"
entry .f.betae -textvariable beta
pack .f.betal -side left
pack .f.betae -side left
button .f.run -text "Run" -command runtest
pack .f.run -side left
proc runtest {} {
set vals [random-gamma $::alpha $::beta 5000]
set remainder 5000
for {set x 0.0} {$x < 20.0} {set x [expr {$x + 0.25}]} {
lappend bins $x
set distvallow [pdf-gamma $::alpha $::beta $x]
set distvalhigh [pdf-gamma $::alpha $::beta [expr {$x + 0.25}]]
# Total trials (5000) * step (0.25) * average value
set distval [expr {int(5000 * 0.25 * 0.5 * ($distvallow + $distvalhigh))}]
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 alpha 1
set beta 0.5
runtest