{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n=========================================================\nGaussian Processes regression: basic introductory example\n=========================================================\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\nmodel as well as its probabilistic nature in the form of a pointwise 95%\nconfidence interval.\n\nNote that the parameter ``alpha`` is applied as a Tikhonov\nregularization of the assumed covariance between the training points.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n\n# Author: Vincent Dubourg <vincent.dubourg@gmail.com>\n#         Jake Vanderplas <vanderplas@astro.washington.edu>\n#         Jan Hendrik Metzen <jhm@informatik.uni-bremen.de>s\n# License: BSD 3 clause\n\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n\nnp.random.seed(1)\n\n\ndef f(x):\n    \"\"\"The function to predict.\"\"\"\n    return x * np.sin(x)\n\n# ----------------------------------------------------------------------\n#  First the noiseless case\nX = np.atleast_2d([1., 3., 5., 6., 7., 8.]).T\n\n# Observations\ny = f(X).ravel()\n\n# Mesh the input space for evaluations of the real function, the prediction and\n# its MSE\nx = np.atleast_2d(np.linspace(0, 10, 1000)).T\n\n# Instantiate a Gaussian Process model\nkernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))\ngp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)\n\n# Fit to data using Maximum Likelihood Estimation of the parameters\ngp.fit(X, y)\n\n# Make the prediction on the meshed x-axis (ask for MSE as well)\ny_pred, sigma = gp.predict(x, return_std=True)\n\n# Plot the function, the prediction and the 95% confidence interval based on\n# the MSE\nplt.figure()\nplt.plot(x, f(x), 'r:', label=u'$f(x) = x\\,\\sin(x)$')\nplt.plot(X, y, 'r.', markersize=10, label=u'Observations')\nplt.plot(x, y_pred, 'b-', label=u'Prediction')\nplt.fill(np.concatenate([x, x[::-1]]),\n         np.concatenate([y_pred - 1.9600 * sigma,\n                        (y_pred + 1.9600 * sigma)[::-1]]),\n         alpha=.5, fc='b', ec='None', label='95% confidence interval')\nplt.xlabel('$x$')\nplt.ylabel('$f(x)$')\nplt.ylim(-10, 20)\nplt.legend(loc='upper left')\n\n# ----------------------------------------------------------------------\n# now the noisy case\nX = np.linspace(0.1, 9.9, 20)\nX = np.atleast_2d(X).T\n\n# Observations and noise\ny = f(X).ravel()\ndy = 0.5 + 1.0 * np.random.random(y.shape)\nnoise = np.random.normal(0, dy)\ny += noise\n\n# Instantiate a Gaussian Process model\ngp = GaussianProcessRegressor(kernel=kernel, alpha=dy ** 2,\n                              n_restarts_optimizer=10)\n\n# Fit to data using Maximum Likelihood Estimation of the parameters\ngp.fit(X, y)\n\n# Make the prediction on the meshed x-axis (ask for MSE as well)\ny_pred, sigma = gp.predict(x, return_std=True)\n\n# Plot the function, the prediction and the 95% confidence interval based on\n# the MSE\nplt.figure()\nplt.plot(x, f(x), 'r:', label=u'$f(x) = x\\,\\sin(x)$')\nplt.errorbar(X.ravel(), y, dy, fmt='r.', markersize=10, label=u'Observations')\nplt.plot(x, y_pred, 'b-', label=u'Prediction')\nplt.fill(np.concatenate([x, x[::-1]]),\n         np.concatenate([y_pred - 1.9600 * sigma,\n                        (y_pred + 1.9600 * sigma)[::-1]]),\n         alpha=.5, fc='b', ec='None', label='95% confidence interval')\nplt.xlabel('$x$')\nplt.ylabel('$f(x)$')\nplt.ylim(-10, 20)\nplt.legend(loc='upper left')\n\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.6.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}