.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/probability/monte_carlo/plot_01_monte_carlo_integration.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_api_gallery_probability_monte_carlo_plot_01_monte_carlo_integration.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_probability_monte_carlo_plot_01_monte_carlo_integration.py:


The Monte Carlo method: estimating an integral by random sampling
=================================================================

Ulam and von Neumann's idea: replace an exact calculation by an average
over random samples. For :math:`U \sim \mathrm{Uniform}(0, 1)`,

.. math::

   \int_0^1 e^x\,dx = E[e^U] \approx \frac1n \sum_{i=1}^n e^{U_i},

with an error that shrinks like :math:`1/\sqrt n` regardless of
dimension. The exact value is :math:`e - 1`. Importance sampling and
control variates then reduce the error at the same sample size.

.. GENERATED FROM PYTHON SOURCE LINES 18-25

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.probability.systems.monte_carlo import control_variates_integrate, importance_sampling_integrate, monte_carlo_integrate

    exact = np.e - 1.0








.. GENERATED FROM PYTHON SOURCE LINES 26-28

Plain Monte Carlo: the estimate converges like 1/sqrt(n)
------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 28-43

.. code-block:: Python


    ns = np.unique(np.logspace(1, 5, 25).astype(int))
    estimates = [monte_carlo_integrate(np.exp, 0.0, 1.0, n=int(n), seed=0) for n in ns]
    values = np.array([r.estimate for r in estimates])
    errors = np.array([r.std_error for r in estimates])

    fig, ax = plt.subplots()
    ax.semilogx(ns, values, "o-", label="Monte Carlo estimate")
    ax.fill_between(ns, values - 2 * errors, values + 2 * errors, alpha=0.3, label=r"$\pm 2$ standard errors")
    ax.axhline(exact, color="k", ls="--", label="exact $e - 1$")
    ax.set_xlabel("number of random samples $n$")
    ax.set_ylabel(r"estimate of $\int_0^1 e^x\,dx$")
    ax.set_title("Monte Carlo integration")
    ax.legend()




.. image-sg:: /api/gallery/probability/monte_carlo/images/sphx_glr_plot_01_monte_carlo_integration_001.png
   :alt: Monte Carlo integration
   :srcset: /api/gallery/probability/monte_carlo/images/sphx_glr_plot_01_monte_carlo_integration_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x11cf70d70>



.. GENERATED FROM PYTHON SOURCE LINES 44-46

Three estimators at the same sample size
------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 46-67

.. code-block:: Python


    n = 20000
    plain = monte_carlo_integrate(np.exp, 0.0, 1.0, n=n, seed=0)
    print(f"plain:               estimate={plain.estimate:.5f}, std_error={plain.std_error:.5f}")

    # importance sampling with q(x) = (1 + x) / 1.5 on [0, 1], drawn by inverting its CDF
    importance = importance_sampling_integrate(
        lambda x: np.exp(x) * ((x >= 0.0) & (x <= 1.0)),
        proposal_sampler=lambda rng, size: -1.0 + np.sqrt(1.0 + 3.0 * rng.uniform(size=size)),
        proposal_pdf=lambda x: (1.0 + x) / 1.5,
        n=n,
        seed=0,
    )
    print(f"importance sampling: estimate={importance.estimate:.5f}, std_error={importance.std_error:.5f}")

    # control variate g(x) = x, with E[U] = 0.5 known exactly
    cv = control_variates_integrate(np.exp, lambda x: x, control_mean=0.5, a=0.0, b=1.0, n=n, seed=0)
    print(f"control variate:     estimate={cv.estimate:.5f}, std_error={cv.std_error:.5f}")

    print(f"\nexact value e - 1 = {exact:.5f}")
    print(f"variance reduction: importance {(plain.std_error / importance.std_error) ** 2:.1f}x, control variate {(plain.std_error / cv.std_error) ** 2:.1f}x")




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

 .. code-block:: none

    plain:               estimate=1.72216, std_error=0.00347
    importance sampling: estimate=1.71966, std_error=0.00116
    control variate:     estimate=1.71775, std_error=0.00044

    exact value e - 1 = 1.71828
    variance reduction: importance 9.0x, control variate 60.8x





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

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


.. _sphx_glr_download_api_gallery_probability_monte_carlo_plot_01_monte_carlo_integration.py:

.. only:: html

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

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

      :download:`Download Jupyter notebook: plot_01_monte_carlo_integration.ipynb <plot_01_monte_carlo_integration.ipynb>`

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

      :download:`Download Python source code: plot_01_monte_carlo_integration.py <plot_01_monte_carlo_integration.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_monte_carlo_integration.zip <plot_01_monte_carlo_integration.zip>`


.. only:: html

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

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