From 318b8da414bdcaad9912a92a5fe92fcc496c8a7a Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.semperelh@protonmail.com>
Date: Sun, 9 Jun 2024 23:08:15 +0200
Subject: [PATCH] add beta sampler and implement beta type conforming to
 distribution

---
 f.go             | 17 +++++++++++++++--
 sample/sample.go |  6 +++++-
 2 files changed, 20 insertions(+), 3 deletions(-)

diff --git a/f.go b/f.go
index 415d83f..ef0ef12 100644
--- a/f.go
+++ b/f.go
@@ -21,21 +21,30 @@ type Distribution interface {
 	Samples() []float64
 }
 
+// Lognormal implementing Distribution
 type Lognormal struct {
 	low  float64
 	high float64
 }
 
-func (l Lognormal) Samples() []float64 {
-	sampler := func(r sample.Src) float64 { return sample.Sample_to(l.low, l.high, r) }
+func (ln Lognormal) Samples() []float64 {
+	sampler := func(r sample.Src) float64 { return sample.Sample_to(ln.low, ln.high, r) }
 	return sample.Sample_parallel(sampler, 1_000_000)
 }
 
+// Beta implementing Distribution
 type Beta struct {
 	a float64
 	b float64
 }
 
+func (beta Beta) Samples() []float64 {
+	sampler := func(r sample.Src) float64 { return sample.Sample_beta(beta.a, beta.b, r) }
+	return sample.Sample_parallel(sampler, 1_000_000)
+
+}
+
+// FilledSamples implementing Distribution
 type FilledSamples struct {
 	xs []float64
 }
@@ -132,6 +141,10 @@ func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
 
 }
 
+func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
+	return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
+}
+
 func joinDists(old_dist Dist, new_dist Dist, op string) (Dist, error) {
 	switch {
 	case old_dist.Type == "Lognormal" && new_dist.Type == "Lognormal" && op == "*":
diff --git a/sample/sample.go b/sample/sample.go
index 26350e3..42ac5a8 100644
--- a/sample/sample.go
+++ b/sample/sample.go
@@ -64,7 +64,6 @@ func Sample_gamma(alpha float64, r Src) float64 {
 		d = alpha - 1.0/3.0
 		c = 1.0 / math.Sqrt(9.0*d)
 
-	OuterLoop:
 		for {
 
 		InnerLoop:
@@ -96,7 +95,12 @@ func Sample_gamma(alpha float64, r Src) float64 {
 	} else {
 		return Sample_gamma(1.0+alpha, r) * math.Pow(Sample_unit_uniform(r), 1.0/alpha)
 	}
+}
 
+func Sample_beta(a float64, b float64, r Src) float64 {
+	gamma_a := Sample_gamma(a, r)
+	gamma_b := Sample_gamma(b, r)
+	return gamma_a / (gamma_a + gamma_b)
 }
 
 func Sample_mixture(fs []func64, weights []float64, r Src) float64 {