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