-
Notifications
You must be signed in to change notification settings - Fork 1.3k
/
block_proba.go
102 lines (91 loc) · 2.03 KB
/
block_proba.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
package messagepool
import (
"math"
"sync"
)
var noWinnersProbCache []float64
var noWinnersProbOnce sync.Once
func noWinnersProb() []float64 {
noWinnersProbOnce.Do(func() {
poissPdf := func(x float64) float64 {
const Mu = 5
lg, _ := math.Lgamma(x + 1)
result := math.Exp((math.Log(Mu) * x) - lg - Mu)
return result
}
out := make([]float64, 0, MaxBlocks)
for i := 0; i < MaxBlocks; i++ {
out = append(out, poissPdf(float64(i)))
}
noWinnersProbCache = out
})
return noWinnersProbCache
}
var noWinnersProbAssumingCache []float64
var noWinnersProbAssumingOnce sync.Once
func noWinnersProbAssumingMoreThanOne() []float64 {
noWinnersProbAssumingOnce.Do(func() {
cond := math.Log(-1 + math.Exp(5))
poissPdf := func(x float64) float64 {
const Mu = 5
lg, _ := math.Lgamma(x + 1)
result := math.Exp((math.Log(Mu) * x) - lg - cond)
return result
}
out := make([]float64, 0, MaxBlocks)
for i := 0; i < MaxBlocks; i++ {
out = append(out, poissPdf(float64(i+1)))
}
noWinnersProbAssumingCache = out
})
return noWinnersProbAssumingCache
}
func binomialCoefficient(n, k float64) float64 {
if k > n {
return math.NaN()
}
r := 1.0
for d := 1.0; d <= k; d++ {
r *= n
r /= d
n--
}
return r
}
func (mp *MessagePool) blockProbabilities(tq float64) []float64 {
noWinners := noWinnersProbAssumingMoreThanOne()
p := 1 - tq
binoPdf := func(x, trials float64) float64 {
// based on https://github.com/atgjack/prob
if x > trials {
return 0
}
if p == 0 {
if x == 0 {
return 1.0
}
return 0.0
}
if p == 1 {
if x == trials {
return 1.0
}
return 0.0
}
coef := binomialCoefficient(trials, x)
pow := math.Pow(p, x) * math.Pow(1-p, trials-x)
if math.IsInf(coef, 0) {
return 0
}
return coef * pow
}
out := make([]float64, 0, MaxBlocks)
for place := 0; place < MaxBlocks; place++ {
var pPlace float64
for otherWinners, pCase := range noWinners {
pPlace += pCase * binoPdf(float64(place), float64(otherWinners))
}
out = append(out, pPlace)
}
return out
}