In weight calculation use unweighted variance from last regression

This fixes a positive feedback where weights could reach inf.
Also change the SD_TO_DIST_RATIO constant to get close to the original
response.

diff --git a/regress.c b/regress.c
index d68526d660a7baabc377754ab65b409717f34865..af69e113ee32e3a0e26d4df4a623f7b5119909c6 100644 (file)
--- a/regress.c
+++ b/regress.c
@@ -232,6 +232,7 @@ RGR_FindBestRegression
  double *b0,                    /* estimated y axis intercept */
  double *b1,                    /* estimated slope */
  double *s2,                    /* estimated variance of data points */
+ double *us2,                   /* estimated unweighted variance of data points */
  
  double *sb0,                   /* estimated standard deviation of
                                    intercept */
@@ -250,7 +251,7 @@ RGR_FindBestRegression
 {
   double P, Q, U, V, W; /* total */
   double resid[MAX_POINTS * REGRESS_RUNS_RATIO];
-  double ss;
+  double ss, uss;
   double a, b, u, ui, aa;
 
   int start, resid_start, nruns, npoints;
@@ -310,17 +311,20 @@ RGR_FindBestRegression
   *b1 = b;
   *b0 = a;
 
-  ss = 0.0;
+  ss = uss = 0.0;
   for (i=start; i<n; i++) {
     ss += resid[i - resid_start]*resid[i - resid_start] / w[i];
+    uss += resid[i - resid_start]*resid[i - resid_start];
   }
 
   npoints = n - start;
-  ss /= (double)(npoints - 2);
+  ss /= npoints - 2;
+  uss /= npoints - 2;
   *sb1 = sqrt(ss / V);
   aa = u * (*sb1);
   *sb0 = sqrt((ss / W) + (aa * aa));
-  *s2 = ss * (double) npoints / W;
+  *s2 = ss * npoints / W;
+  *us2 = uss;
 
   *new_start = start;
   *dof = npoints - 2;

diff --git a/regress.h b/regress.h
index 4c11ca578f5a04d9968f06ee1ecb6a95415276dc..1c618ae8120b50eae0593b09110115464d38b43d 100644 (file)
--- a/regress.h
+++ b/regress.h
@@ -86,6 +86,7 @@ RGR_FindBestRegression
  double *b0,                    /* estimated y axis intercept */
  double *b1,                    /* estimated slope */
  double *s2,                    /* estimated variance of data points */
+ double *us2,                   /* estimated unweighted variance of data points */
  
  double *sb0,                   /* estimated standard deviation of
                                    intercept */

diff --git a/sourcestats.c b/sourcestats.c
index 9560eac413eb90e6c42c23f10a291f9e74549f8e..efd295cf4a5d43fc9831bfc002bebbc659111b2a 100644 (file)
--- a/sourcestats.c
+++ b/sourcestats.c
@@ -105,6 +105,9 @@ struct SST_Stats_Record {
   /* This is the estimated residual variance of the data points */
   double variance;
 
+  /* This is the estimated unweighted variance of the data points */
+  double uvariance;
+
   /* This array contains the sample epochs, in terms of the local
      clock. */
   struct timeval sample_times[MAX_SAMPLES * REGRESS_RUNS_RATIO];
@@ -190,7 +193,7 @@ SST_CreateInstance(unsigned long refid, IPAddr *addr)
   inst->estimated_offset_sd = 86400.0; /* Assume it's at least within a day! */
   inst->offset_time.tv_sec = 0;
   inst->offset_time.tv_usec = 0;
-  inst->variance = 16.0;
+  inst->variance = inst->uvariance = 16.0;
   inst->nruns = 0;
   return inst;
 }
@@ -364,7 +367,7 @@ find_min_delay_sample(SST_Stats inst)
    time.  E.g. a value of 4 means that we think the standard deviation
    is four times the fluctuation  of the peer distance */
 
-#define SD_TO_DIST_RATIO 1.0
+#define SD_TO_DIST_RATIO 1.4
 
 /* ================================================== */
 /* This function runs the linear regression operation on the data.  It
@@ -382,7 +385,7 @@ SST_DoNewRegression(SST_Stats inst)
 
   int degrees_of_freedom;
   int best_start, times_back_start;
-  double est_intercept, est_slope, est_var, est_intercept_sd, est_slope_sd;
+  double est_intercept, est_slope, est_var, est_uvar, est_intercept_sd, est_slope_sd;
   int i, j, nruns;
   double min_distance;
   double sd_weight, sd;
@@ -405,7 +408,7 @@ SST_DoNewRegression(SST_Stats inst)
 
     /* And now, work out the weight vector */
 
-    sd = sqrt(inst->variance);
+    sd = sqrt(inst->uvariance);
     if (sd > min_distance || sd <= 0.0)
       sd = min_distance;
 
@@ -418,7 +421,7 @@ SST_DoNewRegression(SST_Stats inst)
   inst->regression_ok = RGR_FindBestRegression(times_back + inst->runs_samples,
                                          offsets + inst->runs_samples, weights,
                                          inst->n_samples, inst->runs_samples,
-                                         &est_intercept, &est_slope, &est_var,
+                                         &est_intercept, &est_slope, &est_var, &est_uvar,
                                          &est_intercept_sd, &est_slope_sd,
                                          &best_start, &nruns, &degrees_of_freedom);
 
@@ -433,6 +436,7 @@ SST_DoNewRegression(SST_Stats inst)
     inst->offset_time = inst->sample_times[inst->last_sample];
     inst->estimated_offset_sd = est_intercept_sd;
     inst->variance = est_var;
+    inst->uvariance = est_uvar;
     inst->nruns = nruns;
 
     stress = fabs(old_freq - inst->estimated_frequency) / old_skew;