.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/reliability/structural/plot_03_metropolis_ulam_monte_carlo.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_reliability_structural_plot_03_metropolis_ulam_monte_carlo.py: Metropolis and Ulam's Monte Carlo method (1949) =============================================== Stanislaw Ulam, recovering from an illness and playing solitaire, wondered whether it was quicker to estimate the chance of winning by dealing many games than by combinatorics. With Nicholas Metropolis (1949) he described the *Monte Carlo method*: estimate a probability or an integral by random sampling. For a structure the failure probability :math:`P(g(X) < 0)` is the fraction of sampled load and resistance pairs that fail (:func:`~engineeringkit.reliability.monte_carlo_failure_probability`). It needs no linearisation, which makes it the check on the Hasofer-Lind index (:func:`~engineeringkit.reliability.hasofer_lind_index`), but its relative error falls only as :math:`1/\sqrt{N p_f}`. Sampling around the Hasofer-Lind design point instead of the mean (importance sampling) puts half the samples in the failure domain. .. GENERATED FROM PYTHON SOURCE LINES 20-25 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from engineeringkit.reliability import hasofer_lind_index, monte_carlo_failure_probability .. GENERATED FROM PYTHON SOURCE LINES 26-33 A nonlinear limit state ----------------------- A steel tie of yield stress :math:`f_y` ~ N(350, 25) MPa and area :math:`A` ~ N(1000, 50) mm² carries a load :math:`P` ~ N(200, 30) kN: :math:`g = f_y A - P`, in kN. The product makes the limit state curved, so FORM's :math:`\Phi(-\beta)` is not exact. .. GENERATED FROM PYTHON SOURCE LINES 33-45 .. code-block:: Python def g(x): return x[0] * x[1] / 1000 - x[2] mean, std = [350.0, 1000.0, 200.0], [25.0, 50.0, 30.0] form = hasofer_lind_index(g, mean, std) reference = monte_carlo_failure_probability(g, mean, std, samples=4_000_000, seed=1949) print(f"Hasofer-Lind beta = {form.beta:.3f}, FORM p_f = {form.failure_probability:.3e}") print(f"Monte Carlo (4e6 samples) p_f = {reference.failure_probability:.3e} +/- {reference.std_error:.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Hasofer-Lind beta = 3.590, FORM p_f = 1.656e-04 Monte Carlo (4e6 samples) p_f = 1.735e-04 +/- 6.6e-06 .. GENERATED FROM PYTHON SOURCE LINES 46-48 Convergence: crude sampling against importance sampling ------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 48-72 .. code-block:: Python sizes = np.logspace(2.5, 5.5, 13).astype(int) crude = [monte_carlo_failure_probability(g, mean, std, samples=n, seed=n) for n in sizes] weighted = [monte_carlo_failure_probability(g, mean, std, samples=n, seed=n, center=form.design_point) for n in sizes] fig, ax = plt.subplots(figsize=(6, 3.8)) for label, runs in (("crude Monte Carlo", crude), ("importance sampling at the design point", weighted)): pf = np.array([run.failure_probability for run in runs]) error = np.array([run.std_error for run in runs]) ax.errorbar(sizes, pf, yerr=1.96 * error, fmt="o-", capsize=3, label=label) ax.axhline(reference.failure_probability, color="k", ls=":", label="reference") ax.axhline(form.failure_probability, color="C3", ls="--", label="FORM") ax.set(xscale="log", xlabel="samples", ylabel="failure probability", ylim=(0, 3 * reference.failure_probability)) ax.legend(fontsize=8) fig.tight_layout() n = sizes[-1] print(f"{n} samples: crude CoV {crude[-1].cov:.3f}, importance sampling CoV {weighted[-1].cov:.4f}") print(f"samples for 10% crude accuracy: about {100 / reference.failure_probability:.0f}") assert abs(weighted[-1].failure_probability - reference.failure_probability) < 4 * np.hypot(weighted[-1].std_error, reference.std_error) assert weighted[-1].cov < crude[-1].cov / 10 assert reference.confidence_interval(0.999)[0] < reference.failure_probability < reference.confidence_interval(0.999)[1] plt.show() .. image-sg:: /api/gallery/reliability/structural/images/sphx_glr_plot_03_metropolis_ulam_monte_carlo_001.png :alt: plot 03 metropolis ulam monte carlo :srcset: /api/gallery/reliability/structural/images/sphx_glr_plot_03_metropolis_ulam_monte_carlo_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 316227 samples: crude CoV 0.120, importance sampling CoV 0.0036 samples for 10% crude accuracy: about 576369 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.426 seconds) .. _sphx_glr_download_api_gallery_reliability_structural_plot_03_metropolis_ulam_monte_carlo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/reliability/structural/plot_03_metropolis_ulam_monte_carlo.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_03_metropolis_ulam_monte_carlo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_03_metropolis_ulam_monte_carlo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_03_metropolis_ulam_monte_carlo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_