31:T546, md` ## Details The robustness part of the loess algorithm aims to reduce the effect of "outliers" by downweighting observations with relatively large residuals. The downweighting is gradual up to 6x the median absolute residual, when the weighting multiplier becomes 0. The issue is when the curve fits so well that the median residual is practically zero (because over half the smoothed values are very close to the data values), making many other residuals beyond the 6x threshold. Then those areas get ignored by subsequent passes since they have zero weight. This flaw isn't specific to D3. It occurs in any implementation that includes Cleveland's robustness iterations, which includes R's \`\`loess\`\` and \`\`lowess\`\` packages. In creating this test, I discovered an unrelated bug in the D3's [\`\`regressionLoess\`\`](https://github.com/HarryStevens/d3-regression) implementation, where it used the sorted residuals instead of the actual residuals for the robustness adjustments. I use a cloned and corrected version here (which also lets me tune the robustness iterations) and have submitted a pull request to fix the D3 version (now integrated). [This notebook was originally a fork of [LOESS Regression by @harrystevens](https://observablehq.com/@harrystevens/loess-regression), but the linkage was lost among interim versions.] `32:Tcb0,// Adapted from science.js by Jason Davies // License: https://github.com/jasondavies/science.js/blob/master/LICENSE // Source: https://github.com/jasondavies/science.js/blob/master/src/stats/loess.js // Adapted from vega-statistics by Jeffrey Heer // License: https://github.com/vega/vega/blob/f058b099decad9db78301405dd0d2e9d8ba3d51a/LICENSE // Source: https://github.com/vega/vega/blob/f21cb8792b4e0cbe2b1a3fd44b0f5db370dbaadb/packages/vega-statistics/src/regression/loess.js //import { ols } from "./utils/ols"; //import { points } from "./utils/points"; tuneableLoess = function() { let x = d => d[0], y = d => d[1], bandwidth = .3, maxiters = 2, epsilon = 1e-12, minweight = 1e-12; function loess(data, reweightFn = 1) { // 1=median*6, 2 = max(median*6, q90 * 2) const [xv, yv, ux, uy] = points(data, x, y, true), n = xv.length, bw = Math.max(2, ~~(bandwidth * n)), // # nearest neighbors yhat = new Float64Array(n), residuals = new Float64Array(n), robustWeights = new Float64Array(n).fill(1); for (let iter = -1; ++iter <= maxiters; ) { const interval = [0, bw - 1]; for (let i = 0; i < n; ++i) { const dx = xv[i], i0 = interval[0], i1 = interval[1], edge = (dx - xv[i0]) > (xv[i1] - dx) ? i0 : i1; let W = 0, X = 0, Y = 0, XY = 0, X2 = 0, denom = 1 / Math.abs(xv[edge] - dx || 1); // Avoid singularity for (let k = i0; k <= i1; ++k) { const xk = xv[k], yk = yv[k], w = tricube(Math.abs(dx - xk) * denom) * robustWeights[k], xkw = xk * w; W += w; X += xkw; Y += yk * w; XY += yk * xkw; X2 += xk * xkw; } // Linear regression fit const [a, b] = ols(X / W, Y / W, XY / W, X2 / W); yhat[i] = a + b * dx; residuals[i] = Math.abs(yv[i] - yhat[i]); updateInterval(xv, i + 1, interval); } if (iter === maxiters) { break; } const medianResidual = d3.median(residuals); const q90Residual = d3.quantile(residuals, 0.9); // q90(residuals); const cutoff = reweightFn == 1 ? 6 * medianResidual : Math.max(6 * medianResidual, 2 * q90Residual) // if (Math.abs(cutoff) < epsilon) // break; for (let i = 0, arg, w; i < n; ++i){ arg = residuals[i] / cutoff; // Default to epsilon (rather than zero) for large deviations // Keeping weights tiny but non-zero prevents singularites robustWeights[i] = (arg >= 1) ? minweight : ((w = 1 - arg * arg) * w); } } return output(xv, yhat, ux, uy); } loess.bandwidth = function(bw) { return arguments.length ? (bandwidth = bw, loess) : bandwidth; }; loess.minweight = function(e) { return arguments.length ? (minweight = e, loess) : minweight; }; loess.maxiters = function(mi) { return arguments.length ? (maxiters = mi, loess) : maxiters; }; loess.x = function(fn) { return arguments.length ? (x = fn, loess) : x; }; loess.y = function(fn) { return arguments.length ? (y = fn, loess) : y; }; return loess; } 33:[] 2f:{"title":"Loess Regression Robustness Flaw","cells":[{"id":557,"value":"md`# Loess Regression Robustness Flaw\n\nThis notebook demonstrates a minor flaw in robust loess regression which\noccurs when the original data is already very smooth. In that case,\nthe usual robustness passes can create artifacts in the data. Choose a data set.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":638,"value":"viewof datasetname = select({\n// title: \"Data\",\n options: [\"Cleveland Curve\",\"Mächler Spike\",\"Nice Outlier\"],\n value: \"Cleveland Curve\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":661,"value":"html`${datasetname == \"Cleveland Curve\" ? \"Artificial example from Cleveland's Visualizing Data book (1993), constructed by stitching together three smooth functions.\"\n: datasetname == \"Nice Outlier\" ?\n \"A rising sine curve with a couple obvious outliers. Any robustness treatment should eliminate the outlier affect on the smooth curve.\"\n : \"Breakdown case noted by Martin Mächler, but scaled down to fit on our axes. Use a noise scale value of 0.03 or less to get the full effect of an extreme spike. The new rule does better for moderate bandwidth values, but both rules fail at very small bandwidths.\"}`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"viewof bandwidth = slider({\n title: \"Local fit bandwidth\",\n min: .02, \n max: 1, \n step: .01,\n value: .25\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":484,"value":"viewof iters = slider({\n title: \"Robustness iterations\",\n min: 0, \n max: 5, \n step: 1,\n value: 2\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":397,"value":"viewof noise = slider({\n title: \"Noise scale\",\n min: 0, \n max: 1, \n step: .01,\n value: 0.5\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":576,"value":"viewof reweightChoice = radio({\n title: 'Robust reweighting limit',\n //description: 'for robustness iterations',\n options: [\n { label: 'Cleveland: 6 × q₅₀', value: 'q50' },\n { label: 'Adjusted: max(6 × q₅₀, 2 × q₉₀)', value: 'q90' },\n { label: 'Both', value: 'both' },\n ],\n value: 'q50'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":148,"value":"{\n const svg = d3.select(DOM.svg(innerWidth + margin.left + margin.right, innerHeight + margin.top + margin.bottom))\n \n const g = svg.append(\"g\")\n .attr(\"transform\", `translate(${margin.left}, ${margin.top})`);\n\n g.append(\"g\")\n .attr(\"class\", \"axis\")\n .call(xAxisGenerator);\n\n g.append(\"g\")\n .attr(\"class\", \"axis\")\n .attr(\"transform\", `translate(${innerWidth})`)\n .call(yAxisGenerator);\n\n g.selectAll(\"circle\")\n .data(data)\n .enter().append(\"circle\")\n .attr(\"r\", 2)\n .attr(\"cx\", d => xScale(d.x))\n .attr(\"cy\", d => yScale(d.y));\n\n g.append(\"path\")\n .attr(\"class\", \"raw\")\n .attr(\"d\", lineGenerator2(dataRaw));\n\n if (reweightChoice != \"q90\")\n g.append(\"path\")\n .attr(\"class\", \"regression1\")\n .datum(regressionGenerator(data, 1))\n .attr(\"d\", lineGenerator);\n \n if (reweightChoice != \"q50\")\n g.append(\"path\")\n .attr(\"class\", \"regression2\")\n .datum(regressionGenerator(data, 2))\n .attr(\"d\", lineGenerator);\n\n \n return svg.node()\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":567,"value":"md`## Graphics helper funcitons`","pinned":false,"hidden":false,"mode":"ojs"},{"id":34,"value":"regressionGenerator = tuneableLoess()\n .x(d => d.x)\n .y(d => d.y)\n .bandwidth(bandwidth)\n .maxiters(iters)\n //.reweightFn(reweightChoice == \"q50\" ? 1 : 2)\n //.minweight(0)\n;","pinned":false,"hidden":false,"mode":"ojs"},{"id":356,"value":"lineGenerator = d3.line()\n .x(d => xScale(d[0]))\n .y(d => yScale(d[1]));","pinned":false,"hidden":false,"mode":"ojs"},{"id":449,"value":"lineGenerator2 = d3.line()\n .x(d => xScale(d.x))\n .y(d => yScale(d.y));","pinned":false,"hidden":false,"mode":"ojs"},{"id":114,"value":"md`## Data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":536,"value":"md`\\`\\`dataRaw\\`\\` contains the underlying funtion and \\`\\`dataNoise\\`\\` has the unscaled random normal noise. The noise scaling is not linear, so that extremely small noise can be observed (needed for Maechler case).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":400,"value":"data = dataRaw.map(o => { return {\n x: o.x,\n y0:o.y + dataNoise[o.x-1] * (Math.exp(noise) - 1),\n y1:o.y + dataNoise[o.x-1] * Math.pow(noise, 2) * (noise > 0 ? Math.min(1.0, 0.1 / noise) : 1.0),\n y: o.y + dataNoise[o.x-1] * (noise > 0 ? Math.pow(noise, Math.max(2, Math.min(10, 0.3 / noise))) : 0.0)\n};})","pinned":false,"hidden":false,"mode":"ojs"},{"id":436,"value":"dataNoise = Array.from({length: 150}, d3.randomNormal())","pinned":false,"hidden":false,"mode":"ojs"},{"id":382,"value":"dataRaw = datasetname == \"Cleveland Curve\" ? dataClevelandCurve : datasetname == \"Nice Outlier\" ? dataNiceOutlier : dataSpike;\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":596,"value":"dataClevelandCurve = Array.from({length: 150}, (_, i) => {\n const n = 50;\n const z = (i % 50) / 50; // [0..1)\n const k = i + i;\n return {\n x: i + 1,\n z: z,\n y: 5 + (i < n ? -Math.sqrt(1 - z)\n : i < 2 * n ? Math.sin( z * 2 * Math.PI )\n : z * z / 2)\n };\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":645,"value":"dataNiceOutlier = Array.from({length: 150}, (_, i) => {\n const n = 50;\n const z = (i % 50) / 50; // [0..1)\n const k = i + i;\n return {\n x: i + 1,\n z: z,\n y: 4 + i / 150 - Math.sin(i / 20) / 3 + (i == 20 || i == 100 ? 1 : 0)\n };\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":646,"value":"dataSpike\n = Array.from({length: 150}, (_, i) => {\n const n = 50;\n const z = (i % 50) / 50; // [0..1)\n const k = i + i;\n return {\n x: i + 1,\n y: i == 75 ? 6.4 : 3.7\n };\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":150,"value":"md`## Dimensions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":152,"value":"margin = ({left: 24, right: 11, top: 14, bottom: 21})","pinned":false,"hidden":false,"mode":"ojs"},{"id":153,"value":"base = Math.min(width, 600)","pinned":false,"hidden":false,"mode":"ojs"},{"id":156,"value":"innerWidth = base - margin.left - margin.right","pinned":false,"hidden":false,"mode":"ojs"},{"id":158,"value":"innerHeight = Math.floor((base - margin.top - margin.bottom) * 0.8)","pinned":false,"hidden":false,"mode":"ojs"},{"id":136,"value":"md`## Scales\nThese are hard-coded based on the Cleveland Curve data, so the other data sets have to be scaled to fit.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"xScale = d3.scaleLinear()\n .domain([0, 150])\n .range([0, innerWidth])","pinned":false,"hidden":false,"mode":"ojs"},{"id":145,"value":"yScale = d3.scaleLinear()\n .domain([3.5, 6.5])\n .range([innerHeight, 0])","pinned":false,"hidden":false,"mode":"ojs"},{"id":130,"value":"md`## Axes`","pinned":false,"hidden":false,"mode":"ojs"},{"id":131,"value":"xAxisGenerator = d3.axisBottom(xScale)\n .tickSize(innerHeight + 10)\n .tickValues(d3.range(0, 150, 25))\n //.tickFormat(d => d % 100 === 0 ? d : \"`\" + d.toString().slice(2, 4))","pinned":false,"hidden":false,"mode":"ojs"},{"id":133,"value":"yAxisGenerator = d3.axisLeft(yScale)\n .tickSize(innerWidth + 10)\n .tickValues(d3.range(3, 7))\n .tickFormat(d => d)","pinned":false,"hidden":false,"mode":"ojs"},{"id":119,"value":"md`## Styles`","pinned":false,"hidden":false,"mode":"ojs"},{"id":127,"value":"html`\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":515,"value":"md`## Loess regression code`","pinned":false,"hidden":false,"mode":"ojs"},{"id":513,"value":"md`\nThe functions below essentially clone [d3-regression](https://github.com/HarryStevens/d3-regression) for loess with these modifications:\n* Fix bug where \\`\\`median\\`\\` function would modify the input array (affecting robust residuals). [Nov2025 update: now merged into main package] Use \\`\\`d3.median()\\`\\` instead.\n* Make the robust iteration count tuneable.\n* Make the weight for extreme outliers tuneable.\n* Add support for 90th quantile threshold for reweighting residuals.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":477,"value":"// Ordinary Least Squares from vega-statistics by Jeffrey Heer\n// License: https://github.com/vega/vega/blob/f058b099decad9db78301405dd0d2e9d8ba3d51a/LICENSE\n// Source: https://github.com/vega/vega/blob/f058b099decad9db78301405dd0d2e9d8ba3d51a/packages/vega-statistics/src/regression/ols.js\nfunction ols(uX, uY, uXY, uX2) {\n const delta = uX2 - uX * uX,\n slope = Math.abs(delta) < 1e-24 ? 0 : (uXY - uX * uY) / delta,\n intercept = uY - slope * uX;\n\n return [intercept, slope];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":474,"value":"function points(data, x, y, sort) {\n data = data.filter(d => {\n let u = x(d), v = y(d);\n return u != null && isFinite(u) && v != null && isFinite(v);\n });\n\n if (sort) {\n data.sort((a, b) => x(a) - x(b));\n }\n\n const n = data.length,\n X = new Float64Array(n),\n Y = new Float64Array(n);\n\n // extract values, calculate means\n let ux = 0, uy = 0, xv, yv, d;\n for (let i = 0; i < n; ) {\n d = data[i];\n X[i] = xv = +x(d);\n Y[i] = yv = +y(d);\n ++i;\n ux += (xv - ux) / i;\n uy += (yv - uy) / i;\n }\n\n // mean center the data\n for (let i = 0; i < n; ++i) {\n X[i] -= ux;\n Y[i] -= uy;\n }\n\n return [X, Y, ux, uy];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":458,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":462,"value":"\n// Weighting kernel for local regression\nfunction tricube(x) {\n return (x = 1 - x * x * x) * x * x;\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":464,"value":"\n// Advance sliding window interval of nearest neighbors\nfunction updateInterval(xv, i, interval) {\n let val = xv[i],\n left = interval[0],\n right = interval[1] + 1;\n\n if (right >= xv.length) return;\n\n // Step right if distance to new right edge is <= distance to old left edge\n // Step when distance is equal to ensure movement over duplicate x values\n while (i > left && (xv[right] - val) <= (val - xv[left])) {\n interval[0] = ++left;\n interval[1] = right;\n ++right;\n }\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":461,"value":"\n// Generate smoothed output points\n// Average points with repeated x values\nfunction output(xv, yhat, ux, uy) {\n const n = xv.length, out = [];\n let i = 0, cnt = 0, prev = [], v;\n\n for (; i