An extended version of the loess method

Estimates local regression estimates of y and dy/dx at specified x locations.

My contribution here is to provide an extended version of the loess method that allows for estimates at any location within the range of the observed data, and also for estimates of the local slope at those locations. The algorithm also accounts for gaps in the observed data.

The loess method was first introduced by William Cleveland between 1979 and 1996 (see Wikipedia page for details). The method is one of several methods included in the Matlab's smooth function. The method can be used for n-dimensional estimation, but the focus here, and also with the smooth function, is on 1D problems.

Consider a set of observations, yi and xi, where i = 1 to n. The observations of yi are assumed be equal to Yi + ei, where the vector Yi represents "true" values, and the vector ei represents random independent errors.

The loess method estimates the function Y(x) by using a "moving" least-squares fit. The moving fit is applied at specified x locations, and is represented by a 2nd order polynomial fitting function, and tricubic weighting function, both of which are centered on the x location for each estimate. As such, the loess method can be viewed as a moving Taylor-series approximation. A nice aspect of the loess method is that, since it is based on the least squares method, all of the knowledge for that method is directly applicable for understanding and interpreting the loess estimates.

The loess function provided here is generalized to allow for more applications. Usually the loess method only provides estimates at the locations xi for the observations. The function here allows for estimates at any specified location, xq, within the range of the observations. The output includes estimates of Yq and dYq/dx at the specified xq locations. In addition, the code also accounts for gaps in the observations where Yq and dYq/dx cannot be directly estimated.

Note that polynomial least-squares is prone to numerical instabilities, given that the transformation that there are powers of x involved. The least-squares fit in this function is done using the QR decomposition, which should ensure a stable calculation for all problems.

Mark Brandon, Yale University, August 2, 2020