.. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_reliability_sensitivity_reliability_processes_plot_estimate_probability_monte_carlo_process.py: Estimate a process event probability ==================================== The objective of this example is to evaluate the probability of an event based on a stochastic process, using the Monte Carlo estimator. Let :math:`X: \Omega \times \mathcal{D} \rightarrow \mathbb{R}^d` be a stochastic process of dimension :math:`d`, where :math:`\mathcal{D} \in \mathbb{R}^n` is discretized on the mesh :math:`\mathcal{M}`. We define the event :math:`\mathcal{E}` as: .. math:: \begin{aligned} \displaystyle \mathcal{E}(X) = \bigcup_{\underline{t}\in \mathcal{M}}\left\{X_{\underline{t}} \in \mathcal{A} \right\} \end{aligned} where :math:`\mathcal{A}` is a domain of :math:`\mathbb{R}^d`. We estimate the probabilty :math:`p=\mathbb{P}\left(\mathcal{E}(X)\right)` with the Monte Carlo estimator. The Monte Carlo algorithm is manipulated the same way as in the case where the event is based on a random variable independent of time. We illustrate the algorithm on the example of the bidimensionnal white noise process :math:`\varepsilon: \Omega \times \mathcal{D} \rightarrow \mathbb{R}^2` where :math:`\mathcal{D}\in \mathbb{R}`, distributed according to the bidimensionnal standard normal distribution (with zero mean, unit variance and independent marginals). We consider the domain :math:`\mathcal{A} = [1,2] \times [1,2]`. Then the event :math:`\mathcal{E}` writes: .. math:: \begin{aligned} \displaystyle \mathcal{E}(\varepsilon) = \bigcup_{\underline{t}\in \mathcal{M}}\left\{\varepsilon_{t} \in \mathcal{A} \right\} \end{aligned} For all time stamps :math:`t \in \mathcal{M}`, the probability :math:`p_1` that the process enters into the domain :math:`\mathcal{A}` at time :math:`t` writes, using the independence property of the marginals: .. math:: \begin{aligned} p_1 = \mathbb{P}\left(\varepsilon_t \in \mathcal{A}\right) = (\Phi(2) - \Phi(1))^2 \end{aligned} with :math:`\Phi` the cumulative distribution function of the scalar standard *Normal* distribution. As the proces is discretized on a time grid of size :math:`N` and using the independance property of the white noise between two different time stamps and the fact that the white noise follows the same distribution at each time :math:`t`, the final probability :math:`p` writes: .. math:: p = \mathbb{P}\left(\mathcal{E}(\varepsilon)\right) = 1 - (1 - p_1)^{N} With :math:`K=10^4` realizations, using the Monte Carlo estimator, we obtain :math:`p_K = 0.1627`, to be compared to the exact value :math:`p=0.17008` for a time grid of size :math:`N=10`. .. code-block:: default from __future__ import print_function import openturns as ot import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) Create a time grid .. code-block:: default tMin = 0.0 timeStep = 0.1 n = 100 tgrid = ot.RegularGrid(tMin, timeStep, n) Create a normal process .. code-block:: default amplitude = [5.0] scale = [3.0] model = ot.ExponentialModel(scale, amplitude) process = ot.GaussianProcess(model, tgrid) Create the 1-d domain A: [2.,5.] .. code-block:: default lowerBound = [2.0] upperBound = [5.0] domain = ot.Interval(lowerBound, upperBound) Create an event from a Process and a Domain .. code-block:: default event = ot.ProcessEvent(process, domain) Create the Monte-Carlo algorithm .. code-block:: default montecarlo = ot.ProbabilitySimulationAlgorithm(event) # Define the maximum number of simulations montecarlo.setMaximumOuterSampling(1000) # Define the block size montecarlo.setBlockSize(100) # Define the maximum coefficient of variation montecarlo.setMaximumCoefficientOfVariation(0.0025) # Run the algorithm montecarlo.run() # Get the result montecarlo.getResult() .. raw:: html

probabilityEstimate=9.331624e-01 varianceEstimate=5.384142e-06 standard deviation=2.32e-03 coefficient of variation=2.49e-03 confidenceLength(0.95)=9.10e-03 outerSampling=117 blockSize=100



.. code-block:: default graph = montecarlo.drawProbabilityConvergence(0.95) view = viewer.View(graph) plt.show() .. image:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_estimate_probability_monte_carlo_process_001.png :alt: ProbabilitySimulationAlgorithm convergence graph at level 0.95 :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.198 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_processes_plot_estimate_probability_monte_carlo_process.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_probability_monte_carlo_process.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_estimate_probability_monte_carlo_process.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_