test_prob_distr.c 45 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456
  1. /* Copyright (c) 2018-2019, The Tor Project, Inc. */
  2. /* See LICENSE for licensing information */
  3. /**
  4. * \file test_prob_distr.c
  5. * \brief Test probability distributions.
  6. * \detail
  7. *
  8. * For each probability distribution we do two kinds of tests:
  9. *
  10. * a) We do numerical deterministic testing of their cdf/icdf/sf/isf functions
  11. * and the various relationships between them for each distribution. We also
  12. * do deterministic tests on their sampling functions. Test vectors for
  13. * these tests were computed from alternative implementations and were
  14. * eyeballed to make sure they make sense
  15. * (e.g. src/test/prob_distr_mpfr_ref.c computes logit(p) using GNU mpfr
  16. * with 200-bit precision and is then tested in test_logit_logistic()).
  17. *
  18. * b) We do stochastic hypothesis testing (G-test) to ensure that sampling from
  19. * the given distributions is distributed properly. The stochastic tests are
  20. * slow and their false positive rate is not well suited for CI, so they are
  21. * currently disabled-by-default and put into 'tests-slow'.
  22. */
  23. #define PROB_DISTR_PRIVATE
  24. #include "orconfig.h"
  25. #include "test/test.h"
  26. #include "core/or/or.h"
  27. #include "lib/math/prob_distr.h"
  28. #include "lib/math/fp.h"
  29. #include "lib/crypt_ops/crypto_rand.h"
  30. #include <float.h>
  31. #include <math.h>
  32. #include <stdbool.h>
  33. #include <stddef.h>
  34. #include <stdint.h>
  35. #include <stdio.h>
  36. #include <stdlib.h>
  37. /**
  38. * Return floor(d) converted to size_t, as a workaround for complaints
  39. * under -Wbad-function-cast for (size_t)floor(d).
  40. */
  41. static size_t
  42. floor_to_size_t(double d)
  43. {
  44. double integral_d = floor(d);
  45. return (size_t)integral_d;
  46. }
  47. /**
  48. * Return ceil(d) converted to size_t, as a workaround for complaints
  49. * under -Wbad-function-cast for (size_t)ceil(d).
  50. */
  51. static size_t
  52. ceil_to_size_t(double d)
  53. {
  54. double integral_d = ceil(d);
  55. return (size_t)integral_d;
  56. }
  57. /*
  58. * Geometric(p) distribution, supported on {1, 2, 3, ...}.
  59. *
  60. * Compute the probability mass function Geom(n; p) of the number of
  61. * trials before the first success when success has probability p.
  62. */
  63. static double
  64. logpmf_geometric(unsigned n, double p)
  65. {
  66. /* This is actually a check against 1, but we do >= so that the compiler
  67. does not raise a -Wfloat-equal */
  68. if (p >= 1) {
  69. if (n == 1)
  70. return 0;
  71. else
  72. return -HUGE_VAL;
  73. }
  74. return (n - 1)*log1p(-p) + log(p);
  75. }
  76. /**
  77. * Compute the logistic function, translated in output by 1/2:
  78. * logistichalf(x) = logistic(x) - 1/2. Well-conditioned on the entire
  79. * real plane, with maximum condition number 1 at 0.
  80. *
  81. * This implementation gives relative error bounded by 5 eps.
  82. */
  83. static double
  84. logistichalf(double x)
  85. {
  86. /*
  87. * Rewrite this with the identity
  88. *
  89. * 1/(1 + e^{-x}) - 1/2
  90. * = (1 - 1/2 - e^{-x}/2)/(1 + e^{-x})
  91. * = (1/2 - e^{-x}/2)/(1 + e^{-x})
  92. * = (1 - e^{-x})/[2 (1 + e^{-x})]
  93. * = -(e^{-x} - 1)/[2 (1 + e^{-x})],
  94. *
  95. * which we can evaluate by -expm1(-x)/[2 (1 + exp(-x))].
  96. *
  97. * Suppose exp has error d0, + has error d1, expm1 has error
  98. * d2, and / has error d3, so we evaluate
  99. *
  100. * -(1 + d2) (1 + d3) (e^{-x} - 1)
  101. * / [2 (1 + d1) (1 + (1 + d0) e^{-x})].
  102. *
  103. * In the denominator,
  104. *
  105. * 1 + (1 + d0) e^{-x}
  106. * = 1 + e^{-x} + d0 e^{-x}
  107. * = (1 + e^{-x}) (1 + d0 e^{-x}/(1 + e^{-x})),
  108. *
  109. * so the relative error of the numerator is
  110. *
  111. * d' = d2 + d3 + d2 d3,
  112. * and of the denominator,
  113. * d'' = d1 + d0 e^{-x}/(1 + e^{-x}) + d0 d1 e^{-x}/(1 + e^{-x})
  114. * = d1 + d0 L(-x) + d0 d1 L(-x),
  115. *
  116. * where L(-x) is logistic(-x). By Lemma 1 the relative error
  117. * of the quotient is bounded by
  118. *
  119. * 2|d2 + d3 + d2 d3 - d1 - d0 L(x) + d0 d1 L(x)|,
  120. *
  121. * Since 0 < L(x) < 1, this is bounded by
  122. *
  123. * 2|d2| + 2|d3| + 2|d2 d3| + 2|d1| + 2|d0| + 2|d0 d1|
  124. * <= 4 eps + 2 eps^2.
  125. */
  126. if (x < log(DBL_EPSILON/8)) {
  127. /*
  128. * Avoid overflow in e^{-x}. When x < log(eps/4), we
  129. * we further have x < logit(eps/4), so that
  130. * logistic(x) < eps/4. Hence the relative error of
  131. * logistic(x) - 1/2 from -1/2 is bounded by eps/2, and
  132. * so the relative error of -1/2 from logistic(x) - 1/2
  133. * is bounded by eps.
  134. */
  135. return -0.5;
  136. } else {
  137. return -expm1(-x)/(2*(1 + exp(-x)));
  138. }
  139. }
  140. /**
  141. * Compute the log of the sum of the exps. Caller should arrange the
  142. * array in descending order to minimize error because I don't want to
  143. * deal with using temporary space and the one caller in this file
  144. * arranges that anyway.
  145. *
  146. * Warning: This implementation does not handle infinite or NaN inputs
  147. * sensibly, because I don't need that here at the moment. (NaN, or
  148. * -inf and +inf together, should yield NaN; +inf and finite should
  149. * yield +inf; otherwise all -inf should be ignored because exp(-inf) =
  150. * 0.)
  151. */
  152. static double
  153. logsumexp(double *A, size_t n)
  154. {
  155. double maximum, sum;
  156. size_t i;
  157. if (n == 0)
  158. return log(0);
  159. maximum = A[0];
  160. for (i = 1; i < n; i++) {
  161. if (A[i] > maximum)
  162. maximum = A[i];
  163. }
  164. sum = 0;
  165. for (i = n; i --> 0;)
  166. sum += exp(A[i] - maximum);
  167. return log(sum) + maximum;
  168. }
  169. /**
  170. * Compute log(1 - e^x). Defined only for negative x so that e^x < 1.
  171. * This is the complement of a probability in log space.
  172. */
  173. static double
  174. log1mexp(double x)
  175. {
  176. /*
  177. * We want to compute log on [0, 1/2) but log1p on [1/2, +inf),
  178. * so partition x at -log(2) = log(1/2).
  179. */
  180. if (-log(2) < x)
  181. return log(-expm1(x));
  182. else
  183. return log1p(-exp(x));
  184. }
  185. /*
  186. * Tests of numerical errors in computing logit, logistic, and the
  187. * various cdfs, sfs, icdfs, and isfs.
  188. */
  189. #define arraycount(A) (sizeof(A)/sizeof(A[0]))
  190. /** Return relative error between <b>actual</b> and <b>expected</b>.
  191. * Special cases: If <b>expected</b> is zero or infinite, return 1 if
  192. * <b>actual</b> is equal to <b>expected</b> and 0 if not, since the
  193. * usual notion of relative error is undefined but we only use this
  194. * for testing relerr(e, a) <= bound. If either is NaN, return NaN,
  195. * which has the property that NaN <= bound is false no matter what
  196. * bound is.
  197. *
  198. * Beware: if you test !(relerr(e, a) > bound), then then the result
  199. * is true when a is NaN because NaN > bound is false too. See
  200. * CHECK_RELERR for correct use to decide when to report failure.
  201. */
  202. static double
  203. relerr(double expected, double actual)
  204. {
  205. /*
  206. * To silence -Wfloat-equal, we have to test for equality using
  207. * inequalities: we have (fabs(expected) <= 0) iff (expected == 0),
  208. * and (actual <= expected && actual >= expected) iff actual ==
  209. * expected whether expected is zero or infinite.
  210. */
  211. if (fabs(expected) <= 0 || tor_isinf(expected)) {
  212. if (actual <= expected && actual >= expected)
  213. return 0;
  214. else
  215. return 1;
  216. } else {
  217. return fabs((expected - actual)/expected);
  218. }
  219. }
  220. /** Check that relative error of <b>expected</b> and <b>actual</b> is within
  221. * <b>relerr_bound</b>. Caller must arrange to have i and relerr_bound in
  222. * scope. */
  223. #define CHECK_RELERR(expected, actual) do { \
  224. double check_expected = (expected); \
  225. double check_actual = (actual); \
  226. const char *str_expected = #expected; \
  227. const char *str_actual = #actual; \
  228. double check_relerr = relerr(expected, actual); \
  229. if (!(relerr(check_expected, check_actual) <= relerr_bound)) { \
  230. log_warn(LD_GENERAL, "%s:%d: case %u: relerr(%s=%.17e, %s=%.17e)" \
  231. " = %.17e > %.17e\n", \
  232. __func__, __LINE__, (unsigned) i, \
  233. str_expected, check_expected, \
  234. str_actual, check_actual, \
  235. check_relerr, relerr_bound); \
  236. ok = false; \
  237. } \
  238. } while (0)
  239. /* Check that a <= b.
  240. * Caller must arrange to have i in scope. */
  241. #define CHECK_LE(a, b) do { \
  242. double check_a = (a); \
  243. double check_b = (b); \
  244. const char *str_a = #a; \
  245. const char *str_b = #b; \
  246. if (!(check_a <= check_b)) { \
  247. log_warn(LD_GENERAL, "%s:%d: case %u: %s=%.17e > %s=%.17e\n", \
  248. __func__, __LINE__, (unsigned) i, \
  249. str_a, check_a, str_b, check_b); \
  250. ok = false; \
  251. } \
  252. } while (0)
  253. /**
  254. * Test the logit and logistic functions. Confirm that they agree with
  255. * the cdf, sf, icdf, and isf of the standard Logistic distribution.
  256. * Confirm that the sampler for the standard logistic distribution maps
  257. * [0, 1] into the right subinterval for the inverse transform, for
  258. * this implementation.
  259. */
  260. static void
  261. test_logit_logistic(void *arg)
  262. {
  263. (void) arg;
  264. static const struct {
  265. double x; /* x = logit(p) */
  266. double p; /* p = logistic(x) */
  267. double phalf; /* p - 1/2 = logistic(x) - 1/2 */
  268. } cases[] = {
  269. { -HUGE_VAL, 0, -0.5 },
  270. { -1000, 0, -0.5 },
  271. { -710, 4.47628622567513e-309, -0.5 },
  272. { -708, 3.307553003638408e-308, -0.5 },
  273. { -2, .11920292202211755, -.3807970779778824 },
  274. { -1.0000001, .2689414017088022, -.23105859829119776 },
  275. { -1, .2689414213699951, -.23105857863000487 },
  276. { -0.9999999, .26894144103118883, -.2310585589688111 },
  277. /* see src/test/prob_distr_mpfr_ref.c for computation */
  278. { -4.000000000537333e-5, .49999, -1.0000000000010001e-5 },
  279. { -4.000000000533334e-5, .49999, -.00001 },
  280. { -4.000000108916878e-9, .499999999, -1.0000000272292198e-9 },
  281. { -4e-9, .499999999, -1e-9 },
  282. { -4e-16, .5, -1e-16 },
  283. { -4e-300, .5, -1e-300 },
  284. { 0, .5, 0 },
  285. { 4e-300, .5, 1e-300 },
  286. { 4e-16, .5, 1e-16 },
  287. { 3.999999886872274e-9, .500000001, 9.999999717180685e-10 },
  288. { 4e-9, .500000001, 1e-9 },
  289. { 4.0000000005333336e-5, .50001, .00001 },
  290. { 8.000042667076272e-3, .502, .002 },
  291. { 0.9999999, .7310585589688111, .2310585589688111 },
  292. { 1, .7310585786300049, .23105857863000487 },
  293. { 1.0000001, .7310585982911977, .23105859829119774 },
  294. { 2, .8807970779778823, .3807970779778824 },
  295. { 708, 1, .5 },
  296. { 710, 1, .5 },
  297. { 1000, 1, .5 },
  298. { HUGE_VAL, 1, .5 },
  299. };
  300. double relerr_bound = 3e-15; /* >10eps */
  301. size_t i;
  302. bool ok = true;
  303. for (i = 0; i < arraycount(cases); i++) {
  304. double x = cases[i].x;
  305. double p = cases[i].p;
  306. double phalf = cases[i].phalf;
  307. /*
  308. * cdf is logistic, icdf is logit, and symmetry for
  309. * sf/isf.
  310. */
  311. CHECK_RELERR(logistic(x), cdf_logistic(x, 0, 1));
  312. CHECK_RELERR(logistic(-x), sf_logistic(x, 0, 1));
  313. CHECK_RELERR(logit(p), icdf_logistic(p, 0, 1));
  314. CHECK_RELERR(-logit(p), isf_logistic(p, 0, 1));
  315. CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x*2, 0, 2));
  316. CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x*2, 0, 2));
  317. CHECK_RELERR(icdf_logistic(p, 0, 1), icdf_logistic(p, 0, 2)/2);
  318. CHECK_RELERR(isf_logistic(p, 0, 1), isf_logistic(p, 0, 2)/2);
  319. CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x/2, 0, .5));
  320. CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x/2, 0, .5));
  321. CHECK_RELERR(icdf_logistic(p, 0, 1), icdf_logistic(p, 0,.5)*2);
  322. CHECK_RELERR(isf_logistic(p, 0, 1), isf_logistic(p, 0, .5)*2);
  323. CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x*2 + 1, 1, 2));
  324. CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x*2 + 1, 1, 2));
  325. /*
  326. * For p near 0 and p near 1/2, the arithmetic of
  327. * translating by 1 loses precision.
  328. */
  329. if (fabs(p) > DBL_EPSILON && fabs(p) < 0.4) {
  330. CHECK_RELERR(icdf_logistic(p, 0, 1),
  331. (icdf_logistic(p, 1, 2) - 1)/2);
  332. CHECK_RELERR(isf_logistic(p, 0, 1),
  333. (isf_logistic(p, 1, 2) - 1)/2);
  334. }
  335. CHECK_RELERR(p, logistic(x));
  336. CHECK_RELERR(phalf, logistichalf(x));
  337. /*
  338. * On the interior floating-point numbers, either logit or
  339. * logithalf had better give the correct answer.
  340. *
  341. * For probabilities near 0, we can get much finer resolution with
  342. * logit, and for probabilities near 1/2, we can get much finer
  343. * resolution with logithalf by representing them using p - 1/2.
  344. *
  345. * E.g., we can write -.00001 for phalf, and .49999 for p, but the
  346. * difference 1/2 - .00001 gives 1.0000000000010001e-5 in binary64
  347. * arithmetic. So test logit(.49999) which should give the same
  348. * answer as logithalf(-1.0000000000010001e-5), namely
  349. * -4.000000000537333e-5, and also test logithalf(-.00001) which
  350. * gives -4.000000000533334e-5 instead -- but don't expect
  351. * logit(.49999) to give -4.000000000533334e-5 even though it looks
  352. * like 1/2 - .00001.
  353. *
  354. * A naive implementation of logit will just use log(p/(1 - p)) and
  355. * give the answer -4.000000000551673e-05 for .49999, which is
  356. * wrong in a lot of digits, which happens because log is
  357. * ill-conditioned near 1 and thus amplifies whatever relative
  358. * error we made in computing p/(1 - p).
  359. */
  360. if ((0 < p && p < 1) || tor_isinf(x)) {
  361. if (phalf >= p - 0.5 && phalf <= p - 0.5)
  362. CHECK_RELERR(x, logit(p));
  363. if (p >= 0.5 + phalf && p <= 0.5 + phalf)
  364. CHECK_RELERR(x, logithalf(phalf));
  365. }
  366. CHECK_RELERR(-phalf, logistichalf(-x));
  367. if (fabs(phalf) < 0.5 || tor_isinf(x))
  368. CHECK_RELERR(-x, logithalf(-phalf));
  369. if (p < 1 || tor_isinf(x)) {
  370. CHECK_RELERR(1 - p, logistic(-x));
  371. if (p > .75 || tor_isinf(x))
  372. CHECK_RELERR(-x, logit(1 - p));
  373. } else {
  374. CHECK_LE(logistic(-x), 1e-300);
  375. }
  376. }
  377. for (i = 0; i <= 100; i++) {
  378. double p0 = (double)i/100;
  379. CHECK_RELERR(logit(p0/(1 + M_E)), sample_logistic(0, 0, p0));
  380. CHECK_RELERR(-logit(p0/(1 + M_E)), sample_logistic(1, 0, p0));
  381. CHECK_RELERR(logithalf(p0*(0.5 - 1/(1 + M_E))),
  382. sample_logistic(0, 1, p0));
  383. CHECK_RELERR(-logithalf(p0*(0.5 - 1/(1 + M_E))),
  384. sample_logistic(1, 1, p0));
  385. }
  386. if (!ok)
  387. printf("fail logit/logistic / logistic cdf/sf\n");
  388. tt_assert(ok);
  389. done:
  390. ;
  391. }
  392. /**
  393. * Test the cdf, sf, icdf, and isf of the LogLogistic distribution.
  394. */
  395. static void
  396. test_log_logistic(void *arg)
  397. {
  398. (void) arg;
  399. static const struct {
  400. /* x is a point in the support of the LogLogistic distribution */
  401. double x;
  402. /* 'p' is the probability that a random variable X for a given LogLogistic
  403. * probability ditribution will take value less-or-equal to x */
  404. double p;
  405. /* 'np' is the probability that a random variable X for a given LogLogistic
  406. * probability distribution will take value greater-or-equal to x. */
  407. double np;
  408. } cases[] = {
  409. { 0, 0, 1 },
  410. { 1e-300, 1e-300, 1 },
  411. { 1e-17, 1e-17, 1 },
  412. { 1e-15, 1e-15, .999999999999999 },
  413. { .1, .09090909090909091, .90909090909090909 },
  414. { .25, .2, .8 },
  415. { .5, .33333333333333333, .66666666666666667 },
  416. { .75, .42857142857142855, .5714285714285714 },
  417. { .9999, .49997499874993756, .5000250012500626 },
  418. { .99999999, .49999999749999996, .5000000025 },
  419. { .999999999999999, .49999999999999994, .5000000000000002 },
  420. { 1, .5, .5 },
  421. };
  422. double relerr_bound = 3e-15;
  423. size_t i;
  424. bool ok = true;
  425. for (i = 0; i < arraycount(cases); i++) {
  426. double x = cases[i].x;
  427. double p = cases[i].p;
  428. double np = cases[i].np;
  429. CHECK_RELERR(p, cdf_log_logistic(x, 1, 1));
  430. CHECK_RELERR(p, cdf_log_logistic(x/2, .5, 1));
  431. CHECK_RELERR(p, cdf_log_logistic(x*2, 2, 1));
  432. CHECK_RELERR(p, cdf_log_logistic(sqrt(x), 1, 2));
  433. CHECK_RELERR(p, cdf_log_logistic(sqrt(x)/2, .5, 2));
  434. CHECK_RELERR(p, cdf_log_logistic(sqrt(x)*2, 2, 2));
  435. if (2*sqrt(DBL_MIN) < x) {
  436. CHECK_RELERR(p, cdf_log_logistic(x*x, 1, .5));
  437. CHECK_RELERR(p, cdf_log_logistic(x*x/2, .5, .5));
  438. CHECK_RELERR(p, cdf_log_logistic(x*x*2, 2, .5));
  439. }
  440. CHECK_RELERR(np, sf_log_logistic(x, 1, 1));
  441. CHECK_RELERR(np, sf_log_logistic(x/2, .5, 1));
  442. CHECK_RELERR(np, sf_log_logistic(x*2, 2, 1));
  443. CHECK_RELERR(np, sf_log_logistic(sqrt(x), 1, 2));
  444. CHECK_RELERR(np, sf_log_logistic(sqrt(x)/2, .5, 2));
  445. CHECK_RELERR(np, sf_log_logistic(sqrt(x)*2, 2, 2));
  446. if (2*sqrt(DBL_MIN) < x) {
  447. CHECK_RELERR(np, sf_log_logistic(x*x, 1, .5));
  448. CHECK_RELERR(np, sf_log_logistic(x*x/2, .5, .5));
  449. CHECK_RELERR(np, sf_log_logistic(x*x*2, 2, .5));
  450. }
  451. CHECK_RELERR(np, cdf_log_logistic(1/x, 1, 1));
  452. CHECK_RELERR(np, cdf_log_logistic(1/(2*x), .5, 1));
  453. CHECK_RELERR(np, cdf_log_logistic(2/x, 2, 1));
  454. CHECK_RELERR(np, cdf_log_logistic(1/sqrt(x), 1, 2));
  455. CHECK_RELERR(np, cdf_log_logistic(1/(2*sqrt(x)), .5, 2));
  456. CHECK_RELERR(np, cdf_log_logistic(2/sqrt(x), 2, 2));
  457. if (2*sqrt(DBL_MIN) < x && x < 1/(2*sqrt(DBL_MIN))) {
  458. CHECK_RELERR(np, cdf_log_logistic(1/(x*x), 1, .5));
  459. CHECK_RELERR(np, cdf_log_logistic(1/(2*x*x), .5, .5));
  460. CHECK_RELERR(np, cdf_log_logistic(2/(x*x), 2, .5));
  461. }
  462. CHECK_RELERR(p, sf_log_logistic(1/x, 1, 1));
  463. CHECK_RELERR(p, sf_log_logistic(1/(2*x), .5, 1));
  464. CHECK_RELERR(p, sf_log_logistic(2/x, 2, 1));
  465. CHECK_RELERR(p, sf_log_logistic(1/sqrt(x), 1, 2));
  466. CHECK_RELERR(p, sf_log_logistic(1/(2*sqrt(x)), .5, 2));
  467. CHECK_RELERR(p, sf_log_logistic(2/sqrt(x), 2, 2));
  468. if (2*sqrt(DBL_MIN) < x && x < 1/(2*sqrt(DBL_MIN))) {
  469. CHECK_RELERR(p, sf_log_logistic(1/(x*x), 1, .5));
  470. CHECK_RELERR(p, sf_log_logistic(1/(2*x*x), .5, .5));
  471. CHECK_RELERR(p, sf_log_logistic(2/(x*x), 2, .5));
  472. }
  473. CHECK_RELERR(x, icdf_log_logistic(p, 1, 1));
  474. CHECK_RELERR(x/2, icdf_log_logistic(p, .5, 1));
  475. CHECK_RELERR(x*2, icdf_log_logistic(p, 2, 1));
  476. CHECK_RELERR(x, icdf_log_logistic(p, 1, 1));
  477. CHECK_RELERR(sqrt(x)/2, icdf_log_logistic(p, .5, 2));
  478. CHECK_RELERR(sqrt(x)*2, icdf_log_logistic(p, 2, 2));
  479. CHECK_RELERR(sqrt(x), icdf_log_logistic(p, 1, 2));
  480. CHECK_RELERR(x*x/2, icdf_log_logistic(p, .5, .5));
  481. CHECK_RELERR(x*x*2, icdf_log_logistic(p, 2, .5));
  482. if (np < .9) {
  483. CHECK_RELERR(x, isf_log_logistic(np, 1, 1));
  484. CHECK_RELERR(x/2, isf_log_logistic(np, .5, 1));
  485. CHECK_RELERR(x*2, isf_log_logistic(np, 2, 1));
  486. CHECK_RELERR(sqrt(x), isf_log_logistic(np, 1, 2));
  487. CHECK_RELERR(sqrt(x)/2, isf_log_logistic(np, .5, 2));
  488. CHECK_RELERR(sqrt(x)*2, isf_log_logistic(np, 2, 2));
  489. CHECK_RELERR(x*x, isf_log_logistic(np, 1, .5));
  490. CHECK_RELERR(x*x/2, isf_log_logistic(np, .5, .5));
  491. CHECK_RELERR(x*x*2, isf_log_logistic(np, 2, .5));
  492. CHECK_RELERR(1/x, icdf_log_logistic(np, 1, 1));
  493. CHECK_RELERR(1/(2*x), icdf_log_logistic(np, .5, 1));
  494. CHECK_RELERR(2/x, icdf_log_logistic(np, 2, 1));
  495. CHECK_RELERR(1/sqrt(x), icdf_log_logistic(np, 1, 2));
  496. CHECK_RELERR(1/(2*sqrt(x)),
  497. icdf_log_logistic(np, .5, 2));
  498. CHECK_RELERR(2/sqrt(x), icdf_log_logistic(np, 2, 2));
  499. CHECK_RELERR(1/(x*x), icdf_log_logistic(np, 1, .5));
  500. CHECK_RELERR(1/(2*x*x), icdf_log_logistic(np, .5, .5));
  501. CHECK_RELERR(2/(x*x), icdf_log_logistic(np, 2, .5));
  502. }
  503. CHECK_RELERR(1/x, isf_log_logistic(p, 1, 1));
  504. CHECK_RELERR(1/(2*x), isf_log_logistic(p, .5, 1));
  505. CHECK_RELERR(2/x, isf_log_logistic(p, 2, 1));
  506. CHECK_RELERR(1/sqrt(x), isf_log_logistic(p, 1, 2));
  507. CHECK_RELERR(1/(2*sqrt(x)), isf_log_logistic(p, .5, 2));
  508. CHECK_RELERR(2/sqrt(x), isf_log_logistic(p, 2, 2));
  509. CHECK_RELERR(1/(x*x), isf_log_logistic(p, 1, .5));
  510. CHECK_RELERR(1/(2*x*x), isf_log_logistic(p, .5, .5));
  511. CHECK_RELERR(2/(x*x), isf_log_logistic(p, 2, .5));
  512. }
  513. for (i = 0; i <= 100; i++) {
  514. double p0 = (double)i/100;
  515. CHECK_RELERR(0.5*p0/(1 - 0.5*p0), sample_log_logistic(0, p0));
  516. CHECK_RELERR((1 - 0.5*p0)/(0.5*p0),
  517. sample_log_logistic(1, p0));
  518. }
  519. if (!ok)
  520. printf("fail log logistic cdf/sf\n");
  521. tt_assert(ok);
  522. done:
  523. ;
  524. }
  525. /**
  526. * Test the cdf, sf, icdf, isf of the Weibull distribution.
  527. */
  528. static void
  529. test_weibull(void *arg)
  530. {
  531. (void) arg;
  532. static const struct {
  533. /* x is a point in the support of the Weibull distribution */
  534. double x;
  535. /* 'p' is the probability that a random variable X for a given Weibull
  536. * probability ditribution will take value less-or-equal to x */
  537. double p;
  538. /* 'np' is the probability that a random variable X for a given Weibull
  539. * probability distribution will take value greater-or-equal to x. */
  540. double np;
  541. } cases[] = {
  542. { 0, 0, 1 },
  543. { 1e-300, 1e-300, 1 },
  544. { 1e-17, 1e-17, 1 },
  545. { .1, .09516258196404043, .9048374180359595 },
  546. { .5, .3934693402873666, .6065306597126334 },
  547. { .6931471805599453, .5, .5 },
  548. { 1, .6321205588285577, .36787944117144233 },
  549. { 10, .9999546000702375, 4.5399929762484854e-5 },
  550. { 36, .9999999999999998, 2.319522830243569e-16 },
  551. { 37, .9999999999999999, 8.533047625744066e-17 },
  552. { 38, 1, 3.1391327920480296e-17 },
  553. { 100, 1, 3.720075976020836e-44 },
  554. { 708, 1, 3.307553003638408e-308 },
  555. { 710, 1, 4.47628622567513e-309 },
  556. { 1000, 1, 0 },
  557. { HUGE_VAL, 1, 0 },
  558. };
  559. double relerr_bound = 3e-15;
  560. size_t i;
  561. bool ok = true;
  562. for (i = 0; i < arraycount(cases); i++) {
  563. double x = cases[i].x;
  564. double p = cases[i].p;
  565. double np = cases[i].np;
  566. CHECK_RELERR(p, cdf_weibull(x, 1, 1));
  567. CHECK_RELERR(p, cdf_weibull(x/2, .5, 1));
  568. CHECK_RELERR(p, cdf_weibull(x*2, 2, 1));
  569. /* For 0 < x < sqrt(DBL_MIN), x^2 loses lots of bits. */
  570. if (x <= 0 ||
  571. sqrt(DBL_MIN) <= x) {
  572. CHECK_RELERR(p, cdf_weibull(x*x, 1, .5));
  573. CHECK_RELERR(p, cdf_weibull(x*x/2, .5, .5));
  574. CHECK_RELERR(p, cdf_weibull(x*x*2, 2, .5));
  575. }
  576. CHECK_RELERR(p, cdf_weibull(sqrt(x), 1, 2));
  577. CHECK_RELERR(p, cdf_weibull(sqrt(x)/2, .5, 2));
  578. CHECK_RELERR(p, cdf_weibull(sqrt(x)*2, 2, 2));
  579. CHECK_RELERR(np, sf_weibull(x, 1, 1));
  580. CHECK_RELERR(np, sf_weibull(x/2, .5, 1));
  581. CHECK_RELERR(np, sf_weibull(x*2, 2, 1));
  582. CHECK_RELERR(np, sf_weibull(x*x, 1, .5));
  583. CHECK_RELERR(np, sf_weibull(x*x/2, .5, .5));
  584. CHECK_RELERR(np, sf_weibull(x*x*2, 2, .5));
  585. if (x >= 10) {
  586. /*
  587. * exp amplifies the error of sqrt(x)^2
  588. * proportionally to exp(x); for large inputs
  589. * this is significant.
  590. */
  591. double t = -expm1(-x*(2*DBL_EPSILON + DBL_EPSILON));
  592. relerr_bound = t + DBL_EPSILON + t*DBL_EPSILON;
  593. if (relerr_bound < 3e-15)
  594. /*
  595. * The tests are written only to 16
  596. * decimal places anyway even if your
  597. * `double' is, say, i387 binary80, for
  598. * whatever reason.
  599. */
  600. relerr_bound = 3e-15;
  601. CHECK_RELERR(np, sf_weibull(sqrt(x), 1, 2));
  602. CHECK_RELERR(np, sf_weibull(sqrt(x)/2, .5, 2));
  603. CHECK_RELERR(np, sf_weibull(sqrt(x)*2, 2, 2));
  604. }
  605. if (p <= 0.75) {
  606. /*
  607. * For p near 1, not enough precision near 1 to
  608. * recover x.
  609. */
  610. CHECK_RELERR(x, icdf_weibull(p, 1, 1));
  611. CHECK_RELERR(x/2, icdf_weibull(p, .5, 1));
  612. CHECK_RELERR(x*2, icdf_weibull(p, 2, 1));
  613. }
  614. if (p >= 0.25 && !tor_isinf(x) && np > 0) {
  615. /*
  616. * For p near 0, not enough precision in np
  617. * near 1 to recover x. For 0, isf gives inf,
  618. * even if p is precise enough for the icdf to
  619. * work.
  620. */
  621. CHECK_RELERR(x, isf_weibull(np, 1, 1));
  622. CHECK_RELERR(x/2, isf_weibull(np, .5, 1));
  623. CHECK_RELERR(x*2, isf_weibull(np, 2, 1));
  624. }
  625. }
  626. for (i = 0; i <= 100; i++) {
  627. double p0 = (double)i/100;
  628. CHECK_RELERR(3*sqrt(-log(p0/2)), sample_weibull(0, p0, 3, 2));
  629. CHECK_RELERR(3*sqrt(-log1p(-p0/2)),
  630. sample_weibull(1, p0, 3, 2));
  631. }
  632. if (!ok)
  633. printf("fail Weibull cdf/sf\n");
  634. tt_assert(ok);
  635. done:
  636. ;
  637. }
  638. /**
  639. * Test the cdf, sf, icdf, and isf of the generalized Pareto
  640. * distribution.
  641. */
  642. static void
  643. test_genpareto(void *arg)
  644. {
  645. (void) arg;
  646. struct {
  647. /* xi is the 'xi' parameter of the generalized Pareto distribution, and the
  648. * rest are the same as in the above tests */
  649. double xi, x, p, np;
  650. } cases[] = {
  651. { 0, 0, 0, 1 },
  652. { 1e-300, .004, 3.992010656008528e-3, .9960079893439915 },
  653. { 1e-300, .1, .09516258196404043, .9048374180359595 },
  654. { 1e-300, 1, .6321205588285577, .36787944117144233 },
  655. { 1e-300, 10, .9999546000702375, 4.5399929762484854e-5 },
  656. { 1e-200, 1e-16, 9.999999999999999e-17, .9999999999999999 },
  657. { 1e-16, 1e-200, 9.999999999999998e-201, 1 },
  658. { 1e-16, 1e-16, 1e-16, 1 },
  659. { 1e-16, .004, 3.992010656008528e-3, .9960079893439915 },
  660. { 1e-16, .1, .09516258196404043, .9048374180359595 },
  661. { 1e-16, 1, .6321205588285577, .36787944117144233 },
  662. { 1e-16, 10, .9999546000702375, 4.539992976248509e-5 },
  663. { 1e-10, 1e-6, 9.999995000001667e-7, .9999990000005 },
  664. { 1e-8, 1e-8, 9.999999950000001e-9, .9999999900000001 },
  665. { 1, 1e-300, 1e-300, 1 },
  666. { 1, 1e-16, 1e-16, .9999999999999999 },
  667. { 1, .1, .09090909090909091, .9090909090909091 },
  668. { 1, 1, .5, .5 },
  669. { 1, 10, .9090909090909091, .0909090909090909 },
  670. { 1, 100, .9900990099009901, .0099009900990099 },
  671. { 1, 1000, .999000999000999, 9.990009990009992e-4 },
  672. { 10, 1e-300, 1e-300, 1 },
  673. { 10, 1e-16, 9.999999999999995e-17, .9999999999999999 },
  674. { 10, .1, .06696700846319258, .9330329915368074 },
  675. { 10, 1, .21320655780322778, .7867934421967723 },
  676. { 10, 10, .3696701667040189, .6303298332959811 },
  677. { 10, 100, .49886285755007337, .5011371424499267 },
  678. { 10, 1000, .6018968102992647, .3981031897007353 },
  679. };
  680. double xi_array[] = { -1.5, -1, -1e-30, 0, 1e-30, 1, 1.5 };
  681. size_t i, j;
  682. double relerr_bound = 3e-15;
  683. bool ok = true;
  684. for (i = 0; i < arraycount(cases); i++) {
  685. double xi = cases[i].xi;
  686. double x = cases[i].x;
  687. double p = cases[i].p;
  688. double np = cases[i].np;
  689. CHECK_RELERR(p, cdf_genpareto(x, 0, 1, xi));
  690. CHECK_RELERR(p, cdf_genpareto(x*2, 0, 2, xi));
  691. CHECK_RELERR(p, cdf_genpareto(x/2, 0, .5, xi));
  692. CHECK_RELERR(np, sf_genpareto(x, 0, 1, xi));
  693. CHECK_RELERR(np, sf_genpareto(x*2, 0, 2, xi));
  694. CHECK_RELERR(np, sf_genpareto(x/2, 0, .5, xi));
  695. if (p < .5) {
  696. CHECK_RELERR(x, icdf_genpareto(p, 0, 1, xi));
  697. CHECK_RELERR(x*2, icdf_genpareto(p, 0, 2, xi));
  698. CHECK_RELERR(x/2, icdf_genpareto(p, 0, .5, xi));
  699. }
  700. if (np < .5) {
  701. CHECK_RELERR(x, isf_genpareto(np, 0, 1, xi));
  702. CHECK_RELERR(x*2, isf_genpareto(np, 0, 2, xi));
  703. CHECK_RELERR(x/2, isf_genpareto(np, 0, .5, xi));
  704. }
  705. }
  706. for (i = 0; i < arraycount(xi_array); i++) {
  707. for (j = 0; j <= 100; j++) {
  708. double p0 = (j == 0 ? 2*DBL_MIN : (double)j/100);
  709. /* This is actually a check against 0, but we do <= so that the compiler
  710. does not raise a -Wfloat-equal */
  711. if (fabs(xi_array[i]) <= 0) {
  712. /*
  713. * When xi == 0, the generalized Pareto
  714. * distribution reduces to an
  715. * exponential distribution.
  716. */
  717. CHECK_RELERR(-log(p0/2),
  718. sample_genpareto(0, p0, 0));
  719. CHECK_RELERR(-log1p(-p0/2),
  720. sample_genpareto(1, p0, 0));
  721. } else {
  722. CHECK_RELERR(expm1(-xi_array[i]*log(p0/2))/xi_array[i],
  723. sample_genpareto(0, p0, xi_array[i]));
  724. CHECK_RELERR((j == 0 ? DBL_MIN :
  725. expm1(-xi_array[i]*log1p(-p0/2))/xi_array[i]),
  726. sample_genpareto(1, p0, xi_array[i]));
  727. }
  728. CHECK_RELERR(isf_genpareto(p0/2, 0, 1, xi_array[i]),
  729. sample_genpareto(0, p0, xi_array[i]));
  730. CHECK_RELERR(icdf_genpareto(p0/2, 0, 1, xi_array[i]),
  731. sample_genpareto(1, p0, xi_array[i]));
  732. }
  733. }
  734. tt_assert(ok);
  735. done:
  736. ;
  737. }
  738. /**
  739. * Test the deterministic sampler for uniform distribution on [a, b].
  740. *
  741. * This currently only tests whether the outcome lies within [a, b].
  742. */
  743. static void
  744. test_uniform_interval(void *arg)
  745. {
  746. (void) arg;
  747. struct {
  748. /* Sample from a uniform distribution with parameters 'a' and 'b', using
  749. * 't' as the sampling index. */
  750. double t, a, b;
  751. } cases[] = {
  752. { 0, 0, 0 },
  753. { 0, 0, 1 },
  754. { 0, 1.0000000000000007, 3.999999999999995 },
  755. { 0, 4000, 4000 },
  756. { 0.42475836677491291, 4000, 4000 },
  757. { 0, -DBL_MAX, DBL_MAX },
  758. { 0.25, -DBL_MAX, DBL_MAX },
  759. { 0.5, -DBL_MAX, DBL_MAX },
  760. };
  761. size_t i = 0;
  762. bool ok = true;
  763. for (i = 0; i < arraycount(cases); i++) {
  764. double t = cases[i].t;
  765. double a = cases[i].a;
  766. double b = cases[i].b;
  767. CHECK_LE(a, sample_uniform_interval(t, a, b));
  768. CHECK_LE(sample_uniform_interval(t, a, b), b);
  769. CHECK_LE(a, sample_uniform_interval(1 - t, a, b));
  770. CHECK_LE(sample_uniform_interval(1 - t, a, b), b);
  771. CHECK_LE(sample_uniform_interval(t, -b, -a), -a);
  772. CHECK_LE(-b, sample_uniform_interval(t, -b, -a));
  773. CHECK_LE(sample_uniform_interval(1 - t, -b, -a), -a);
  774. CHECK_LE(-b, sample_uniform_interval(1 - t, -b, -a));
  775. }
  776. tt_assert(ok);
  777. done:
  778. ;
  779. }
  780. /********************** Stochastic tests ****************************/
  781. /*
  782. * Psi test, sometimes also called G-test. The psi test statistic,
  783. * suitably scaled, has chi^2 distribution, but the psi test tends to
  784. * have better statistical power in practice to detect deviations than
  785. * the chi^2 test does. (The chi^2 test statistic is the first term of
  786. * the Taylor expansion of the psi test statistic.) The psi test is
  787. * generic, for any CDF; particular distributions might have higher-
  788. * power tests to distinguish them from predictable deviations or bugs.
  789. *
  790. * We choose the psi critical value so that a single psi test has
  791. * probability below alpha = 1% of spuriously failing even if all the
  792. * code is correct. But the false positive rate for a suite of n tests
  793. * is higher: 1 - Binom(0; n, alpha) = 1 - (1 - alpha)^n. For n = 10,
  794. * this is about 10%, and for n = 100 it is well over 50%.
  795. *
  796. * Given that these tests will run with every CI job, we want to drive down the
  797. * false positive rate. We can drive it down by running each test four times,
  798. * and accepting it if it passes at least once; in that case, it is as if we
  799. * used Binom(4; 2, alpha) = alpha^4 as the false positive rate for each test,
  800. * and for n = 10 tests, it would be 9.99999959506e-08. If each CI build has 14
  801. * jobs, then the chance of a CI build failing is 1.39999903326e-06, which
  802. * means that a CI build will break with probability 50% after about 495106
  803. * builds.
  804. *
  805. * The critical value for a chi^2 distribution with 100 degrees of
  806. * freedom and false positive rate alpha = 1% was taken from:
  807. *
  808. * NIST/SEMATECH e-Handbook of Statistical Methods, Section
  809. * 1.3.6.7.4 `Critical Values of the Chi-Square Distribution',
  810. * <http://www.itl.nist.gov/div898/handbook/eda/section3/eda3674.htm>,
  811. * retrieved 2018-10-28.
  812. */
  813. static const size_t NSAMPLES = 100000;
  814. /* Number of chances we give to the test to succeed. */
  815. static const unsigned NTRIALS = 4;
  816. /* Number of times we want the test to pass per NTRIALS. */
  817. static const unsigned NPASSES_MIN = 1;
  818. #define PSI_DF 100 /* degrees of freedom */
  819. static const double PSI_CRITICAL = 135.807; /* critical value, alpha = .01 */
  820. /**
  821. * Perform a psi test on an array of sample counts, C, adding up to N
  822. * samples, and an array of log expected probabilities, logP,
  823. * representing the null hypothesis for the distribution of samples
  824. * counted. Return false if the psi test rejects the null hypothesis,
  825. * true if otherwise.
  826. */
  827. static bool
  828. psi_test(const size_t C[PSI_DF], const double logP[PSI_DF], size_t N)
  829. {
  830. double psi = 0;
  831. double c = 0; /* Kahan compensation */
  832. double t, u;
  833. size_t i;
  834. for (i = 0; i < PSI_DF; i++) {
  835. /*
  836. * c*log(c/(n*p)) = (1/n) * f*log(f/p) where f = c/n is
  837. * the frequency, and f*log(f/p) ---> 0 as f ---> 0, so
  838. * this is a reasonable choice. Further, any mass that
  839. * _fails_ to turn up in this bin will inflate another
  840. * bin instead, so we don't really lose anything by
  841. * ignoring empty bins even if they have high
  842. * probability.
  843. */
  844. if (C[i] == 0)
  845. continue;
  846. t = C[i]*(log((double)C[i]/N) - logP[i]) - c;
  847. u = psi + t;
  848. c = (u - psi) - t;
  849. psi = u;
  850. }
  851. psi *= 2;
  852. return psi <= PSI_CRITICAL;
  853. }
  854. static bool
  855. test_stochastic_geometric_impl(double p)
  856. {
  857. const struct geometric geometric = {
  858. .base = GEOMETRIC(geometric),
  859. .p = p,
  860. };
  861. double logP[PSI_DF] = {0};
  862. unsigned ntry = NTRIALS, npass = 0;
  863. unsigned i;
  864. size_t j;
  865. /* Compute logP[i] = Geom(i + 1; p). */
  866. for (i = 0; i < PSI_DF - 1; i++)
  867. logP[i] = logpmf_geometric(i + 1, p);
  868. /* Compute logP[n-1] = log (1 - (P[0] + P[1] + ... + P[n-2])). */
  869. logP[PSI_DF - 1] = log1mexp(logsumexp(logP, PSI_DF - 1));
  870. while (ntry --> 0) {
  871. size_t C[PSI_DF] = {0};
  872. for (j = 0; j < NSAMPLES; j++) {
  873. double n_tmp = dist_sample(&geometric.base);
  874. /* Must be an integer. (XXX -Wfloat-equal) */
  875. tor_assert(ceil(n_tmp) <= n_tmp && ceil(n_tmp) >= n_tmp);
  876. /* Must be a positive integer. */
  877. tor_assert(n_tmp >= 1);
  878. /* Probability of getting a value in the billions is negligible. */
  879. tor_assert(n_tmp <= (double)UINT_MAX);
  880. unsigned n = (unsigned) n_tmp;
  881. if (n > PSI_DF)
  882. n = PSI_DF;
  883. C[n - 1]++;
  884. }
  885. if (psi_test(C, logP, NSAMPLES)) {
  886. if (++npass >= NPASSES_MIN)
  887. break;
  888. }
  889. }
  890. if (npass >= NPASSES_MIN) {
  891. /* printf("pass %s sampler\n", "geometric"); */
  892. return true;
  893. } else {
  894. printf("fail %s sampler\n", "geometric");
  895. return false;
  896. }
  897. }
  898. /**
  899. * Divide the support of <b>dist</b> into histogram bins in <b>logP</b>. Start
  900. * at the 1st percentile and ending at the 99th percentile. Pick the bin
  901. * boundaries using linear interpolation so that they are uniformly spaced.
  902. *
  903. * In each bin logP[i] we insert the expected log-probability that a sampled
  904. * value will fall into that bin. We will use this as the null hypothesis of
  905. * the psi test.
  906. *
  907. * Set logP[i] = log(CDF(x_i) - CDF(x_{i-1})), where x_-1 = -inf, x_n =
  908. * +inf, and x_i = i*(hi - lo)/(n - 2).
  909. */
  910. static void
  911. bin_cdfs(const struct dist *dist, double lo, double hi, double *logP, size_t n)
  912. {
  913. #define CDF(x) dist_cdf(dist, x)
  914. #define SF(x) dist_sf(dist, x)
  915. const double w = (hi - lo)/(n - 2);
  916. double halfway = dist_icdf(dist, 0.5);
  917. double x_0, x_1;
  918. size_t i;
  919. size_t n2 = ceil_to_size_t((halfway - lo)/w);
  920. tor_assert(lo <= halfway);
  921. tor_assert(halfway <= hi);
  922. tor_assert(n2 <= n);
  923. x_1 = lo;
  924. logP[0] = log(CDF(x_1) - 0); /* 0 = CDF(-inf) */
  925. for (i = 1; i < n2; i++) {
  926. x_0 = x_1;
  927. /* do the linear interpolation */
  928. x_1 = (i <= n/2 ? lo + i*w : hi - (n - 2 - i)*w);
  929. /* set the expected log-probability */
  930. logP[i] = log(CDF(x_1) - CDF(x_0));
  931. }
  932. x_0 = hi;
  933. logP[n - 1] = log(SF(x_0) - 0); /* 0 = SF(+inf) = 1 - CDF(+inf) */
  934. /* In this loop we are filling out the high part of the array. We are using
  935. * SF because in these cases the CDF is near 1 where precision is lower. So
  936. * instead we are using SF near 0 where the precision is higher. We have
  937. * SF(t) = 1 - CDF(t). */
  938. for (i = 1; i < n - n2; i++) {
  939. x_1 = x_0;
  940. /* do the linear interpolation */
  941. x_0 = (i <= n/2 ? hi - i*w : lo + (n - 2 - i)*w);
  942. /* set the expected log-probability */
  943. logP[n - i - 1] = log(SF(x_0) - SF(x_1));
  944. }
  945. #undef SF
  946. #undef CDF
  947. }
  948. /**
  949. * Draw NSAMPLES samples from dist, counting the number of samples x in
  950. * the ith bin C[i] if x_{i-1} <= x < x_i, where x_-1 = -inf, x_n =
  951. * +inf, and x_i = i*(hi - lo)/(n - 2).
  952. */
  953. static void
  954. bin_samples(const struct dist *dist, double lo, double hi, size_t *C, size_t n)
  955. {
  956. const double w = (hi - lo)/(n - 2);
  957. size_t i;
  958. for (i = 0; i < NSAMPLES; i++) {
  959. double x = dist_sample(dist);
  960. size_t bin;
  961. if (x < lo)
  962. bin = 0;
  963. else if (x < hi)
  964. bin = 1 + floor_to_size_t((x - lo)/w);
  965. else
  966. bin = n - 1;
  967. tor_assert(bin < n);
  968. C[bin]++;
  969. }
  970. }
  971. /**
  972. * Carry out a Psi test on <b>dist</b>.
  973. *
  974. * Sample NSAMPLES from dist, putting them in bins from -inf to lo to
  975. * hi to +inf, and apply up to two psi tests. True if at least one psi
  976. * test passes; false if not. False positive rate should be bounded by
  977. * 0.01^2 = 0.0001.
  978. */
  979. static bool
  980. test_psi_dist_sample(const struct dist *dist)
  981. {
  982. double logP[PSI_DF] = {0};
  983. unsigned ntry = NTRIALS, npass = 0;
  984. double lo = dist_icdf(dist, 1/(double)(PSI_DF + 2));
  985. double hi = dist_isf(dist, 1/(double)(PSI_DF + 2));
  986. /* Create the null hypothesis in logP */
  987. bin_cdfs(dist, lo, hi, logP, PSI_DF);
  988. /* Now run the test */
  989. while (ntry --> 0) {
  990. size_t C[PSI_DF] = {0};
  991. bin_samples(dist, lo, hi, C, PSI_DF);
  992. if (psi_test(C, logP, NSAMPLES)) {
  993. if (++npass >= NPASSES_MIN)
  994. break;
  995. }
  996. }
  997. /* Did we fail or succeed? */
  998. if (npass >= NPASSES_MIN) {
  999. /* printf("pass %s sampler\n", dist_name(dist));*/
  1000. return true;
  1001. } else {
  1002. printf("fail %s sampler\n", dist_name(dist));
  1003. return false;
  1004. }
  1005. }
  1006. /* This is the seed of the deterministic randomness */
  1007. static uint8_t rng_seed[16];
  1008. static crypto_xof_t *rng_xof = NULL;
  1009. /** Initialize the seed of the deterministic randomness. */
  1010. static void
  1011. init_deterministic_rand(void)
  1012. {
  1013. crypto_rand((char*)rng_seed, sizeof(rng_seed));
  1014. crypto_xof_free(rng_xof);
  1015. rng_xof = crypto_xof_new();
  1016. crypto_xof_add_bytes(rng_xof, rng_seed, sizeof(rng_seed));
  1017. }
  1018. static void
  1019. teardown_deterministic_rand(void)
  1020. {
  1021. crypto_xof_free(rng_xof);
  1022. }
  1023. static void
  1024. dump_seed(void)
  1025. {
  1026. printf("\n"
  1027. "NOTE: This is a stochastic test, and we expect it to fail from\n"
  1028. "time to time, with some low probability. If you see it fail more\n"
  1029. "than one trial in 100, though, please tell us.\n\n"
  1030. "Seed: %s\n",
  1031. hex_str((const char*)rng_seed, sizeof(rng_seed)));
  1032. }
  1033. /** Produce deterministic randomness for the stochastic tests using the global
  1034. * deterministic_rand_counter seed
  1035. *
  1036. * This function produces deterministic data over multiple calls iff it's
  1037. * called in the same call order with the same 'n' parameter (which is the
  1038. * case for the psi test). If not, outputs will deviate. */
  1039. static void
  1040. crypto_rand_deterministic(char *out, size_t n)
  1041. {
  1042. /* Use a XOF to squeeze bytes out of that silly counter */
  1043. tor_assert(rng_xof);
  1044. crypto_xof_squeeze_bytes(rng_xof, (uint8_t*)out, n);
  1045. }
  1046. static void
  1047. test_stochastic_uniform(void *arg)
  1048. {
  1049. (void) arg;
  1050. const struct uniform uniform01 = {
  1051. .base = UNIFORM(uniform01),
  1052. .a = 0,
  1053. .b = 1,
  1054. };
  1055. const struct uniform uniform_pos = {
  1056. .base = UNIFORM(uniform_pos),
  1057. .a = 1.23,
  1058. .b = 4.56,
  1059. };
  1060. const struct uniform uniform_neg = {
  1061. .base = UNIFORM(uniform_neg),
  1062. .a = -10,
  1063. .b = -1,
  1064. };
  1065. const struct uniform uniform_cross = {
  1066. .base = UNIFORM(uniform_cross),
  1067. .a = -1.23,
  1068. .b = 4.56,
  1069. };
  1070. const struct uniform uniform_subnormal = {
  1071. .base = UNIFORM(uniform_subnormal),
  1072. .a = 4e-324,
  1073. .b = 4e-310,
  1074. };
  1075. const struct uniform uniform_subnormal_cross = {
  1076. .base = UNIFORM(uniform_subnormal_cross),
  1077. .a = -4e-324,
  1078. .b = 4e-310,
  1079. };
  1080. bool ok = true, tests_failed = true;
  1081. init_deterministic_rand();
  1082. MOCK(crypto_rand, crypto_rand_deterministic);
  1083. ok &= test_psi_dist_sample(&uniform01.base);
  1084. ok &= test_psi_dist_sample(&uniform_pos.base);
  1085. ok &= test_psi_dist_sample(&uniform_neg.base);
  1086. ok &= test_psi_dist_sample(&uniform_cross.base);
  1087. ok &= test_psi_dist_sample(&uniform_subnormal.base);
  1088. ok &= test_psi_dist_sample(&uniform_subnormal_cross.base);
  1089. tt_assert(ok);
  1090. tests_failed = false;
  1091. done:
  1092. if (tests_failed) {
  1093. dump_seed();
  1094. }
  1095. teardown_deterministic_rand();
  1096. UNMOCK(crypto_rand);
  1097. }
  1098. static bool
  1099. test_stochastic_logistic_impl(double mu, double sigma)
  1100. {
  1101. const struct logistic dist = {
  1102. .base = LOGISTIC(dist),
  1103. .mu = mu,
  1104. .sigma = sigma,
  1105. };
  1106. /* XXX Consider some fancier logistic test. */
  1107. return test_psi_dist_sample(&dist.base);
  1108. }
  1109. static bool
  1110. test_stochastic_log_logistic_impl(double alpha, double beta)
  1111. {
  1112. const struct log_logistic dist = {
  1113. .base = LOG_LOGISTIC(dist),
  1114. .alpha = alpha,
  1115. .beta = beta,
  1116. };
  1117. /* XXX Consider some fancier log logistic test. */
  1118. return test_psi_dist_sample(&dist.base);
  1119. }
  1120. static bool
  1121. test_stochastic_weibull_impl(double lambda, double k)
  1122. {
  1123. const struct weibull dist = {
  1124. .base = WEIBULL(dist),
  1125. .lambda = lambda,
  1126. .k = k,
  1127. };
  1128. /*
  1129. * XXX Consider applying a Tiku-Singh test:
  1130. *
  1131. * M.L. Tiku and M. Singh, `Testing the two-parameter
  1132. * Weibull distribution', Communications in Statistics --
  1133. * Theory and Methods A10(9), 1981, 907--918.
  1134. *https://www.tandfonline.com/doi/pdf/10.1080/03610928108828082?needAccess=true
  1135. */
  1136. return test_psi_dist_sample(&dist.base);
  1137. }
  1138. static bool
  1139. test_stochastic_genpareto_impl(double mu, double sigma, double xi)
  1140. {
  1141. const struct genpareto dist = {
  1142. .base = GENPARETO(dist),
  1143. .mu = mu,
  1144. .sigma = sigma,
  1145. .xi = xi,
  1146. };
  1147. /* XXX Consider some fancier GPD test. */
  1148. return test_psi_dist_sample(&dist.base);
  1149. }
  1150. static void
  1151. test_stochastic_genpareto(void *arg)
  1152. {
  1153. bool ok = 0;
  1154. bool tests_failed = true;
  1155. (void) arg;
  1156. init_deterministic_rand();
  1157. MOCK(crypto_rand, crypto_rand_deterministic);
  1158. ok = test_stochastic_genpareto_impl(0, 1, -0.25);
  1159. tt_assert(ok);
  1160. ok = test_stochastic_genpareto_impl(0, 1, -1e-30);
  1161. tt_assert(ok);
  1162. ok = test_stochastic_genpareto_impl(0, 1, 0);
  1163. tt_assert(ok);
  1164. ok = test_stochastic_genpareto_impl(0, 1, 1e-30);
  1165. tt_assert(ok);
  1166. ok = test_stochastic_genpareto_impl(0, 1, 0.25);
  1167. tt_assert(ok);
  1168. ok = test_stochastic_genpareto_impl(-1, 1, -0.25);
  1169. tt_assert(ok);
  1170. ok = test_stochastic_genpareto_impl(1, 2, 0.25);
  1171. tt_assert(ok);
  1172. tests_failed = false;
  1173. done:
  1174. if (tests_failed) {
  1175. dump_seed();
  1176. }
  1177. teardown_deterministic_rand();
  1178. UNMOCK(crypto_rand);
  1179. }
  1180. static void
  1181. test_stochastic_geometric(void *arg)
  1182. {
  1183. bool ok = 0;
  1184. bool tests_failed = true;
  1185. (void) arg;
  1186. init_deterministic_rand();
  1187. MOCK(crypto_rand, crypto_rand_deterministic);
  1188. ok = test_stochastic_geometric_impl(0.1);
  1189. tt_assert(ok);
  1190. ok = test_stochastic_geometric_impl(0.5);
  1191. tt_assert(ok);
  1192. ok = test_stochastic_geometric_impl(0.9);
  1193. tt_assert(ok);
  1194. ok = test_stochastic_geometric_impl(1);
  1195. tt_assert(ok);
  1196. tests_failed = false;
  1197. done:
  1198. if (tests_failed) {
  1199. dump_seed();
  1200. }
  1201. teardown_deterministic_rand();
  1202. UNMOCK(crypto_rand);
  1203. }
  1204. static void
  1205. test_stochastic_logistic(void *arg)
  1206. {
  1207. bool ok = 0;
  1208. bool tests_failed = true;
  1209. (void) arg;
  1210. init_deterministic_rand();
  1211. MOCK(crypto_rand, crypto_rand_deterministic);
  1212. ok = test_stochastic_logistic_impl(0, 1);
  1213. tt_assert(ok);
  1214. ok = test_stochastic_logistic_impl(0, 1e-16);
  1215. tt_assert(ok);
  1216. ok = test_stochastic_logistic_impl(1, 10);
  1217. tt_assert(ok);
  1218. ok = test_stochastic_logistic_impl(-10, 100);
  1219. tt_assert(ok);
  1220. tests_failed = false;
  1221. done:
  1222. if (tests_failed) {
  1223. dump_seed();
  1224. }
  1225. teardown_deterministic_rand();
  1226. UNMOCK(crypto_rand);
  1227. }
  1228. static void
  1229. test_stochastic_log_logistic(void *arg)
  1230. {
  1231. bool ok = 0;
  1232. bool tests_failed = true;
  1233. (void) arg;
  1234. init_deterministic_rand();
  1235. MOCK(crypto_rand, crypto_rand_deterministic);
  1236. ok = test_stochastic_log_logistic_impl(1, 1);
  1237. tt_assert(ok);
  1238. ok = test_stochastic_log_logistic_impl(1, 10);
  1239. tt_assert(ok);
  1240. ok = test_stochastic_log_logistic_impl(M_E, 1e-1);
  1241. tt_assert(ok);
  1242. ok = test_stochastic_log_logistic_impl(exp(-10), 1e-2);
  1243. tt_assert(ok);
  1244. tests_failed = false;
  1245. done:
  1246. if (tests_failed) {
  1247. dump_seed();
  1248. }
  1249. teardown_deterministic_rand();
  1250. UNMOCK(crypto_rand);
  1251. }
  1252. static void
  1253. test_stochastic_weibull(void *arg)
  1254. {
  1255. bool ok = 0;
  1256. bool tests_failed = true;
  1257. (void) arg;
  1258. init_deterministic_rand();
  1259. MOCK(crypto_rand, crypto_rand_deterministic);
  1260. ok = test_stochastic_weibull_impl(1, 0.5);
  1261. tt_assert(ok);
  1262. ok = test_stochastic_weibull_impl(1, 1);
  1263. tt_assert(ok);
  1264. ok = test_stochastic_weibull_impl(1, 1.5);
  1265. tt_assert(ok);
  1266. ok = test_stochastic_weibull_impl(1, 2);
  1267. tt_assert(ok);
  1268. ok = test_stochastic_weibull_impl(10, 1);
  1269. tt_assert(ok);
  1270. tests_failed = false;
  1271. done:
  1272. if (tests_failed) {
  1273. dump_seed();
  1274. }
  1275. teardown_deterministic_rand();
  1276. UNMOCK(crypto_rand);
  1277. }
  1278. struct testcase_t prob_distr_tests[] = {
  1279. { "logit_logistics", test_logit_logistic, TT_FORK, NULL, NULL },
  1280. { "log_logistic", test_log_logistic, TT_FORK, NULL, NULL },
  1281. { "weibull", test_weibull, TT_FORK, NULL, NULL },
  1282. { "genpareto", test_genpareto, TT_FORK, NULL, NULL },
  1283. { "uniform_interval", test_uniform_interval, TT_FORK, NULL, NULL },
  1284. END_OF_TESTCASES
  1285. };
  1286. struct testcase_t slow_stochastic_prob_distr_tests[] = {
  1287. { "stochastic_genpareto", test_stochastic_genpareto, TT_FORK, NULL, NULL },
  1288. { "stochastic_geometric", test_stochastic_geometric, TT_FORK, NULL, NULL },
  1289. { "stochastic_uniform", test_stochastic_uniform, TT_FORK, NULL, NULL },
  1290. { "stochastic_logistic", test_stochastic_logistic, TT_FORK, NULL, NULL },
  1291. { "stochastic_log_logistic", test_stochastic_log_logistic, TT_FORK, NULL,
  1292. NULL },
  1293. { "stochastic_weibull", test_stochastic_weibull, TT_FORK, NULL, NULL },
  1294. END_OF_TESTCASES
  1295. };