{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Gaussian Processes regression: basic introductory example\n\nA simple one-dimensional regression example computed in two different ways:\n\n1. A noise-free case\n2. A noisy case with known noise-level per datapoint\n\nIn both cases, the kernel's parameters are estimated using the maximum\nlikelihood principle.\n\nThe figures illustrate the interpolating property of the Gaussian Process model\nas well as its probabilistic nature in the form of a pointwise 95% confidence\ninterval.\n\nNote that `alpha` is a parameter to control the strength of the Tikhonov\nregularization on the assumed training points' covariance matrix.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The scikit-learn developers\n# SPDX-License-Identifier: BSD-3-Clause"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Dataset generation\n\nWe will start by generating a synthetic dataset. The true generative process\nis defined as $f(x) = x \\sin(x)$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nX = np.linspace(start=0, stop=10, num=1_000).reshape(-1, 1)\ny = np.squeeze(X * np.sin(X))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nplt.plot(X, y, label=r\"$f(x) = x \\sin(x)$\", linestyle=\"dotted\")\nplt.legend()\nplt.xlabel(\"$x$\")\nplt.ylabel(\"$f(x)$\")\n_ = plt.title(\"True generative process\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will use this dataset in the next experiment to illustrate how Gaussian\nProcess regression is working.\n\n## Example with noise-free target\n\nIn this first example, we will use the true generative process without\nadding any noise. For training the Gaussian Process regression, we will only\nselect few samples.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rng = np.random.RandomState(1)\ntraining_indices = rng.choice(np.arange(y.size), size=6, replace=False)\nX_train, y_train = X[training_indices], y[training_indices]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we fit a Gaussian process on these few training data samples. We will\nuse a radial basis function (RBF) kernel and a constant parameter to fit the\namplitude.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF\n\nkernel = 1 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2))\ngaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)\ngaussian_process.fit(X_train, y_train)\ngaussian_process.kernel_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After fitting our model, we see that the hyperparameters of the kernel have\nbeen optimized. Now, we will use our kernel to compute the mean prediction\nof the full dataset and plot the 95% confidence interval.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_prediction, std_prediction = gaussian_process.predict(X, return_std=True)\n\nplt.plot(X, y, label=r\"$f(x) = x \\sin(x)$\", linestyle=\"dotted\")\nplt.scatter(X_train, y_train, label=\"Observations\")\nplt.plot(X, mean_prediction, label=\"Mean prediction\")\nplt.fill_between(\n    X.ravel(),\n    mean_prediction - 1.96 * std_prediction,\n    mean_prediction + 1.96 * std_prediction,\n    alpha=0.5,\n    label=r\"95% confidence interval\",\n)\nplt.legend()\nplt.xlabel(\"$x$\")\nplt.ylabel(\"$f(x)$\")\n_ = plt.title(\"Gaussian process regression on noise-free dataset\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We see that for a prediction made on a data point close to the one from the\ntraining set, the 95% confidence has a small amplitude. Whenever a sample\nfalls far from training data, our model's prediction is less accurate and the\nmodel prediction is less precise (higher uncertainty).\n\n## Example with noisy targets\n\nWe can repeat a similar experiment adding an additional noise to the target\nthis time. It will allow seeing the effect of the noise on the fitted model.\n\nWe add some random Gaussian noise to the target with an arbitrary\nstandard deviation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "noise_std = 0.75\ny_train_noisy = y_train + rng.normal(loc=0.0, scale=noise_std, size=y_train.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create a similar Gaussian process model. In addition to the kernel, this\ntime, we specify the parameter `alpha` which can be interpreted as the\nvariance of a Gaussian noise.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gaussian_process = GaussianProcessRegressor(\n    kernel=kernel, alpha=noise_std**2, n_restarts_optimizer=9\n)\ngaussian_process.fit(X_train, y_train_noisy)\nmean_prediction, std_prediction = gaussian_process.predict(X, return_std=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the mean prediction and the uncertainty region as before.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(X, y, label=r\"$f(x) = x \\sin(x)$\", linestyle=\"dotted\")\nplt.errorbar(\n    X_train,\n    y_train_noisy,\n    noise_std,\n    linestyle=\"None\",\n    color=\"tab:blue\",\n    marker=\".\",\n    markersize=10,\n    label=\"Observations\",\n)\nplt.plot(X, mean_prediction, label=\"Mean prediction\")\nplt.fill_between(\n    X.ravel(),\n    mean_prediction - 1.96 * std_prediction,\n    mean_prediction + 1.96 * std_prediction,\n    color=\"tab:orange\",\n    alpha=0.5,\n    label=r\"95% confidence interval\",\n)\nplt.legend()\nplt.xlabel(\"$x$\")\nplt.ylabel(\"$f(x)$\")\n_ = plt.title(\"Gaussian process regression on a noisy dataset\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The noise affects the predictions close to the training samples: the\npredictive uncertainty near to the training samples is larger because we\nexplicitly model a given level target noise independent of the input\nvariable.\n\n"
      ]
    }
  ],
  "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
}