| 119854 | // Source: https://github.com/jasondavies/science.js/blob/master/src/stats/loess.js |
| 119855 | // License: https://github.com/jasondavies/science.js/blob/master/LICENSE |
| 119856 | function loess(data, x, y, bandwidth) { |
| 119857 | const [xv, yv, ux, uy] = points(data, x, y, true), n = xv.length, bw = Math.max(2, ~~(bandwidth * n)), // # nearest neighbors |
| 119858 | yhat = new Float64Array(n), residuals = new Float64Array(n), robustWeights = new Float64Array(n).fill(1); |
| 119859 | for(let iter = -1; ++iter <= maxiters;){ |
| 119860 | const interval = [ |
| 119861 | 0, |
| 119862 | bw - 1 |
| 119863 | ]; |
| 119864 | for(let i = 0; i < n; ++i){ |
| 119865 | const dx = xv[i], i0 = interval[0], i1 = interval[1], edge = dx - xv[i0] > xv[i1] - dx ? i0 : i1; |
| 119866 | let W = 0, X = 0, Y = 0, XY = 0, X2 = 0; |
| 119867 | const denom = 1 / Math.abs(xv[edge] - dx || 1); // avoid singularity! |
| 119868 | for(let k = i0; k <= i1; ++k){ |
| 119869 | const xk = xv[k], yk = yv[k], w = tricube(Math.abs(dx - xk) * denom) * robustWeights[k], xkw = xk * w; |
| 119870 | W += w; |
| 119871 | X += xkw; |
| 119872 | Y += yk * w; |
| 119873 | XY += yk * xkw; |
| 119874 | X2 += xk * xkw; |
| 119875 | } // linear regression fit |
| 119876 | const [a, b] = ols(X / W, Y / W, XY / W, X2 / W); |
| 119877 | yhat[i] = a + b * dx; |
| 119878 | residuals[i] = Math.abs(yv[i] - yhat[i]); |
| 119879 | updateInterval(xv, i + 1, interval); |
| 119880 | } |
| 119881 | if (iter === maxiters) break; |
| 119882 | const medianResidual = (0, _d3Array.median)(residuals); |
| 119883 | if (Math.abs(medianResidual) < epsilon) break; |
| 119884 | for(let i1 = 0, arg, w; i1 < n; ++i1){ |
| 119885 | arg = residuals[i1] / (6 * medianResidual); // default to epsilon (rather than zero) for large deviations |
| 119886 | // keeping weights tiny but non-zero prevents singularites |
| 119887 | robustWeights[i1] = arg >= 1 ? epsilon : (w = 1 - arg * arg) * w; |
| 119888 | } |
| 119889 | } |
| 119890 | return output(xv, yhat, ux, uy); |
| 119891 | } // weighting kernel for local regression |
| 119892 | function tricube(x) { |
| 119893 | return (x = 1 - x * x * x) * x * x; |
| 119894 | } // advance sliding window interval of nearest neighbors |