.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\bayesian-optimization.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_bayesian-optimization.py>`
        to download the full example code or to run this example in your browser via Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_bayesian-optimization.py:


==================================
Bayesian optimization with `skopt`
==================================

Gilles Louppe, Manoj Kumar July 2016.
Reformatted by Holger Nahrstaedt 2020

.. currentmodule:: skopt

Problem statement
-----------------

We are interested in solving

.. math::
    x^* = arg \\min_x f(x)

under the constraints that

- :math:`f` is a black box for which no closed form is known
  (nor its gradients);
- :math:`f` is expensive to evaluate;
- and evaluations of :math:`y = f(x)` may be noisy.

**Disclaimer.** If you do not have these constraints, then there
is certainly a better optimization algorithm than Bayesian optimization.

This example uses :class:`plots.plot_gaussian_process` which is available
since version 0.8.

Bayesian optimization loop
--------------------------

For :math:`t=1:T`:

1. Given observations :math:`(x_i, y_i=f(x_i))` for :math:`i=1:t`, build a
   probabilistic model for the objective :math:`f`. Integrate out all
   possible true functions, using Gaussian process regression.

2. optimize a cheap acquisition/utility function :math:`u` based on the
   posterior distribution for sampling the next point.
   :math:`x_{t+1} = arg \\min_x u(x)`
   Exploit uncertainty to balance exploration against exploitation.

3. Sample the next observation :math:`y_{t+1}` at :math:`x_{t+1}`.


Acquisition functions
---------------------

Acquisition functions :math:`u(x)` specify which sample :math:`x`: should be
tried next:

- Expected improvement (default):
  :math:`-EI(x) = -\\mathbb{E} [f(x) - f(x_t^+)]`
- Lower confidence bound: :math:`LCB(x) = \\mu_{GP}(x) + \\kappa \\sigma_{GP}(x)`
- Probability of improvement: :math:`-PI(x) = -P(f(x) \\geq f(x_t^+) + \\kappa)`

where :math:`x_t^+` is the best point observed so far.

In most cases, acquisition functions provide knobs (e.g., :math:`\\kappa`) for
controlling the exploration-exploitation trade-off.
- Search in regions where :math:`\\mu_{GP}(x)` is high (exploitation)
- Probe regions where uncertainty :math:`\\sigma_{GP}(x)` is high (exploration)

.. GENERATED FROM PYTHON SOURCE LINES 67-77

.. code-block:: Python


    print(__doc__)

    import numpy as np

    np.random.seed(237)
    import matplotlib.pyplot as plt

    from skopt.plots import plot_gaussian_process








.. GENERATED FROM PYTHON SOURCE LINES 78-82

Toy example
-----------

Let assume the following noisy function :math:`f`:

.. GENERATED FROM PYTHON SOURCE LINES 82-90

.. code-block:: Python


    noise_level = 0.1


    def f(x, noise_level=noise_level):
        return np.sin(5 * x[0]) * (1 - np.tanh(x[0] ** 2)) + np.random.randn() * noise_level









.. GENERATED FROM PYTHON SOURCE LINES 91-94

**Note.** In `skopt`, functions :math:`f` are assumed to take as input a 1D
vector :math:`x`: represented as an array-like and to return a scalar
:math:`f(x)`:.

.. GENERATED FROM PYTHON SOURCE LINES 94-115

.. code-block:: Python


    # Plot f(x) + contours
    x = np.linspace(-2, 2, 400).reshape(-1, 1)
    fx = [f(x_i, noise_level=0.0) for x_i in x]
    plt.plot(x, fx, "r--", label="True (unknown)")
    plt.fill(
        np.concatenate([x, x[::-1]]),
        np.concatenate(
            (
                [fx_i - 1.9600 * noise_level for fx_i in fx],
                [fx_i + 1.9600 * noise_level for fx_i in fx[::-1]],
            )
        ),
        alpha=0.2,
        fc="r",
        ec="None",
    )
    plt.legend()
    plt.grid()
    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_bayesian-optimization_001.png
   :alt: bayesian optimization
   :srcset: /auto_examples/images/sphx_glr_bayesian-optimization_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 116-118

Bayesian optimization based on gaussian process regression is implemented in
:class:`gp_minimize` and can be carried out as follows:

.. GENERATED FROM PYTHON SOURCE LINES 118-131

.. code-block:: Python


    from skopt import gp_minimize

    res = gp_minimize(
        f,  # the function to minimize
        [(-2.0, 2.0)],  # the bounds on each dimension of x
        acq_func="EI",  # the acquisition function
        n_calls=15,  # the number of evaluations of f
        n_random_starts=5,  # the number of random initialization points
        noise=0.1**2,  # the noise level (optional)
        random_state=1234,
    )  # the random seed








.. GENERATED FROM PYTHON SOURCE LINES 132-133

Accordingly, the approximated minimum is found to be:

.. GENERATED FROM PYTHON SOURCE LINES 133-136

.. code-block:: Python


    f"x^*={res.x[0]:.4f}, f(x^*)={res.fun:.4f}"





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    'x^*=-0.3552, f(x^*)=-1.0079'



.. GENERATED FROM PYTHON SOURCE LINES 137-148

For further inspection of the results, attributes of the `res` named tuple
provide the following information:

- `x` [float]: location of the minimum.
- `fun` [float]: function value at the minimum.
- `models`: surrogate models used for each iteration.
- `x_iters` [array]:
  location of function evaluation for each iteration.
- `func_vals` [array]: function value for each iteration.
- `space` [Space]: the optimization space.
- `specs` [dict]: parameters passed to the function.

.. GENERATED FROM PYTHON SOURCE LINES 148-151

.. code-block:: Python


    print(res)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

              fun: -1.007919274002016
                x: [-0.35518414273753307]
        func_vals: [ 3.716e-02  6.739e-03 ...  8.157e-03 -7.976e-01]
          x_iters: [[-0.009345334109402526], [1.2713537644662787], [0.4484475787090836], [1.0854396754496047], [1.4426790855107496], [0.9579248468740365], [-0.4515808656811222], [-0.6859481043850504], [-0.35518414273753307], [-0.29315377717222235], [-0.32099415298782463], [-2.0], [2.0], [-1.3373742019079444], [-0.24784228664930108]]
           models: [GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775), GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
                                            n_restarts_optimizer=2, noise=0.010000000000000002,
                                            normalize_y=True, random_state=822569775)]
            space: Space([Real(low=-2.0, high=2.0, prior='uniform', transform='normalize')])
     random_state: RandomState(MT19937)
            specs:     args:                    func: <function f at 0x0000020BD78EB060>
                                          dimensions: Space([Real(low=-2.0, high=2.0, prior='uniform', transform='normalize')])
                                      base_estimator: GaussianProcessRegressor(kernel=1**2 * Matern(length_scale=1, nu=2.5),
                                                                               n_restarts_optimizer=2, noise=0.010000000000000002,
                                                                               normalize_y=True, random_state=822569775)
                                             n_calls: 15
                                     n_random_starts: 5
                                    n_initial_points: 10
                             initial_point_generator: random
                                            acq_func: EI
                                       acq_optimizer: auto
                                                  x0: None
                                                  y0: None
                                        random_state: RandomState(MT19937)
                                             verbose: False
                                            callback: None
                                            n_points: 10000
                                n_restarts_optimizer: 5
                                                  xi: 0.01
                                               kappa: 1.96
                                              n_jobs: 1
                                    model_queue_size: None
                                    space_constraint: None
                   function: base_minimize




.. GENERATED FROM PYTHON SOURCE LINES 152-155

Together these attributes can be used to visually inspect the results of the
minimization, such as the convergence trace or the acquisition function at
the last iteration:

.. GENERATED FROM PYTHON SOURCE LINES 155-160

.. code-block:: Python


    from skopt.plots import plot_convergence

    plot_convergence(res)




.. image-sg:: /auto_examples/images/sphx_glr_bayesian-optimization_002.png
   :alt: Convergence plot
   :srcset: /auto_examples/images/sphx_glr_bayesian-optimization_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: title={'center': 'Convergence plot'}, xlabel='Number of calls $n$', ylabel='$\\min f(x)$ after $n$ calls'>



.. GENERATED FROM PYTHON SOURCE LINES 161-165

Let us now visually examine

1. The approximation of the fit gp model to the original function.
2. The acquisition values that determine the next point to be queried.

.. GENERATED FROM PYTHON SOURCE LINES 165-173

.. code-block:: Python


    plt.rcParams["figure.figsize"] = (8, 14)


    def f_wo_noise(x):
        return f(x, noise_level=0)









.. GENERATED FROM PYTHON SOURCE LINES 174-175

Plot the 5 iterations following the 5 random points

.. GENERATED FROM PYTHON SOURCE LINES 175-214

.. code-block:: Python


    for n_iter in range(5):
        # Plot true function.
        plt.subplot(5, 2, 2 * n_iter + 1)

        if n_iter == 0:
            show_legend = True
        else:
            show_legend = False

        ax = plot_gaussian_process(
            res,
            n_calls=n_iter,
            objective=f_wo_noise,
            noise_level=noise_level,
            show_legend=show_legend,
            show_title=False,
            show_next_point=False,
            show_acq_func=False,
        )
        ax.set_ylabel("")
        ax.set_xlabel("")
        # Plot EI(x)
        plt.subplot(5, 2, 2 * n_iter + 2)
        ax = plot_gaussian_process(
            res,
            n_calls=n_iter,
            show_legend=show_legend,
            show_title=False,
            show_mu=False,
            show_acq_func=True,
            show_observations=False,
            show_next_point=True,
        )
        ax.set_ylabel("")
        ax.set_xlabel("")

    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_bayesian-optimization_003.png
   :alt: bayesian optimization
   :srcset: /auto_examples/images/sphx_glr_bayesian-optimization_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 215-232

The first column shows the following:

1. The true function.
2. The approximation to the original function by the gaussian process model
3. How sure the GP is about the function.

The second column shows the acquisition function values after every
surrogate model is fit. It is possible that we do not choose the global
minimum but a local minimum depending on the minimizer used to minimize
the acquisition function.

At the points closer to the points previously evaluated at, the variance
dips to zero.

Finally, as we increase the number of points, the GP model approaches
the actual function. The final few points are clustered around the minimum
because the GP does not gain anything more by further exploration:

.. GENERATED FROM PYTHON SOURCE LINES 232-239

.. code-block:: Python


    plt.rcParams["figure.figsize"] = (6, 4)

    # Plot f(x) + contours
    _ = plot_gaussian_process(res, objective=f_wo_noise, noise_level=noise_level)

    plt.show()



.. image-sg:: /auto_examples/images/sphx_glr_bayesian-optimization_004.png
   :alt: x* = -0.3552, f(x*) = -1.0079
   :srcset: /auto_examples/images/sphx_glr_bayesian-optimization_004.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 3.907 seconds)


.. _sphx_glr_download_auto_examples_bayesian-optimization.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/holgern/scikit-optimize/master?urlpath=lab/tree/notebooks/auto_examples/bayesian-optimization.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: bayesian-optimization.ipynb <bayesian-optimization.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: bayesian-optimization.py <bayesian-optimization.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_