156 lines
3.3 KiB
Go
156 lines
3.3 KiB
Go
package main
|
|
|
|
import (
|
|
"fmt"
|
|
"git.nunosempere.com/NunoSempere/probppl/choose"
|
|
"math"
|
|
rand "math/rand/v2"
|
|
)
|
|
|
|
type src = *rand.Rand
|
|
type pplKnownDistrib = map[int64]float64
|
|
|
|
func generatePeopleKnownDistribution(r src) map[int64]float64 {
|
|
mapping := make(map[int64]float64)
|
|
sum := 0.0
|
|
|
|
// Consider zero case separately
|
|
/*p0 := r.Float64()
|
|
mapping[0.0] = p0
|
|
sum += p0
|
|
*/
|
|
|
|
// Consider successive exponents of 1.5
|
|
num := 16.0
|
|
base := 2.0
|
|
for i := 1; i < 8; i++ {
|
|
num = num * base
|
|
p := r.Float64()
|
|
mapping[int64(num)] = p
|
|
sum += p
|
|
}
|
|
|
|
for key, value := range mapping {
|
|
mapping[key] = value / sum
|
|
}
|
|
|
|
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 pp <= sum {
|
|
return int64(i)
|
|
}
|
|
|
|
}
|
|
fmt.Printf("%f, %f\n", sum, pp)
|
|
fmt.Println(d)
|
|
panic("Probabilities should sum up to 1")
|
|
}
|
|
|
|
func aboutEq(a int64, b int64) bool {
|
|
h := int64(3)
|
|
return ((-h) <= (a - b)) && ((a - b) <= h)
|
|
}
|
|
|
|
func draw148PplFromDistributionAndCheck(d pplKnownDistrib, r src, show bool) 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) {
|
|
if show {
|
|
// fmt.Println(count)
|
|
}
|
|
if aboutEq(count[0], 69) && aboutEq(count[1], 46) && aboutEq(count[2], 19) && aboutEq(count[3], 14) {
|
|
return 1
|
|
} else {
|
|
return 0
|
|
}
|
|
}
|
|
|
|
func getUnnormalizedBayesianUpdateForDistribution(d pplKnownDistrib, r src) int64 {
|
|
var sum int64 = 0
|
|
n := 100
|
|
for i := 0; i < n; i++ {
|
|
/* if i%1000 == 0 {
|
|
fmt.Println(i)
|
|
} */
|
|
draw_result := draw148PplFromDistributionAndCheck(d, r, i == 0)
|
|
// fmt.Println(draw_result)
|
|
sum += draw_result
|
|
}
|
|
return sum // float64(sum) / float64(n)
|
|
}
|
|
|
|
func main() {
|
|
|
|
var r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
|
|
|
|
sum := int64(0)
|
|
for i := 0; i < 1000; i++ {
|
|
|
|
people_known_distribution := generatePeopleKnownDistribution(r)
|
|
// fmt.Println(people_known_distribution)
|
|
result := getUnnormalizedBayesianUpdateForDistribution(people_known_distribution, r)
|
|
fmt.Println(i)
|
|
if result > 0 {
|
|
fmt.Println(people_known_distribution)
|
|
fmt.Println(result)
|
|
}
|
|
sum += result
|
|
// fmt.Println(result)
|
|
}
|
|
fmt.Println(sum)
|
|
|
|
}
|