{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Illustration of prior and posterior Gaussian process for different kernels\n\nThis example illustrates the prior and posterior of a\n:class:`~sklearn.gaussian_process.GaussianProcessRegressor` with different\nkernels. Mean, standard deviation, and 5 samples are shown for both prior\nand posterior distributions.\n\nHere, we only give some illustration. To know more about kernels' formulation,\nrefer to the `User Guide <gp_kernels>`.\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": [
        "## Helper function\n\nBefore presenting each individual kernel available for Gaussian processes,\nwe will define an helper function allowing us plotting samples drawn from\nthe Gaussian process.\n\nThis function will take a\n:class:`~sklearn.gaussian_process.GaussianProcessRegressor` model and will\ndrawn sample from the Gaussian process. If the model was not fit, the samples\nare drawn from the prior distribution while after model fitting, the samples are\ndrawn from the posterior distribution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\n\ndef plot_gpr_samples(gpr_model, n_samples, ax):\n    \"\"\"Plot samples drawn from the Gaussian process model.\n\n    If the Gaussian process model is not trained then the drawn samples are\n    drawn from the prior distribution. Otherwise, the samples are drawn from\n    the posterior distribution. Be aware that a sample here corresponds to a\n    function.\n\n    Parameters\n    ----------\n    gpr_model : `GaussianProcessRegressor`\n        A :class:`~sklearn.gaussian_process.GaussianProcessRegressor` model.\n    n_samples : int\n        The number of samples to draw from the Gaussian process distribution.\n    ax : matplotlib axis\n        The matplotlib axis where to plot the samples.\n    \"\"\"\n    x = np.linspace(0, 5, 100)\n    X = x.reshape(-1, 1)\n\n    y_mean, y_std = gpr_model.predict(X, return_std=True)\n    y_samples = gpr_model.sample_y(X, n_samples)\n\n    for idx, single_prior in enumerate(y_samples.T):\n        ax.plot(\n            x,\n            single_prior,\n            linestyle=\"--\",\n            alpha=0.7,\n            label=f\"Sampled function #{idx + 1}\",\n        )\n    ax.plot(x, y_mean, color=\"black\", label=\"Mean\")\n    ax.fill_between(\n        x,\n        y_mean - y_std,\n        y_mean + y_std,\n        alpha=0.1,\n        color=\"black\",\n        label=r\"$\\pm$ 1 std. dev.\",\n    )\n    ax.set_xlabel(\"x\")\n    ax.set_ylabel(\"y\")\n    ax.set_ylim([-3, 3])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Dataset and Gaussian process generation\nWe will create a training dataset that we will use in the different sections.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rng = np.random.RandomState(4)\nX_train = rng.uniform(0, 5, 10).reshape(-1, 1)\ny_train = np.sin((X_train[:, 0] - 2.5) ** 2)\nn_samples = 5"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Kernel cookbook\n\nIn this section, we illustrate some samples drawn from the prior and posterior\ndistributions of the Gaussian process with different kernels.\n\n### Radial Basis Function kernel\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.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0))\ngpr = GaussianProcessRegressor(kernel=kernel, random_state=0)\n\nfig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))\n\n# plot prior\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])\naxs[0].set_title(\"Samples from prior distribution\")\n\n# plot posterior\ngpr.fit(X_train, y_train)\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])\naxs[1].scatter(X_train[:, 0], y_train, color=\"red\", zorder=10, label=\"Observations\")\naxs[1].legend(bbox_to_anchor=(1.05, 1.5), loc=\"upper left\")\naxs[1].set_title(\"Samples from posterior distribution\")\n\nfig.suptitle(\"Radial Basis Function kernel\", fontsize=18)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Kernel parameters before fit:\\n{kernel})\")\nprint(\n    f\"Kernel parameters after fit: \\n{gpr.kernel_} \\n\"\n    f\"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Rational Quadratic kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.gaussian_process.kernels import RationalQuadratic\n\nkernel = 1.0 * RationalQuadratic(length_scale=1.0, alpha=0.1, alpha_bounds=(1e-5, 1e15))\ngpr = GaussianProcessRegressor(kernel=kernel, random_state=0)\n\nfig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))\n\n# plot prior\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])\naxs[0].set_title(\"Samples from prior distribution\")\n\n# plot posterior\ngpr.fit(X_train, y_train)\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])\naxs[1].scatter(X_train[:, 0], y_train, color=\"red\", zorder=10, label=\"Observations\")\naxs[1].legend(bbox_to_anchor=(1.05, 1.5), loc=\"upper left\")\naxs[1].set_title(\"Samples from posterior distribution\")\n\nfig.suptitle(\"Rational Quadratic kernel\", fontsize=18)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Kernel parameters before fit:\\n{kernel})\")\nprint(\n    f\"Kernel parameters after fit: \\n{gpr.kernel_} \\n\"\n    f\"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Exp-Sine-Squared kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.gaussian_process.kernels import ExpSineSquared\n\nkernel = 1.0 * ExpSineSquared(\n    length_scale=1.0,\n    periodicity=3.0,\n    length_scale_bounds=(0.1, 10.0),\n    periodicity_bounds=(1.0, 10.0),\n)\ngpr = GaussianProcessRegressor(kernel=kernel, random_state=0)\n\nfig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))\n\n# plot prior\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])\naxs[0].set_title(\"Samples from prior distribution\")\n\n# plot posterior\ngpr.fit(X_train, y_train)\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])\naxs[1].scatter(X_train[:, 0], y_train, color=\"red\", zorder=10, label=\"Observations\")\naxs[1].legend(bbox_to_anchor=(1.05, 1.5), loc=\"upper left\")\naxs[1].set_title(\"Samples from posterior distribution\")\n\nfig.suptitle(\"Exp-Sine-Squared kernel\", fontsize=18)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Kernel parameters before fit:\\n{kernel})\")\nprint(\n    f\"Kernel parameters after fit: \\n{gpr.kernel_} \\n\"\n    f\"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Dot-product kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.gaussian_process.kernels import ConstantKernel, DotProduct\n\nkernel = ConstantKernel(0.1, (0.01, 10.0)) * (\n    DotProduct(sigma_0=1.0, sigma_0_bounds=(0.1, 10.0)) ** 2\n)\ngpr = GaussianProcessRegressor(kernel=kernel, random_state=0, normalize_y=True)\n\nfig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))\n\n# plot prior\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])\naxs[0].set_title(\"Samples from prior distribution\")\n\n# plot posterior\ngpr.fit(X_train, y_train)\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])\naxs[1].scatter(X_train[:, 0], y_train, color=\"red\", zorder=10, label=\"Observations\")\naxs[1].legend(bbox_to_anchor=(1.05, 1.5), loc=\"upper left\")\naxs[1].set_title(\"Samples from posterior distribution\")\n\nfig.suptitle(\"Dot-product kernel\", fontsize=18)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Kernel parameters before fit:\\n{kernel})\")\nprint(\n    f\"Kernel parameters after fit: \\n{gpr.kernel_} \\n\"\n    f\"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Mat\u00e9rn kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sklearn.gaussian_process.kernels import Matern\n\nkernel = 1.0 * Matern(length_scale=1.0, length_scale_bounds=(1e-1, 10.0), nu=1.5)\ngpr = GaussianProcessRegressor(kernel=kernel, random_state=0)\n\nfig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))\n\n# plot prior\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])\naxs[0].set_title(\"Samples from prior distribution\")\n\n# plot posterior\ngpr.fit(X_train, y_train)\nplot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])\naxs[1].scatter(X_train[:, 0], y_train, color=\"red\", zorder=10, label=\"Observations\")\naxs[1].legend(bbox_to_anchor=(1.05, 1.5), loc=\"upper left\")\naxs[1].set_title(\"Samples from posterior distribution\")\n\nfig.suptitle(\"Mat\u00e9rn kernel\", fontsize=18)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Kernel parameters before fit:\\n{kernel})\")\nprint(\n    f\"Kernel parameters after fit: \\n{gpr.kernel_} \\n\"\n    f\"Log-likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.3f}\"\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
}