From 874488843a8ab78659ada598127fd2aba5559130 Mon Sep 17 00:00:00 2001 From: Agnish Ghosh Date: Thu, 25 Jul 2024 17:58:29 +0530 Subject: [PATCH] add: hypergeom cdf --- beacon_chain/spec/helpers.nim | 23 +++++ tests/test_helpers.nim | 187 ++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+) diff --git a/beacon_chain/spec/helpers.nim b/beacon_chain/spec/helpers.nim index b4521250a..b76c5d120 100644 --- a/beacon_chain/spec/helpers.nim +++ b/beacon_chain/spec/helpers.nim @@ -527,3 +527,26 @@ proc blockToBlockHeader*(blck: ForkyBeaconBlock): ExecutionBlockHeader = proc compute_execution_block_hash*(blck: ForkyBeaconBlock): Eth2Digest = rlpHash blockToBlockHeader(blck) + +from std/math import exp, ln +from std/sequtils import foldl + +func ln_binomial(n, k: int): float64 = + if k > n: + low(float64) + else: + template ln_factorial(n: int): float64 = + (2 .. n).foldl(a + ln(b.float64), 0.0) + ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k) + +func hypergeom_cdf*(k: int, population: int, successes: int, draws: int): + float64 = + if k < draws + successes - population: + 0.0 + elif k >= min(successes, draws): + 1.0 + else: + let ln_denom = ln_binomial(population, draws) + (0 .. k).foldl(a + exp( + ln_binomial(successes, b) + + ln_binomial(population - successes, draws - b) - ln_denom), 0.0) \ No newline at end of file diff --git a/tests/test_helpers.nim b/tests/test_helpers.nim index 1893e7a5b..8366bf35d 100644 --- a/tests/test_helpers.nim +++ b/tests/test_helpers.nim @@ -67,3 +67,190 @@ suite "Spec helpers": process(fieldVar, i shl childDepth) i += 1 process(state, state.numLeaves) + + test "hypergeom_cdf": + # Generated with SciPy's hypergeom.cdf() function + const tests = [ + ( 0, 2, 1, 1, 0.5), + ( 8, 200, 162, 9, 0.85631007588636132), + ( 2, 20, 11, 5, 0.39551083591331271), + ( 2, 5, 4, 3, 0.59999999999999987), + ( 16, 100, 71, 28, 0.050496322336354399), + ( 1, 5, 2, 2, 0.90000000000000002), + ( 0, 5, 4, 1, 0.20000000000000004), + ( 27, 200, 110, 54, 0.24032479119039216), + ( 0, 10, 2, 5, 0.22222222222222224), + ( 3, 50, 27, 5, 0.77138514980460271), + ( 2, 50, 24, 8, 0.15067269856977925), + ( 4, 20, 16, 7, 0.10113519091847264), + ( 13, 500, 408, 15, 0.79686197891279686), + ( 0, 5, 3, 1, 0.40000000000000008), + ( 0, 20, 14, 2, 0.078947368421052627), + ( 49, 100, 62, 79, 0.6077614986362827), + ( 2, 10, 3, 6, 0.83333333333333337), + ( 0, 50, 31, 2, 0.13959183673469389), + ( 2, 5, 4, 3, 0.59999999999999987), + ( 4, 50, 21, 8, 0.81380887468704521), + ( 0, 10, 7, 2, 0.066666666666666652), + ( 0, 10, 1, 4, 0.59999999999999987), + ( 0, 20, 4, 2, 0.63157894736842102), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 39, 500, 427, 51, 0.05047757656076568), + ( 2, 100, 6, 21, 0.89490672557682871), + ( 5, 20, 11, 9, 0.68904501071683733), + ( 0, 2, 1, 1, 0.5), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 14, 50, 27, 30, 0.16250719969887772), + ( 0, 5, 4, 1, 0.20000000000000004), + ( 0, 5, 4, 1, 0.20000000000000004), + ( 2, 10, 8, 4, 0.13333333333333333), + ( 1, 5, 3, 2, 0.69999999999999996), + ( 25, 100, 77, 31, 0.79699287800204943), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 7, 20, 15, 8, 0.94891640866873062), + ( 3, 50, 26, 7, 0.45339412360688952), + ( 1, 10, 8, 2, 0.37777777777777771), + ( 40, 200, 61, 134, 0.4491054454532335), + ( 1, 5, 2, 4, 0.40000000000000008), + ( 0, 10, 6, 1, 0.39999999999999991), + ( 1, 50, 10, 13, 0.19134773839560071), + ( 0, 2, 1, 1, 0.5), + ( 1, 20, 5, 2, 0.94736842105263153), + ( 7, 50, 12, 30, 0.57532691212157849), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 6, 10, 7, 9, 0.69999999999999996), + ( 0, 20, 2, 1, 0.90000000000000002), + ( 2, 10, 5, 3, 0.91666666666666663), + ( 0, 10, 8, 1, 0.19999999999999998), + (258, 500, 372, 347, 0.53219975096883698), + ( 1, 3, 2, 2, 0.66666666666666674), + ( 45, 200, 129, 68, 0.69415691010446789), + ( 1, 10, 8, 2, 0.37777777777777771), + ( 0, 10, 2, 1, 0.80000000000000004), + ( 1, 10, 4, 5, 0.26190476190476192), + ( 3, 50, 36, 4, 0.74422492401215801), + ( 0, 20, 6, 1, 0.69999999999999996), + ( 0, 5, 2, 3, 0.10000000000000002), + ( 1, 200, 47, 9, 0.33197417194852796), + ( 20, 50, 32, 30, 0.78323921453982637), + ( 16, 50, 21, 34, 0.9149336897131396), + ( 17, 50, 38, 22, 0.69599001425795692), + ( 0, 5, 2, 3, 0.10000000000000002), + ( 1, 5, 3, 2, 0.69999999999999996), + ( 0, 10, 9, 1, 0.10000000000000001), + ( 0, 5, 2, 3, 0.10000000000000002), + ( 2, 10, 5, 6, 0.26190476190476192), + ( 0, 5, 2, 1, 0.59999999999999987), + ( 7, 20, 16, 9, 0.62538699690402466), + ( 1, 100, 27, 2, 0.92909090909090908), + ( 27, 100, 58, 50, 0.271780848715515), + ( 47, 100, 96, 51, 0.063730084348641039), + ( 1, 20, 6, 2, 0.92105263157894735), + ( 1, 10, 6, 2, 0.66666666666666674), + ( 0, 2, 1, 1, 0.5), + ( 0, 20, 11, 1, 0.45000000000000001), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 0, 2, 1, 1, 0.5), + ( 0, 10, 1, 7, 0.29999999999999999), + ( 0, 2, 1, 1, 0.5), + ( 0, 100, 36, 1, 0.64000000000000001), + ( 1, 100, 68, 2, 0.53979797979797983), + ( 13, 200, 79, 29, 0.80029860188814683), + ( 0, 10, 5, 1, 0.49999999999999994), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 13, 100, 64, 21, 0.5065368728909565), + ( 1, 10, 6, 4, 0.11904761904761905), + ( 0, 2, 1, 1, 0.5), + ( 0, 5, 1, 2, 0.59999999999999987), + ( 0, 2, 1, 1, 0.5), + ( 1, 5, 4, 2, 0.40000000000000008), + ( 14, 50, 41, 17, 0.65850372332742224), + ( 0, 2, 1, 1, 0.5), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 1, 100, 2, 62, 0.61797979797979785), + ( 0, 2, 1, 1, 0.5), + ( 0, 2, 1, 1, 0.5), + ( 12, 500, 312, 16, 0.91020698917397613), + ( 0, 20, 2, 6, 0.47894736842105257), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 1, 10, 3, 4, 0.66666666666666674), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 6, 50, 20, 14, 0.72026241648862666), + ( 3, 20, 14, 6, 0.22523219814241485), + ( 0, 2, 1, 1, 0.5), + ( 4, 100, 72, 7, 0.30429108474790234), + ( 0, 5, 1, 2, 0.59999999999999987), + ( 0, 10, 4, 1, 0.59999999999999998), + ( 1, 3, 2, 2, 0.66666666666666674), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 22, 50, 46, 24, 0.66413373860182379), + ( 1, 5, 2, 4, 0.40000000000000008), + ( 62, 100, 80, 79, 0.3457586020522983), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 0, 10, 2, 7, 0.066666666666666666), + ( 0, 2, 1, 1, 0.5), + ( 0, 5, 2, 1, 0.59999999999999987), + ( 42, 200, 145, 57, 0.65622325663713577), + ( 1, 20, 12, 3, 0.34385964912280703), + ( 0, 2, 1, 1, 0.5), + ( 2, 10, 4, 7, 0.33333333333333331), + ( 1, 5, 3, 2, 0.69999999999999996), + ( 0, 10, 6, 2, 0.1333333333333333), + ( 2, 10, 6, 5, 0.26190476190476192), + ( 0, 5, 2, 1, 0.59999999999999987), + ( 1, 3, 2, 2, 0.66666666666666674), + ( 0, 50, 25, 2, 0.24489795918367349), + ( 0, 50, 39, 1, 0.22), + ( 2, 5, 3, 3, 0.90000000000000002), + ( 9, 50, 46, 10, 0.60316977854971765), + ( 0, 5, 2, 1, 0.59999999999999987), + ( 72, 500, 324, 112, 0.49074275180525029), + ( 0, 50, 9, 7, 0.22507959200836167), + ( 0, 5, 2, 2, 0.30000000000000004), + ( 17, 100, 35, 60, 0.067474411926413541), + ( 15, 100, 83, 17, 0.83718038506483827), + ( 0, 10, 7, 1, 0.29999999999999999), + ( 28, 200, 87, 77, 0.071226044946921765), + (154, 500, 361, 212, 0.61327756805578304), + ( 1, 10, 2, 3, 0.93333333333333335), + ( 0, 10, 4, 4, 0.071428571428571425), + ( 0, 5, 1, 1, 0.79999999999999993), + ( 2, 5, 3, 4, 0.59999999999999987), + ( 0, 10, 4, 1, 0.59999999999999998), + ( 0, 3, 2, 1, 0.33333333333333331), + ( 0, 10, 3, 1, 0.69999999999999996), + ( 0, 50, 10, 1, 0.80000000000000004), + ( 0, 2, 1, 1, 0.5), + ( 0, 10, 1, 3, 0.69999999999999996), + ( 2, 20, 12, 4, 0.53457172342621262), + ( 0, 5, 4, 1, 0.20000000000000004), + ( 4, 20, 9, 7, 0.89821981424148611), + ( 2, 200, 188, 3, 0.17021775544388609), + (132, 500, 298, 215, 0.78880271135040059), + ( 2, 5, 4, 3, 0.59999999999999987), + ( 0, 2, 1, 1, 0.5), + ( 2, 10, 6, 5, 0.26190476190476192), + ( 0, 3, 1, 1, 0.66666666666666674), + (156, 200, 128, 174, 1), + ( 1, 20, 6, 4, 0.65737874097007221), + ( 0, 5, 0, 0, 1), + (488, 500, 198, 500, 1), + (143, 500, 8, 371, 1), + ( 2, 10, 6, 5, 0.26190476190476192), + ( 1, 5, 2, 4, 0.40000000000000008), + ( 0, 3, 2, 0, 1), + ( 12, 50, 7, 17, 1), + (129, 200, 43, 133, 1), + ( 0, 5, 3, 0, 1), + ( 0, 2, 1, 1, 0.5), + ( 5, 20, 20, 17, 0), + ( 4, 10, 4, 8, 1), + ( 46, 500, 478, 58, 5.1715118817799218e-07), + ( 0, 3, 2, 3, 0), + ( 0, 3, 1, 1, 0.66666666666666674), + ( 76, 500, 0, 120, 1), + ( 1, 100, 41, 12, 0.011989696504564528), + ] + for (k, population, successes, draws, val) in tests: + check: abs(hypergeom_cdf(k, population, successes, draws) - val) < 1e-11 \ No newline at end of file