{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Polynomial and Spline interpolation\n\nThis example demonstrates how to approximate a function with polynomials up to\ndegree ``degree`` by using ridge regression. We show two different ways given\n``n_samples`` of 1d points ``x_i``:\n\n- :class:`~sklearn.preprocessing.PolynomialFeatures` generates all monomials\n up to ``degree``. This gives us the so called Vandermonde matrix with\n ``n_samples`` rows and ``degree + 1`` columns::\n\n [[1, x_0, x_0 ** 2, x_0 ** 3, ..., x_0 ** degree],\n [1, x_1, x_1 ** 2, x_1 ** 3, ..., x_1 ** degree],\n ...]\n\n Intuitively, this matrix can be interpreted as a matrix of pseudo features\n (the points raised to some power). The matrix is akin to (but different from)\n the matrix induced by a polynomial kernel.\n\n- :class:`~sklearn.preprocessing.SplineTransformer` generates B-spline basis\n functions. A basis function of a B-spline is a piece-wise polynomial function\n of degree ``degree`` that is non-zero only between ``degree+1`` consecutive\n knots. Given ``n_knots`` number of knots, this results in matrix of\n ``n_samples`` rows and ``n_knots + degree - 1`` columns::\n\n [[basis_1(x_0), basis_2(x_0), ...],\n [basis_1(x_1), basis_2(x_1), ...],\n ...]\n\nThis example shows that these two transformers are well suited to model\nnon-linear effects with a linear model, using a pipeline to add non-linear\nfeatures. Kernel methods extend this idea and can induce very high (even\ninfinite) dimensional feature spaces.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Author: Mathieu Blondel\n# Jake Vanderplas\n# Christian Lorentzen\n# Malte Londschien\n# License: BSD 3 clause\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom sklearn.linear_model import Ridge\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import PolynomialFeatures, SplineTransformer"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start by defining a function that we intend to approximate and prepare\nplotting it.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def f(x):\n \"\"\"Function to be approximated by polynomial interpolation.\"\"\"\n return x * np.sin(x)\n\n\n# whole range we want to plot\nx_plot = np.linspace(-1, 11, 100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make it interesting, we only give a small subset of points to train on.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x_train = np.linspace(0, 10, 100)\nrng = np.random.RandomState(0)\nx_train = np.sort(rng.choice(x_train, size=20, replace=False))\ny_train = f(x_train)\n\n# create 2D-array versions of these arrays to feed to transformers\nX_train = x_train[:, np.newaxis]\nX_plot = x_plot[:, np.newaxis]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready to create polynomial features and splines, fit on the\ntraining points and show how well they interpolate.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# plot function\nlw = 2\nfig, ax = plt.subplots()\nax.set_prop_cycle(\n color=[\"black\", \"teal\", \"yellowgreen\", \"gold\", \"darkorange\", \"tomato\"]\n)\nax.plot(x_plot, f(x_plot), linewidth=lw, label=\"ground truth\")\n\n# plot training points\nax.scatter(x_train, y_train, label=\"training points\")\n\n# polynomial features\nfor degree in [3, 4, 5]:\n model = make_pipeline(PolynomialFeatures(degree), Ridge(alpha=1e-3))\n model.fit(X_train, y_train)\n y_plot = model.predict(X_plot)\n ax.plot(x_plot, y_plot, label=f\"degree {degree}\")\n\n# B-spline with 4 + 3 - 1 = 6 basis functions\nmodel = make_pipeline(SplineTransformer(n_knots=4, degree=3), Ridge(alpha=1e-3))\nmodel.fit(X_train, y_train)\n\ny_plot = model.predict(X_plot)\nax.plot(x_plot, y_plot, label=\"B-spline\")\nax.legend(loc=\"lower center\")\nax.set_ylim(-20, 10)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This shows nicely that higher degree polynomials can fit the data better. But\nat the same time, too high powers can show unwanted oscillatory behaviour\nand are particularly dangerous for extrapolation beyond the range of fitted\ndata. This is an advantage of B-splines. They usually fit the data as well as\npolynomials and show very nice and smooth behaviour. They have also good\noptions to control the extrapolation, which defaults to continue with a\nconstant. Note that most often, you would rather increase the number of knots\nbut keep ``degree=3``.\n\nIn order to give more insights into the generated feature bases, we plot all\ncolumns of both transformers separately.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(ncols=2, figsize=(16, 5))\npft = PolynomialFeatures(degree=3).fit(X_train)\naxes[0].plot(x_plot, pft.transform(X_plot))\naxes[0].legend(axes[0].lines, [f\"degree {n}\" for n in range(4)])\naxes[0].set_title(\"PolynomialFeatures\")\n\nsplt = SplineTransformer(n_knots=4, degree=3).fit(X_train)\naxes[1].plot(x_plot, splt.transform(X_plot))\naxes[1].legend(axes[1].lines, [f\"spline {n}\" for n in range(6)])\naxes[1].set_title(\"SplineTransformer\")\n\n# plot knots of spline\nknots = splt.bsplines_[0].t\naxes[1].vlines(knots[3:-3], ymin=0, ymax=0.8, linestyles=\"dashed\")\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the left plot, we recognize the lines corresponding to simple monomials\nfrom ``x**0`` to ``x**3``. In the right figure, we see the six B-spline\nbasis functions of ``degree=3`` and also the four knot positions that were\nchosen during ``fit``. Note that there are ``degree`` number of additional\nknots each to the left and to the right of the fitted interval. These are\nthere for technical reasons, so we refrain from showing them. Every basis\nfunction has local support and is continued as a constant beyond the fitted\nrange. This extrapolating behaviour could be changed by the argument\n``extrapolation``.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Periodic Splines\nIn the previous example we saw the limitations of polynomials and splines for\nextrapolation beyond the range of the training observations. In some\nsettings, e.g. with seasonal effects, we expect a periodic continuation of\nthe underlying signal. Such effects can be modelled using periodic splines,\nwhich have equal function value and equal derivatives at the first and last\nknot. In the following case we show how periodic splines provide a better fit\nboth within and outside of the range of training data given the additional\ninformation of periodicity. The splines period is the distance between\nthe first and last knot, which we specify manually.\n\nPeriodic splines can also be useful for naturally periodic features (such as\nday of the year), as the smoothness at the boundary knots prevents a jump in\nthe transformed values (e.g. from Dec 31st to Jan 1st). For such naturally\nperiodic features or more generally features where the period is known, it is\nadvised to explicitly pass this information to the `SplineTransformer` by\nsetting the knots manually.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def g(x):\n \"\"\"Function to be approximated by periodic spline interpolation.\"\"\"\n return np.sin(x) - 0.7 * np.cos(x * 3)\n\n\ny_train = g(x_train)\n\n# Extend the test data into the future:\nx_plot_ext = np.linspace(-1, 21, 200)\nX_plot_ext = x_plot_ext[:, np.newaxis]\n\nlw = 2\nfig, ax = plt.subplots()\nax.set_prop_cycle(color=[\"black\", \"tomato\", \"teal\"])\nax.plot(x_plot_ext, g(x_plot_ext), linewidth=lw, label=\"ground truth\")\nax.scatter(x_train, y_train, label=\"training points\")\n\nfor transformer, label in [\n (SplineTransformer(degree=3, n_knots=10), \"spline\"),\n (\n SplineTransformer(\n degree=3,\n knots=np.linspace(0, 2 * np.pi, 10)[:, None],\n extrapolation=\"periodic\",\n ),\n \"periodic spline\",\n ),\n]:\n model = make_pipeline(transformer, Ridge(alpha=1e-3))\n model.fit(X_train, y_train)\n y_plot_ext = model.predict(X_plot_ext)\n ax.plot(x_plot_ext, y_plot_ext, label=label)\n\nax.legend()\nfig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\nknots = np.linspace(0, 2 * np.pi, 4)\nsplt = SplineTransformer(knots=knots[:, None], degree=3, extrapolation=\"periodic\").fit(\n X_train\n)\nax.plot(x_plot_ext, splt.transform(X_plot_ext))\nax.legend(ax.lines, [f\"spline {n}\" for n in range(3)])\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.18"
}
},
"nbformat": 4,
"nbformat_minor": 0
}