File: plot_bayes_distribution.py

package info (click to toggle)
openturns 1.26-4
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 67,708 kB
  • sloc: cpp: 261,605; python: 67,030; ansic: 4,378; javascript: 406; sh: 185; xml: 164; makefile: 101
file content (60 lines) | stat: -rw-r--r-- 1,555 bytes parent folder | download
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
r"""
Create a Joint by Conditioning  distribution
============================================
"""

# %%
# In this example we are going to build the distribution of the random vector:
#
# .. math::
#    (Y, \vect{X}|\vect{\Theta})
#
# with :math:`\vect{X}` conditioned by the random vector :math:`\vect{\Theta}` obtained with the random variable :math:`Y` through a function :math:`f`:
#
# .. math::
#    \vect{\Theta} = f(\vect{Y})
#

# %%
import openturns as ot
import openturns.viewer as otv


# %%
# We consider the following case: :math:`X|\vect{\Theta} \sim \cN(\vect{\Theta})`
# with :math:`\vect{\Theta} = (Y, 0.1 + Y^2)` and :math:`Y \sim \cN(0,1)`.
#
# We first create the :math:`Y` distribution:
YDist = ot.Normal(0.0, 1.0)

# %%
# Then we create the link function :math:`f: y \rightarrow (y, 0.1+y^2)`:
f = ot.SymbolicFunction(["y"], ["y", "0.1 + y^2"])

# %%
# Then, we create the :math:`\vect{X}|\vect{\Theta}` distribution:
XgivenThetaDist = ot.Normal()

# %%
# At last, we create the distribution of :math:`(Y,X)`:
XDist = ot.JointByConditioningDistribution(XgivenThetaDist, YDist, f)
XDist.setDescription(["Y", r"$X|\mathbf{\boldsymbol{\Theta}} = f(Y)$"])
XDist

# %%
# Get a sample:
sample = XDist.getSample(100)

# %%
# Draw the PDF:
ot.ResourceMap.SetAsString("Contour-DefaultColorMapNorm", "rank")
graph = XDist.drawPDF(sample.getMin(), sample.getMax(), [256] * 2)
graph.setTitle(r"$(Y,X)$ iso-PDF")
cloud = ot.Cloud(sample)
cloud.setColor("red")
cloud.setLegend("sample")
graph.add(cloud)
view = otv.View(graph)

# %%
view.ShowAll()