Smoothers

The greybox.smoothers module provides two non-parametric smoothers that reproduce R’s stats::lowess() and stats::supsmu() to machine precision. Each wraps a native pybind11 extension and returns a dictionary with sorted abscissa and smoothed ordinate.

LOWESS

greybox.smoothers.lowess(x, y=None, f=0.6666666666666666, iter=3, delta=None)[source]

LOWESS smoother (Locally Weighted Scatterplot Smoothing).

Performs locally weighted polynomial regression using Cleveland’s LOWESS algorithm with a tricube weight function and iterative reweighting for robustness to outliers. The implementation matches R’s stats::lowess() to machine precision.

Parameters:

x (array_like) – The abscissa values. May also be a 2-D array with two columns, in which case the first column is used as x and the second as y.
y (array_like, optional) – The ordinate values. Required unless x already contains both columns.
f (float, default=2/3) – Smoother span – the fraction of points in the local neighbourhood used for each fit. Larger values produce smoother curves.
iter (int, default=3) – Number of robustifying iterations. Each iteration downweights observations with large residuals from the previous fit.
delta (float, optional) – Distance threshold for interpolation. Within delta of an evaluated point, the fit is linearly interpolated rather than recomputed. If None, defaults to 0.01 * (max(x) - min(x)).

Returns:

Dictionary with two keys:

"x" : numpy.ndarray of sorted abscissa values.
"y" : numpy.ndarray of smoothed ordinate values, aligned with "x".

Return type:

dict

Notes

The smoother fits a weighted linear regression at each point using a tricube weight function

\[w(u) = (1 - |u|^3)^3\]

applied to neighbours within the local span. Robustness iterations downweight outliers using a bisquare weight on the residuals.

References

Cleveland, W. S. (1979). “Robust Locally Weighted Regression and Smoothing Scatterplots”. Journal of the American Statistical Association, 74(368), 829-836. DOI: https://doi.org/10.1080/01621459.1979.10481038

Examples

Smooth a noisy sine wave:

>>> import numpy as np
>>> from greybox import lowess
>>> rng = np.random.default_rng(0)
>>> x = np.linspace(0, 6, 60)
>>> y = np.sin(x) + rng.normal(0, 0.2, 60)
>>> out = lowess(x, y, f=0.4)
>>> out["y"].shape
(60,)

Pass both columns as a 2-D array:

>>> xy = np.column_stack([x, y])
>>> out = lowess(xy)

Example:

import numpy as np
from greybox import lowess

rng = np.random.default_rng(0)
x = np.linspace(0, 6, 80)
y = np.sin(x) + rng.normal(0, 0.2, 80)

smoothed = lowess(x, y, f=0.4)
# smoothed["x"] are the sorted x values; smoothed["y"] is the smoothed curve.

Reference

Cleveland, W. S. (1979). “Robust Locally Weighted Regression and Smoothing Scatterplots”. Journal of the American Statistical Association, 74(368), 829-836. DOI: 10.1080/01621459.1979.10481038.

SuperSmoother (supsmu)

greybox.smoothers.supsmu(x, y, wt=None, span='cv', periodic=False, bass=0.0)[source]

Friedman’s variable-span super-smoother (SuperSmoother).

Smooths a scatter-plot using Friedman’s super-smoother algorithm. The smoother evaluates three running linear smoothers (“tweeters”) with spans 0.05, 0.2, and 0.5 of the sample size, then chooses the best span at each abscissa value via cross-validated residuals. The implementation is a direct port of R’s stats::supsmu() (FORTRAN supsmu from ppr.f) and matches R to machine precision.

Parameters:

x (array_like) – The abscissa values. Will be sorted internally if not already ascending.
y (array_like) – The ordinate values, same length as x.
wt (array_like, optional) – Per-observation weights. Default is uniform.
span (float or {"cv"}, default="cv") – Smoother span. Pass a float in (0, 1] for a fixed span; "cv" (or 0) selects the span automatically via leave-one-out cross-validation at each point.
periodic (bool, default=False) – If True, treat x as a periodic variable in [0, 1].
bass (float, default=0.0) – Bass-tone control in [0, 10]. Larger values shift the cross-validated span towards the smoother end of the range (more smoothing in noisy regions). Values outside the range disable the adjustment.

Returns:

Dictionary with two keys:

"x" : numpy.ndarray of sorted abscissa values.
"y" : numpy.ndarray of smoothed ordinate values, aligned with "x".

Return type:

dict

Notes

For small samples (n < 40) or data with substantial serial correlation between observations close in x, a prespecified fixed span (span=0.2 to span=0.4) is recommended over cross-validated selection.

The cross-validation step uses leave-one-out residuals from each of the three tweeters, smooths them with the medium-span smoother, then picks the smallest residual at each abscissa. A final smooth with the smallest span produces the output.

References

Friedman, J. H. (1984). “A Variable Span Smoother”. Technical Report 5 (SLAC-PUB-3477; STAN-LCS-005), Laboratory for Computational Statistics, Department of Statistics, Stanford University. https://www.osti.gov/biblio/1447470

Examples

Smooth count-data noise with the default cross-validated span:

>>> import numpy as np
>>> from greybox import supsmu
>>> rng = np.random.default_rng(1)
>>> x = np.arange(100, dtype=float)
>>> y = 0.05 * x + rng.normal(0, 1, 100)
>>> out = supsmu(x, y)
>>> out["y"].shape
(100,)

Force a fixed span of 0.3:

>>> out_fixed = supsmu(x, y, span=0.3)

Example:

import numpy as np
from greybox import supsmu

rng = np.random.default_rng(1)
x = np.arange(100, dtype=float)
y = 0.05 * x + rng.normal(0, 1.0, 100)

# Default: cross-validated span selection
cv = supsmu(x, y)

# Fixed span (0 < span <= 1)
fixed = supsmu(x, y, span=0.3)

# Bass-tone control increases smoothness in noisy regions
smoother = supsmu(x, y, bass=5.0)

Reference

Friedman, J. H. (1984). “A Variable Span Smoother”. Technical Report 5 (SLAC-PUB-3477; STAN-LCS-005), Laboratory for Computational Statistics, Department of Statistics, Stanford University. OSTI 1447470.

Plotting smoother output

Both lowess() and supsmu() return plain dictionaries, mirroring R’s list(x, y). They can be overlaid on a scatter plot with matplotlib using the standard pattern:

import matplotlib.pyplot as plt
import numpy as np
from greybox import lowess, supsmu

rng = np.random.default_rng(0)
x = np.linspace(0, 6, 80)
y = np.sin(x) + rng.normal(0, 0.2, 80)

plt.scatter(x, y, s=10, alpha=0.5, label="data")
lo = lowess(x, y, f=0.4)
sm = supsmu(x, y)
plt.plot(lo["x"], lo["y"], color="red", label="LOWESS")
plt.plot(sm["x"], sm["y"], color="blue", label="SuperSmoother")
plt.legend()
plt.show()

The same shape works as an overlay on an existing Axes: ax.plot(out["x"], out["y"]).