Tests: Derive check_histogram's tolerance from the sampling distribution

libeigen/eigen!2981

Closes #3130

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

test/rand.cpp

diff --git a/test/rand.cpp b/test/rand.cpp
index 99df7b9..4ed1a88 100644
--- a/test/rand.cpp
+++ b/test/rand.cpp
@@ -131,6 +131,24 @@
   RangeType bin_width_;
 };
 
+// A bin of probability p receives Binomial(count, p) samples, so the normalized hist(i) has
+// standard deviation sqrt((1 - p) / (count * p)) about 1. Bonferroni over the bins: k solving
+// erfc(k / sqrt(2)) = kFalseFailureRate / bins caps at kFalseFailureRate the probability that
+// any bin strays past k sigma.
+inline double histogram_k_sigma(int bins) {
+  constexpr double kFalseFailureRate = 1e-9;
+  const double target = kFalseFailureRate / double(bins);
+  double lo = 0.0, hi = 40.0;
+  for (int iter = 0; iter < 100; ++iter) {
+    const double mid = 0.5 * (lo + hi);
+    if (std::erfc(mid / std::sqrt(2.0)) > target)
+      lo = mid;
+    else
+      hi = mid;
+  }
+  return hi;
+}
+
 template <typename Scalar>
 void check_histogram(Scalar x, Scalar y, int bins) {
   constexpr int repeats = 10000;
@@ -143,12 +161,16 @@
       int bin = hist_helper.bin(r);
       hist(bin)++;
     }
-  //  Normalize bins by probability.
+  //  Normalize bins by probability, and bound each by its own sampling standard deviation.
   hist /= count;
+  const double k = histogram_k_sigma(bins);
+  Eigen::VectorXd tol(bins);
   for (int i = 0; i < bins; ++i) {
-    hist(i) = hist(i) / hist_helper.uniform_bin_probability(i);
+    const double p = hist_helper.uniform_bin_probability(i);
+    hist(i) = hist(i) / p;
+    tol(i) = k * std::sqrt((1.0 - p) / (count * p));
   }
-  VERIFY(((hist.array() - 1.0).abs() < 0.05).all());
+  VERIFY(((hist.array() - 1.0).abs() < tol.array()).all());
 }
 
 template <typename Scalar>
@@ -163,12 +185,16 @@
       int bin = hist_helper.bin(r);
       hist(bin)++;
     }
-  //  Normalize bins by probability.
+  //  Normalize bins by probability, and bound each by its own sampling standard deviation.
   hist /= count;
+  const double k = histogram_k_sigma(bins);
+  Eigen::VectorXd tol(bins);
   for (int i = 0; i < bins; ++i) {
-    hist(i) = hist(i) / hist_helper.uniform_bin_probability(i);
+    const double p = hist_helper.uniform_bin_probability(i);
+    hist(i) = hist(i) / p;
+    tol(i) = k * std::sqrt((1.0 - p) / (count * p));
   }
-  VERIFY(((hist.array() - 1.0).abs() < 0.05).all());
+  VERIFY(((hist.array() - 1.0).abs() < tol.array()).all());
 }
 
 template <>