From 78d0c7da87b400c7cad1dc62a47b6193c330b546 Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Sun, 25 Feb 2024 17:31:22 -0300
Subject: [PATCH] get mvp running!

---
 probppl.go | 116 ++++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 92 insertions(+), 24 deletions(-)

diff --git a/probppl.go b/probppl.go
index 057b002..43e2a29 100644
--- a/probppl.go
+++ b/probppl.go
@@ -1,19 +1,15 @@
 package main
 
 import "fmt"
-import rand "math/rand/v2"
-
-/*
+import "git.nunosempere.com/NunoSempere/probppl/choose"
 import "math"
-import "sync"
 import rand "math/rand/v2"
 
-*/
-
 type src = *rand.Rand
+type pplKnownDistrib = map[int64]float64
 
-func generateRandomMapping(r src) map[int]float64 {
-	mapping := make(map[int]float64)
+func generatePeopleKnownDistribution(r src) map[int64]float64 {
+	mapping := make(map[int64]float64)
 	sum := 0.0
 
 	// Consider zero case separately
@@ -27,7 +23,7 @@ func generateRandomMapping(r src) map[int]float64 {
 	for i := 1; i < 22; i++ {
 		l = l * m
 		p := r.Float64()
-		mapping[int(l)] = p
+		mapping[int64(l)] = p
 		sum += p
 	}
 
@@ -35,27 +31,99 @@ func generateRandomMapping(r src) map[int]float64 {
 		mapping[key] = value / sum
 	}
 
-	/*
-		for i := 0; i <= 1000; i++ {
-			if current_probability < 1.0 {
-				num_possible_steps := int((1-current_probability)/0.001) + 1
-				fmt.Println("num possible steps: %d", num_possible_steps)
-				step := float64(r.IntN(num_possible_steps)) * 0.001
-				current_probability += step
-			}
-			mapping[i] = current_probability
-		}
-	*/
 	return mapping
 }
 
+func chooseWrapper(n int64, k int64) int64 {
+	if n < k {
+		return 0
+	} else {
+		return choose.Choose(n, k)
+	}
+}
+
+func getProbabilityOfKBirthdayMatchesGivenNPeopleKnown(n int64, k int64) float64 {
+	return float64(chooseWrapper(n, k)) * math.Pow(1.0/365.0, float64(k)) * math.Pow(1.0-(1.0/365.0), float64(n-k))
+}
+
+func getMatchesDrawGivenNPeopleKnown(n int64, r src) int64 {
+	p0 := getProbabilityOfKBirthdayMatchesGivenNPeopleKnown(n, 0)
+	p1 := getProbabilityOfKBirthdayMatchesGivenNPeopleKnown(n, 1)
+	p2 := getProbabilityOfKBirthdayMatchesGivenNPeopleKnown(n, 2)
+
+	p := r.Float64()
+	if p < p0 {
+		return 0
+	} else if p < (p0 + p1) {
+		return 1
+	} else if p < (p0 + p1 + p2) {
+		return 2
+	} else {
+		return 3 // stands for 'greater than 3'
+	}
+}
+
+/*
+  Draw 148 times
+  How many people do you know that were born in the same day of the year as you?
+	0: 46.6% | 69
+	1: 31.1% | 46
+	2: 12.8% | 19
+	≥3: 9.5% | 14
+*/
+
+func drawFromDistributionWithReplacement(d pplKnownDistrib, r src) int64 {
+	pp := r.Float64()
+	sum := 0.0
+	for i, p := range d {
+		sum += p
+		if p <= sum {
+			return int64(i)
+		}
+
+	}
+	fmt.Printf("%f, %f\n", sum, pp)
+	fmt.Println(d)
+	panic("Probabilities should sum up to 1")
+}
+
+func draw148PplFromDistributionAndCheck(d pplKnownDistrib, r src) int64 {
+
+	count := make(map[int64]int64)
+	count[0] = 0
+	count[1] = 0
+	count[2] = 0
+	count[3] = 0
+	for i := 0; i < 148; i++ {
+		person_i_ppl_known := drawFromDistributionWithReplacement(d, r)
+		person_i_num_birthday_matches := getMatchesDrawGivenNPeopleKnown(person_i_ppl_known, r)
+		count[person_i_num_birthday_matches]++
+	}
+	if (count[0] == 69) && (count[1] == 46) && (count[2] == 19) && (count[3] == 14) {
+		return 1
+	} else {
+		return 0
+	}
+}
+
+func getUnnormalizedBayesianUpdateForDistribution(d pplKnownDistrib, r src) float64 {
+	var sum int64 = 0
+	n := 1_000
+	for i := 0; i < n; i++ {
+		sum += draw148PplFromDistributionAndCheck(d, r)
+	}
+	return float64(sum) / float64(n)
+}
+
 func main() {
 	fmt.Println("Hello world")
 
 	var r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
-	randomMapping := generateRandomMapping(r)
-	for i, value := range randomMapping {
-		fmt.Printf("%d: %.4f\n", i, value)
-	}
+	people_known_distribution := generatePeopleKnownDistribution(r)
+	/* for i, p := range people_known_distribution {
+		fmt.Print64f("%d: %.4f\n", i, p)
+	} */
+	result := getUnnormalizedBayesianUpdateForDistribution(people_known_distribution, r)
+	fmt.Println(result)
 
 }