1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
.. _monte_carlo_simulation:
Monte Carlo simulation
----------------------
Using the probability distribution the probability distribution of a random
vector :math:`\inputRV`, we seek to evaluate the following probability:
.. math::
P_f = \Prob{g\left( \inputRV \right) \leq 0}
Here, :math:`\inputRV` is a random vector, :math:`\model` the function known as *limit state function*
which enables the definition of the event :math:`\cD_f = \{\vect{x} \in \Rset^{\inputDim} \, / \, \model(\inputRV) \le 0\}`.
If we have the set :math:`\left\{ \vect{x}_1,\ldots,\vect{x}_\sampleSize \right\}` of :math:`\sampleSize`
independent samples of the random vector :math:`\inputRV`,
we can estimate :math:`\widehat{P}_f` as follows:
.. math::
\widehat{P}_f = \frac{1}{\sampleSize} \sum_{i=1}^\sampleSize \mathbf{1}_{ \left\{ \model(\vect{x}_i) \leq 0 \right\} }
where :math:`\mathbf{1}_{ \left\{ \model(\vect{x}_i) \leq 0 \right\} }`
describes the indicator function equal to 1 if :math:`\model(\vect{x}_i) \leq 0`
and equal to 0 otherwise; the idea here is in fact to estimate the required
probability by the proportion of cases, among the :math:`\sampleSize` samples of :math:`\inputRV`,
for which the event :math:`\cD_f` occurs.
By the law of large numbers, we know that this estimation converges to the
required value :math:`P_f` as the sample size :math:`\sampleSize` tends to infinity.
The Central Limit Theorem allows one to build an asymptotic confidence interval
using the normal limit distribution as follows:
.. math::
\lim_{\sampleSize\rightarrow\infty}\Prob{P_f\in[\widehat{P}_{f,\inf},\widehat{P}_{f,\sup}]} = \alpha
with
.. math::
\widehat{P}_{f,\inf} = \widehat{P}_f - q_{\alpha}\sqrt{\frac{\widehat{P}_f(1-\widehat{P}_f)}{\sampleSize}},
\widehat{P}_{f,\sup} = \widehat{P}_f + q_{\alpha}\sqrt{\frac{\widehat{P}_f(1-\widehat{P}_f)}{\sampleSize}}
and :math:`q_\alpha` is the :math:`(1+\alpha)/2`-quantile of the standard normal distribution.
One can also use a convergence indicator that is independent of the confidence
level :math:`\alpha`: the coefficient of variation, which is the ratio between the
asymptotic standard deviation of the estimate and its mean value.
It is a relative measure of dispersion given by:
.. math::
\textrm{CV}_{\widehat{P}_f}=\sqrt{ \frac{1-\widehat{P}_f}{\sampleSize \widehat{P}_f}}\simeq\frac{1}{\sqrt{\sampleSize\widehat{P}_f}}\mbox{ for }\widehat{P}_f\ll 1
It must be emphasized that these results are *asymptotic* and as such needs that :math:`\sampleSize\gg 1`.
The convergence to the standard normal distribution is dominated by the skewness
of :math:`\mathbf{1}_{ \left\{ \model(\vect{x}_i) \leq 0 \right\} }`
divided by the sample size :math:`\sampleSize`, it means :math:`\frac{1-2P_f}{\sampleSize\sqrt{P_f(1-P_f)}}`.
In the usual case :math:`P_f\ll 1`, the distribution is highly skewed and this
term is approximately equal to :math:`\frac{1}{\sampleSize\sqrt{P_f}}\simeq\textrm{CV}_{\widehat{P}_f}/\sqrt{\sampleSize}`.
A rule of thumb is that if :math:`\textrm{CV}_{\widehat{P}_f}<0.1`
with :math:`\sampleSize>10^2`, the asymptotic nature of the Central Limit Theorem is not problematic.
.. plot::
import openturns as ot
from matplotlib import pyplot as plt
from openturns.viewer import View
f = ot.SymbolicFunction(['x'], ['17-exp(0.1*(x-1.0))'])
graph = f.draw(0.0, 12.0)
dist = ot.Normal([5.0, 15.0], [1.0, 0.25], ot.IdentityMatrix(2))
N = 1000
sample = dist.getSample(N)
sample1 = ot.Sample(0, 2)
sample2 = ot.Sample(0, 2)
for X in sample:
x, y = X
if f([x])[0] > y:
sample1.add(X)
else:
sample2.add(X)
cloud = ot.Cloud(sample1)
cloud.setColor('green')
cloud.setPointStyle('square')
graph.add(cloud)
cloud = ot.Cloud(sample2)
cloud.setColor('red')
cloud.setPointStyle('square')
graph.add(cloud)
graph.setTitle('Monte Carlo simulation (Pf=0.048, N=1000)')
graph.setLegends(['domain Df', 'simulations'])
graph.setLegendPosition('upper right')
View(graph)
The method is also referred to as Direct sampling, Crude Monte Carlo method, Classical Monte Carlo integration.
.. topic:: API:
- See :class:`~openturns.ProbabilitySimulationAlgorithm`
.. topic:: Examples:
- See :doc:`/auto_reliability_sensitivity/reliability/plot_estimate_probability_monte_carlo`
.. topic:: References:
- Robert C.P., Casella G. (2004). Monte-Carlo Statistical Methods, Springer, ISBN 0-387-21239-6, 2nd ed.
- Rubinstein R.Y. (1981). Simulation and The Monte-Carlo methods, John Wiley \& Sons
|
