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
|
"""
Create a threshold event
========================
"""
# %%
# Abstract
# --------
#
# We present in this example the creation and the use of a :class:`~openturns.ThresholdEvent` to estimate a simple integral.
#
import openturns as ot
import openturns.viewer as otv
# %%
# We consider a standard Gaussian random vector :math:`X` and build a random vector from this distribution.
distX = ot.Normal()
vecX = ot.RandomVector(distX)
# %%
# We consider the simple model :math:`f:x \mapsto |x|` and consider the output random variable :math:`Y = f(X)`.
f = ot.SymbolicFunction(["x1"], ["abs(x1)"])
vecY = ot.CompositeRandomVector(f, vecX)
# %%
# We define a very simple :class:`~openturns.ThresholdEvent` which happpens whenever :math:`|X| < 1.0` :
thresholdEvent = ot.ThresholdEvent(vecY, ot.Less(), 1.0)
# %%
# For the normal distribution, it is a well-known fact that the values lower than one standard deviation (here exactly 1)
# away from the mean (here 0) account roughly for 68.27% of the set.
# So the probability of the event is:
#
print("Probability of the event : %.4f" % 0.6827)
# %%
# We can also use a basic estimator to get the probability of the event by drawing samples from the initial distribution `distX` and counting those which realize the event:
print(
"Probability of the event (event sampling) : %.4f"
% thresholdEvent.getSample(1000).computeMean()[0]
)
# %%
# The geometric interpretation is simply the area under the PDF of the standard normal distribution for :math:`x \in [-1,1]` which we draw thereafter.
# %%
def linearSample(xmin, xmax, npoints):
"""
Returns a sample created from a regular grid
from xmin to xmax with npoints points.
"""
step = (xmax - xmin) / (npoints - 1)
rg = ot.RegularGrid(xmin, step, npoints)
vertices = rg.getVertices()
return vertices
# %%
# The boundary of the event are the lines :math:`x = -1.0` and :math:`x = 1.0`
a, b = -1, 1
# %%
nplot = 100 # Number of points in the plot
x = linearSample(a, b, nplot)
y = distX.computePDF(x)
vLow = [0.0 for i in range(nplot)]
vUp = [y[i, 0] for i in range(nplot)]
area = distX.computeCDF(b) - distX.computeCDF(a)
boundsPoly = ot.Polygon.FillBetween(x.asPoint(), vLow, vUp)
# %%
# We add the colored area to the PDF graph.
graph = distX.drawPDF()
graph.add(boundsPoly)
graph.setTitle("Probability of the event E = %.4f" % (area))
graph.setLegends([""])
view = otv.View(graph)
# %%
# Display all figures
otv.View.ShowAll()
|
