From d57300e8ce1154a715b96712898a5564051b9b66 Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.semperelh@protonmail.com>
Date: Sun, 9 Jun 2024 23:01:49 +0200
Subject: [PATCH] add implementation for gamma distribution

---
 sample/sample.go | 52 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)

diff --git a/sample/sample.go b/sample/sample.go
index b82ca64..26350e3 100644
--- a/sample/sample.go
+++ b/sample/sample.go
@@ -47,6 +47,58 @@ func Sample_to(low float64, high float64, r Src) float64 {
 	return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r))
 }
 
+func Sample_gamma(alpha float64, r Src) float64 {
+
+	// a simple method for generating gamma variables, marsaglia and wan tsang, 2001
+	// https://dl.acm.org/doi/pdf/10.1145/358407.358414
+	// see also the references/ folder
+	// note that the wikipedia page for the gamma distribution includes a scaling parameter
+	// k or beta
+	// https://en.wikipedia.org/wiki/gamma_distribution
+	// such that gamma_k(alpha, k) = k * gamma(alpha)
+	// or gamma_beta(alpha, beta) = gamma(alpha) / beta
+	// so far i have not needed to use this, and thus the second parameter is by default 1.
+
+	if alpha >= 1 {
+		var d, c, x, v, u float64
+		d = alpha - 1.0/3.0
+		c = 1.0 / math.Sqrt(9.0*d)
+
+	OuterLoop:
+		for {
+
+		InnerLoop:
+			for {
+				x = Sample_unit_normal(r)
+				v = 1.0 + c*x
+				if v > 0.0 {
+					break InnerLoop
+				}
+			}
+
+			v = v * v * v
+			u = Sample_unit_uniform(r)
+
+			if u < 1.0-0.0331*(x*x*x*x) { // Condition 1
+				// the 0.0331 doesn't inspire much confidence
+				// however, this isn't the whole story
+				// by knowing that Condition 1 implies condition 2
+				// we realize that this is just a way of making the algorithm faster
+				// i.e., of not using the logarithms
+				return d * v
+			}
+			if math.Log(u) < 0.5*(x*x)+d*(1.0-v+math.Log(v)) { // Condition 2
+				return d * v
+			}
+
+		}
+
+	} else {
+		return Sample_gamma(1.0+alpha, r) * math.Pow(Sample_unit_uniform(r), 1.0/alpha)
+	}
+
+}
+
 func Sample_mixture(fs []func64, weights []float64, r Src) float64 {
 
 	// fmt.Println("weights initially: ", weights)