in TestApplications/Web/js/lib/components/d3-3.4.13/lib/science/science.stats.js [460:584]
function smooth(xval, yval, weights) {
var n = xval.length,
i;
if (n !== yval.length) throw {error: "Mismatched array lengths"};
if (n == 0) throw {error: "At least one point required."};
if (arguments.length < 3) {
weights = [];
i = -1; while (++i < n) weights[i] = 1;
}
science_stats_loessFiniteReal(xval);
science_stats_loessFiniteReal(yval);
science_stats_loessFiniteReal(weights);
science_stats_loessStrictlyIncreasing(xval);
if (n == 1) return [yval[0]];
if (n == 2) return [yval[0], yval[1]];
var bandwidthInPoints = Math.floor(bandwidth * n);
if (bandwidthInPoints < 2) throw {error: "Bandwidth too small."};
var res = [],
residuals = [],
robustnessWeights = [];
// Do an initial fit and 'robustnessIters' robustness iterations.
// This is equivalent to doing 'robustnessIters+1' robustness iterations
// starting with all robustness weights set to 1.
i = -1; while (++i < n) {
res[i] = 0;
residuals[i] = 0;
robustnessWeights[i] = 1;
}
var iter = -1;
while (++iter <= robustnessIters) {
var bandwidthInterval = [0, bandwidthInPoints - 1];
// At each x, compute a local weighted linear regression
var x;
i = -1; while (++i < n) {
x = xval[i];
// Find out the interval of source points on which
// a regression is to be made.
if (i > 0) {
science_stats_loessUpdateBandwidthInterval(xval, weights, i, bandwidthInterval);
}
var ileft = bandwidthInterval[0],
iright = bandwidthInterval[1];
// Compute the point of the bandwidth interval that is
// farthest from x
var edge = (xval[i] - xval[ileft]) > (xval[iright] - xval[i]) ? ileft : iright;
// Compute a least-squares linear fit weighted by
// the product of robustness weights and the tricube
// weight function.
// See http://en.wikipedia.org/wiki/Linear_regression
// (section "Univariate linear case")
// and http://en.wikipedia.org/wiki/Weighted_least_squares
// (section "Weighted least squares")
var sumWeights = 0,
sumX = 0,
sumXSquared = 0,
sumY = 0,
sumXY = 0,
denom = Math.abs(1 / (xval[edge] - x));
for (var k = ileft; k <= iright; ++k) {
var xk = xval[k],
yk = yval[k],
dist = k < i ? x - xk : xk - x,
w = science_stats_loessTricube(dist * denom) * robustnessWeights[k] * weights[k],
xkw = xk * w;
sumWeights += w;
sumX += xkw;
sumXSquared += xk * xkw;
sumY += yk * w;
sumXY += yk * xkw;
}
var meanX = sumX / sumWeights,
meanY = sumY / sumWeights,
meanXY = sumXY / sumWeights,
meanXSquared = sumXSquared / sumWeights;
var beta = (Math.sqrt(Math.abs(meanXSquared - meanX * meanX)) < accuracy)
? 0 : ((meanXY - meanX * meanY) / (meanXSquared - meanX * meanX));
var alpha = meanY - beta * meanX;
res[i] = beta * x + alpha;
residuals[i] = Math.abs(yval[i] - res[i]);
}
// No need to recompute the robustness weights at the last
// iteration, they won't be needed anymore
if (iter === robustnessIters) {
break;
}
// Recompute the robustness weights.
// Find the median residual.
var sortedResiduals = residuals.slice();
sortedResiduals.sort();
var medianResidual = sortedResiduals[Math.floor(n / 2)];
if (Math.abs(medianResidual) < accuracy)
break;
var arg,
w;
i = -1; while (++i < n) {
arg = residuals[i] / (6 * medianResidual);
robustnessWeights[i] = (arg >= 1) ? 0 : ((w = 1 - arg * arg) * w);
}
}
return res;
}