yaole
116 lines
2.4 KiB
Go
116 lines
2.4 KiB
Go
package ope
|
|
|
|
import (
|
|
"math"
|
|
)
|
|
|
|
// uniformRand returns a random number in [0, 1]
|
|
// len(coins) >= 32
|
|
func uniformRand(coins []byte) float64 {
|
|
b := int64(0)
|
|
for i := 0; i < 32; i++ {
|
|
b <<= 1
|
|
b |= int64(coins[i])
|
|
}
|
|
return float64(b) / float64(uint64(1)<<32-1)
|
|
}
|
|
|
|
// hypergeometric10 returns a sample in hypergeometric distribution.
|
|
// Pr[X=x] = choose(y, x) * choose(N-y, M-x) / choose(N, M)
|
|
func hypergeometricSmall(M int64, N int64, y int64, coins []byte) int64 {
|
|
d1 := N - y
|
|
d2 := min(M, N-M)
|
|
|
|
Y := d2
|
|
K := y
|
|
for Y > 0.0 && K != 0 {
|
|
U := uniformRand(coins)
|
|
Y -= int64(U + float64(Y)/float64(d1+K))
|
|
K--
|
|
}
|
|
Z := d2 - Y
|
|
if N-M < M {
|
|
Z = y - Z
|
|
}
|
|
return Z
|
|
}
|
|
|
|
func loggam[T int | int64 | float32 | float64](x T) float64 {
|
|
return math.Log(math.Gamma(float64(x)))
|
|
}
|
|
|
|
func hypergeometric(M int64, N int64, y int64, coins []byte) int64 {
|
|
if y > 10 {
|
|
return hypergeometricSmall(M, N, y, coins)
|
|
}
|
|
return hypergeometricBig(M, N, y, coins)
|
|
}
|
|
|
|
func hypergeometricBig(M int64, N int64, y int64, coins []byte) int64 {
|
|
D1 := 1.7155277699214135
|
|
D2 := 0.8989161620588988
|
|
// # long mingoodbad, maxgoodbad, popsize, m, d9;
|
|
// # double d4, d5, d6, d7, d8, d10, d11;
|
|
// # long Z;
|
|
// # double T, W, X, Y;
|
|
good := M
|
|
bad := N - M
|
|
sample := y
|
|
|
|
mingoodbad := min(good, bad)
|
|
popsize := N
|
|
maxgoodbad := max(good, bad)
|
|
m := min(sample, popsize-sample)
|
|
d4 := float64(mingoodbad) / float64(popsize)
|
|
d5 := 1.0 - d4
|
|
d6 := float64(m)*d4 + 0.5
|
|
d7 := math.Sqrt(float64(popsize-m)*float64(sample)*d4*d5/float64(popsize-1) + 0.5)
|
|
d8 := D1*d7 + D2
|
|
d9 := (m + 1) * (mingoodbad + 1) / (popsize + 2)
|
|
d10 := loggam(d9+1) + loggam(mingoodbad-d9+1) + loggam(m-d9+1) + loggam(maxgoodbad-m+d9+1)
|
|
d11 := min(float64(min(m, mingoodbad)+1), math.Floor(d6+16*d7))
|
|
// # 16 for 16-decimal-digit precision in D1 and D2
|
|
|
|
var Z int64
|
|
for {
|
|
X := uniformRand(coins)
|
|
Y := uniformRand(coins)
|
|
W := d6 + d8*(Y-0.5)/X
|
|
|
|
// # fast rejection:
|
|
if W < 0.0 || W >= d11 {
|
|
continue
|
|
}
|
|
|
|
Z = int64(math.Floor(W))
|
|
T := d10 - (loggam(Z+1) + loggam(mingoodbad-Z+1) + loggam(m-Z+1) + loggam(maxgoodbad-m+Z+1))
|
|
|
|
// # fast acceptance:
|
|
if (X*(4.0-X) - 3.0) <= T {
|
|
break
|
|
}
|
|
|
|
// # fast rejection:
|
|
if X*(X-T) >= 1 {
|
|
continue
|
|
}
|
|
|
|
// # acceptance:
|
|
if 2.0*math.Log(X) <= T {
|
|
break
|
|
}
|
|
}
|
|
|
|
// # this is a correction to HRUA* by Ivan Frohne in rv.py
|
|
if good > bad {
|
|
Z = m - Z
|
|
}
|
|
|
|
// # another fix from rv.py to allow sample to exceed popsize/2
|
|
if m < sample {
|
|
Z = good - Z
|
|
}
|
|
|
|
return Z
|
|
}
|