diff --git a/regress.c b/regress.c index d68526d..af69e11 100644 --- 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; iestimated_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, °rees_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;