Hello,

I'm submitting code for gamma-distributed random variables. It uses the file "incgamma.tcl" for the incomplete gamma function that I recently submitted.

The provided functions are:
  cdf-gamma
  pdf-gamma
  random-gamma

I've tested some values for the cdf & pdf against the corresponding functions pgamma & dgamma in R. For the random deviates, I looked at the histogram output (also provided in the package) and eyeballed the correspondence with plots of the gamma distribution.

Eric

Logged In: YES 
user_id=400048
Originator: NO

Incorporated this code

Logged In: YES 
user_id=816411
Originator: YES

Note that the zip file contains all of the fixes mentioned in previous posts. Also, it includes the incomplete gamma function. Finally, it puts everything in a "gamma" namespace (which I expect you'll move over to the math::statistics namespace).

File Deleted - 257779:

File Added - 260324: gamma.zip

Logged In: YES 
user_id=816411
Originator: YES

This is not a change to the basic code. Instead, it's just rearranging the testing code a bit. I'm attaching a zip file with all of the needed files to run the tests. The code itself is (almost) stand-alone: just a few calls to "simple test" (also included).

File Added: gamma.zip

Logged In: YES 
user_id=816411
Originator: YES

Also, adding a better way to test the random deviates. Now for the testing routine, add the pdf to the histogram of values. Here's the revised testing code:

##########################################################################################
##
## 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

Logged In: YES 
user_id=816411
Originator: YES

There was an error in random-gamma. Here is the corrected version:

# 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
}

File Added - 257779: gammadist.tcl