Beta Distribution in Nimrod

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