I needed to implement the Beta distribution in Nimrod, which also involved implementing the Gamma distribution. I based (copied) the implementation on the Python random module.

I imported Nimrod random module using nimble, the Nimrod package manager. I’d never used nimble before, but I was able to follow the instruction here very easily.

import math
import random

proc gamma_variate( alpha, beta: float ): float =
    if( alpha <= 0.0 )or( beta <= 0.0 ):
        raise
    if alpha > 1.0:
        var u1, u2, v, x, z, r: float
        var cond1, cond2: bool
        let ainv = math.sqrt( 2.0 * alpha - 1.0 )
        let bbb = alpha - math.ln( 4.0 )
        let ccc = alpha + ainv
        let SG_MAGICCONST = 1.0 + math.ln( 4.5 )
        while true:
            u1 = random.random()
            if( u1 < 1e-7 )or( u1 > 0.9999999 ):
                continue
            u2 = 1.0 - random.random()
            v = math.ln( u1 / ( 1.0 - u1 ) ) / ainv
            x = alpha * math.exp( v )
            z = u1 * u1 * u2
            r = bbb + ccc*v - x
            cond1 = r + SG_MAGICCONST - 4.5*z >= 0.0
            cond2 = r >= math.ln( z )
            if cond1 or cond2:
                return x * beta
    if alpha == 1.0:
        var u = random.random()
        while u < 1e-7:
            u = random.random()
        return -math.ln(u) * beta
    if alpha < 1.0:
        var u, b, p, x, u1: float
        while true:
            u = random.random()
            b = ( math.E + alpha ) / math.E
            p = b * u
            if p <= 1.0:
                x = math.pow( p, ( 1.0 / alpha ) )
            else:
                x = -math.ln( ( b - p ) / alpha )
            u1 = random.random()
            if p > 1.0:
                if u1 <= math.pow( x, ( alpha - 1.0 ) ):
                    break
            elif u1 <= math.exp(-x):
                break
        return x * beta

proc beta_variate( alpha, beta: float ): float =
    let y = gamma_variate( alpha, 1.0 )
    if y == 0.0:
        return 0.0
    else:
        return y / ( y + gamma_variate( beta, 1.0 ) )