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