{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Robust linear estimator fitting\n\nHere a sine function is fit with a polynomial of order 3, for values\nclose to zero.\n\nRobust fitting is demonstrated in different situations:\n\n- No measurement errors, only modelling errors (fitting a sine with a\n  polynomial)\n\n- Measurement errors in X\n\n- Measurement errors in y\n\nThe median absolute deviation to non corrupt new data is used to judge\nthe quality of the prediction.\n\nWhat we can see that:\n\n- RANSAC is good for strong outliers in the y direction\n\n- TheilSen is good for small outliers, both in direction X and y, but has\n  a break point above which it performs worse than OLS.\n\n- The scores of HuberRegressor may not be compared directly to both TheilSen\n  and RANSAC because it does not attempt to completely filter the outliers\n  but lessen their effect.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The scikit-learn developers\n# SPDX-License-Identifier: BSD-3-Clause\n\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\nfrom sklearn.linear_model import (\n    HuberRegressor,\n    LinearRegression,\n    RANSACRegressor,\n    TheilSenRegressor,\n)\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import PolynomialFeatures\n\nnp.random.seed(42)\n\nX = np.random.normal(size=400)\ny = np.sin(X)\n# Make sure that it X is 2D\nX = X[:, np.newaxis]\n\nX_test = np.random.normal(size=200)\ny_test = np.sin(X_test)\nX_test = X_test[:, np.newaxis]\n\ny_errors = y.copy()\ny_errors[::3] = 3\n\nX_errors = X.copy()\nX_errors[::3] = 3\n\ny_errors_large = y.copy()\ny_errors_large[::3] = 10\n\nX_errors_large = X.copy()\nX_errors_large[::3] = 10\n\nestimators = [\n    (\"OLS\", LinearRegression()),\n    (\"Theil-Sen\", TheilSenRegressor(random_state=42)),\n    (\"RANSAC\", RANSACRegressor(random_state=42)),\n    (\"HuberRegressor\", HuberRegressor()),\n]\ncolors = {\n    \"OLS\": \"turquoise\",\n    \"Theil-Sen\": \"gold\",\n    \"RANSAC\": \"lightgreen\",\n    \"HuberRegressor\": \"black\",\n}\nlinestyle = {\"OLS\": \"-\", \"Theil-Sen\": \"-.\", \"RANSAC\": \"--\", \"HuberRegressor\": \"--\"}\nlw = 3\n\nx_plot = np.linspace(X.min(), X.max())\nfor title, this_X, this_y in [\n    (\"Modeling Errors Only\", X, y),\n    (\"Corrupt X, Small Deviants\", X_errors, y),\n    (\"Corrupt y, Small Deviants\", X, y_errors),\n    (\"Corrupt X, Large Deviants\", X_errors_large, y),\n    (\"Corrupt y, Large Deviants\", X, y_errors_large),\n]:\n    plt.figure(figsize=(5, 4))\n    plt.plot(this_X[:, 0], this_y, \"b+\")\n\n    for name, estimator in estimators:\n        model = make_pipeline(PolynomialFeatures(3), estimator)\n        model.fit(this_X, this_y)\n        mse = mean_squared_error(model.predict(X_test), y_test)\n        y_plot = model.predict(x_plot[:, np.newaxis])\n        plt.plot(\n            x_plot,\n            y_plot,\n            color=colors[name],\n            linestyle=linestyle[name],\n            linewidth=lw,\n            label=\"%s: error = %.3f\" % (name, mse),\n        )\n\n    legend_title = \"Error of Mean\\nAbsolute Deviation\\nto Non-corrupt Data\"\n    legend = plt.legend(\n        loc=\"upper right\", frameon=False, title=legend_title, prop=dict(size=\"x-small\")\n    )\n    plt.xlim(-4, 10.2)\n    plt.ylim(-2, 10.2)\n    plt.title(title)\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.21"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}