## Gamma Distributions

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.

Here's a screenshot showing a combination of the pdf and the random deviates, to see how well they agree: ```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```