quantiles: force step update with stable input values

The algorithm was designed for estimating quantiles in streams of
integer values. When the estimate is equal to the input value, the
step state variable does not change. This causes problems for the
floating-point adaptation used for measurents of delay in chrony.

One problem is numerical instability due to the strict comparison of
the input value and the current estimate.

Another problem is with signals that are so stable that the nanosecond
resolution of the system functions becomes the limitation. There is a
large difference in the value of the step state variable, which
determines how quickly the estimate will adapt to a new distribution,
between signals that are constant in the nanosecond resolution and
signals that can move in two nanoseconds.

Change the estimate update to never consider the input value equal to
the current estimate and don't set the estimate exactly to the input
value. Keep it off by a quarter of the minimum step to force jumping
around the input value if it's constant and decreasing the step variable
to negative values. Also fix the initial adjustment to step at least by
the minimum step (the original algorithm is described with ceil(), not
fabs()).

diff --git a/quantiles.c b/quantiles.c
index a579b29785e4f3af259bd913fad56c815b34515a..f4783d4d5ab9b2533821476fd704485ab071c05c 100644 (file)
--- a/quantiles.c
+++ b/quantiles.c
@@ -138,22 +138,26 @@ static void
 update_estimate(struct Quantile *quantile, double value, double p, double rand,
                 double min_step)
 {
-  if (value > quantile->est && rand > (1.0 - p)) {
+  if (value >= quantile->est) {
+    if (rand < (1.0 - p))
+        return;
     quantile->step += quantile->sign > 0 ? min_step : -min_step;
-    quantile->est += quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    quantile->est += quantile->step > min_step ? quantile->step : min_step;
     if (quantile->est > value) {
       quantile->step += value - quantile->est;
-      quantile->est = value;
+      quantile->est = value + min_step / 4.0;
     }
     if (quantile->sign < 0 && quantile->step > min_step)
       quantile->step = min_step;
     quantile->sign = 1;
-  } else if (value < quantile->est && rand > p) {
+  } else {
+    if (rand < p)
+      return;
     quantile->step += quantile->sign < 0 ? min_step : -min_step;
-    quantile->est -= quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    quantile->est -= quantile->step > min_step ? quantile->step : min_step;
     if (quantile->est < value) {
       quantile->step += quantile->est - value;
-      quantile->est = value;
+      quantile->est = value - min_step / 4.0;
     }
     if (quantile->sign > 0 && quantile->step > min_step)
       quantile->step = min_step;

diff --git a/test/unit/quantiles.c b/test/unit/quantiles.c
index 86907808e5a6824389df1f94625afa6eb19ff2ee..e08801e18aec9c82abf36468d6516bb3b4db2483 100644 (file)
--- a/test/unit/quantiles.c
+++ b/test/unit/quantiles.c
@@ -28,7 +28,7 @@ test_unit(void)
 {
   int i, j, k, min_k, max_k, q, r, in_order, out_order;
   QNT_Instance inst;
-  double x;
+  double x, x2;
 
   in_order = out_order = 0;
 
@@ -63,6 +63,32 @@ test_unit(void)
     QNT_Reset(inst);
     TEST_CHECK(inst->n_set == 0);
 
+    for (j = 0; j < 3; j++) {
+      QNT_Reset(inst);
+      x = (i + 950) / 1e9;
+      if (random() % 10)
+        x2 = (i + 951) / 1e9;
+      else
+        x2 = x;
+      for (k = 0; k < 1000; k++)
+        QNT_Accumulate(inst, random() % 2 ? x : x2);
+      for (k = 0; k < inst->n_quants; k++) {
+        TEST_CHECK(inst->quants[k].est > x - 0.4e-9);
+        TEST_CHECK(inst->quants[k].est < x2 + 0.4e-9);
+        TEST_CHECK(inst->quants[k].step < -15e-9);
+        TEST_CHECK(inst->quants[k].step > -1000e-9);
+        if (min_k * 2 == q && k < inst->repeat) {
+          if (x == x2) {
+            TEST_CHECK(inst->quants[k].step < -750e-9);
+            TEST_CHECK(inst->quants[k].step > -1000e-9);
+          } else {
+            TEST_CHECK(inst->quants[k].step < -350e-9);
+            TEST_CHECK(inst->quants[k].step > -600e-9);
+          }
+        }
+      }
+    }
+
     QNT_DestroyInstance(inst);
   }