robustpolyfit 0.0.6
Outlier-robust regression for polynomials 1 <= degree <= 3

# outlier_robust_polyfit

- [Summary](#summary)
- [Installation](#installation)
- [Usage](#usage)
- [Algorithm details, comparison with alternatives](#algorithm)

## Summary

This package fits linear, quadratic or cubic polynomials with Student's 
T-distributed errors using the EM algorithm. It's 2x as fast as
Scipy's siegelslopes for fitting robust linear models to small datasets
(10 - 20 datapoints) and an order of magnitude faster for larger 
datasets (>=100 datapoints). It yields a robust fit for datasets
where 25% of datapoints are outliers (use with caution for datasets
where the "outlier" population may be greater than 25%).

![example](https://raw.githubusercontent.com/jlparki/outlier_robust_polyfit/main/resources/example_1.png)

![example](https://raw.githubusercontent.com/jlparki/outlier_robust_polyfit/main/resources/example_2.png)

## Installation

    pip install robustpolyfit

## Usage

    from RobustPolyfit import robust_polyfit
    class robust_polyfit(int max_iter=500, float tol=1e-2,
    int polyorder = 1, int df = 1)

- Parameters
    - max_iter: The maximum number of iterations
    - tol: Tolerance for convergence (change in lower bound).
Setting this to a higher value reduces the number of iterations, leading
to faster convergence (and in some cases a poorer fit).
    - polyorder: The degree of the polynomial. Allowed values are 1, 2 and 3.
    - df: The degrees of freedom of the Student's t-distribution. A 
smaller value leads to a higher level of tolerance for outliers. df=1
is the smallest allowed value. Large values for df, e.g. >> 5, are essentially
equivalent to normally distributed error and do not offer any benefit compared
to standard least squares fitting.


    fit(x, y)

    - Fits the model using the x and y data supplied.
Both x and y must be numpy np.float64 one-dimensional arrays. If there are
any problems or the fit does not converge, a ValueError will result.

   predict(x)

    - Generates predictions for the x numpy array supplied.

   get_coefs()

    - Gets the coefficients from the fit. They are returned as a numpy array
ordered from lowest to highest degree terms (intercept, first degree term,
second degree term etc).

#### Example

    from RobustPolyfit import robust_polyfit
    from scipy.stats import norm
    import matplotlib.pyplot as plt, numpy as np

    #Create a model using default settings, except for polyorder, which specifies
    #the degree of the polynomial.
    robmodel = robust_polyfit(polyorder=2)

    #Create some data to which to fit this model. We'll create a noisy quadratic
    #with 80 datapoints, then add another 20 generated by something completely different.
    x1 = np.linspace(-10, 10, 80)
    y1 = x1**2 + 5*x1 - 2.5

    percent_dev = norm.rvs(loc=1, scale=0.2, size=80, random_state = 123)
    y1 = y1 * percent_dev

    x2 = np.linspace(-10, 10, 20)
    y2 = 25.0*x2 - 10.5

    x = np.concatenate([x1, x2])
    y = np.concatenate([y1, y2])
    idx = np.argsort(x)
    x, y = x[idx], y[idx]

    #Fit using the robustpolyfit model, then using numpy's polyfit.
    robmodel.fit(x, y)
    rob_preds = robmodel.predict(x)
    rob_coefs = robmodel.get_coefs()

    lstsq_coefs = np.polyfit(x, y, deg=2)
    lstsq_preds = lstsq_coefs[0]*x**2 + lstsq_coefs[1]*x + lstsq_coefs[2]

    plt.scatter(x, y, s=5, label="Raw data")
    plt.plot(x, rob_preds, linestyle="dashed", color="black", label="Robust regression")
    plt.plot(x, lstsq_preds, linestyle="dashed", color="red", label="Least squares")
    plt.title("Robust regression vs least squares for a quadratic")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.show()
    
![example](https://raw.githubusercontent.com/jlparki/outlier_robust_polyfit/main/resources/example_code.png)


## Algorithm

Least squares assumes errors are normally distributed and thus misbehaves badly in the presence of outliers. If there are outliers, a more appropriate assumption is the heavy-tailed Student's t-distribution. Linear models with Student's t-distributed errors do not fit in closed form, so this package fits using the expectation-maximization or EM algorithm. We can do so because an equivalent form for a Student t-distribution is a Gaussian scale mixture:

![equation 1](https://raw.githubusercontent.com/jlparki/outlier_robust_polyfit/main/resources/eq1.png)

Then zi is the hidden variable, and we can use the usual EM approach (see [1] among many others for details). On the E-step, we update the zi values, and on the M step we update the mean using weighted linear regression with the zis as weights and update the variance. Because the zi values in the E-step are
determined using the Mahalanobis distance from the mean of the distribution, outliers are
iteratively down-weighted.

While the EM algorithm has good convergence guarantees, this is obviously much slower than
least squares, because for least squares you perform a QR decomposition once, whereas here
you do so on each iteration. _Usually_ the model converges in five iterations or less, but
there are cases where this is not true, and even then you can anticipate this will be 3-5
times slower than least squares for an equivalent number of datapoints. Also, the EM algorithm is sensitive to choice of starting point. Here we are selecting a starting point by using a simple least squares fit. This works well for typical use cases but may break down if a large fraction of the data lies well outside the expected region.

Compared to Scipy and Scikit-learn's robust linear regression algorithms, however, this one is much faster for any reasonable number of datapoints (partly thanks to some Cython-based optimization, and partly thanks to the choice of algorithm). For 25 datapoints, RobustPolyfit is 2x faster than Scipy's siegelslopes, and siegelslopes exhibits O(N^2) scaling:

![equation 1](https://raw.githubusercontent.com/jlparki/outlier_robust_polyfit/main/resources/time_comp.png)

The use-case is for fitting a linear, quadratic or cubic where outliers may be present but are anticipated to represent a modest fraction of the data. You should be aware that this packge does not zero-center the data and you may want to do that for numerical stability if you are fitting a quadratic or cubic with large x-values.

[1] Kevin P. Murphy. _Machine Learning: A Probabilistic Perspective._ Cambridge, MA: The MIT Press, 2012.