[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".

    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} {
        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 x [expr {$alpha * pow((1.0/$r1 - 1.0), 1.0/$beta)}]
                # Apply acceptance criterion
                if {log(4.0*$r1*$r1*$r2) < ($alpha - $beta) * log($x/(1.0*$alpha)) + $alpha - $x} {
                    break
                }
            }
            lappend retval $x
        }
        
        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]
        for {set x 0.0} {$x < 20.0} {set x [expr {$x + 0.25}]} {
            lappend bins $x
        }
        set bincounts [::math::statistics::histogram $bins $vals]
        .c delete hist
        math::statistics::plot-scale .c 0 20 0 [math::statistics::max $bincounts]
        math::statistics::plot-histogram .c $bincounts $bins hist
    }

    set alpha 1
    set beta 0.5
    runtest


----
!!!!!!
[Category Statistics] | [Category Mathematics]
!!!!!!