.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user/_auto_howtos/RandomVariablesAndAlgebra.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_user__auto_howtos_RandomVariablesAndAlgebra.py: Random Variables and Algebra in CUQIpy ====================================== CUQIpy provides a simple algebraic framework for defining and manipulating random variables. In this example, we demonstrate how to define random variables, apply algebraic operations on them, and finally use them in Bayesian Problems. .. GENERATED FROM PYTHON SOURCE LINES 13-15 Defining Random Variables ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-26 .. code-block:: Python # Random variables can be defined by either initialising the RandomVariable class # with a distribution object or by retrieving the `rv` attribute of a distribution. # The distribution object can be any distribution from the `cuqi.distribution` module. from cuqi.distribution import Normal from cuqi.experimental.algebra import RandomVariable x = RandomVariable(Normal(0, 1)) y = Normal(0, 1).rv .. GENERATED FROM PYTHON SOURCE LINES 27-32 Recording Algebraic Operations ------------------------------ We can now perform some algebraic operations on the random variables. The operations are recorded in a computational graph, which can be evaluated later. .. GENERATED FROM PYTHON SOURCE LINES 32-39 .. code-block:: Python print("Basic operations: \n") print(f"x + y yields:\n{x + y}\n") print(f"x - y yields:\n{x - y}\n") print(f"x * y yields:\n{x * y}\n") print(f"x / y yields:\n{x / y}\n") .. rst-class:: sphx-glr-script-out .. code-block:: none Basic operations: x + y yields: Transformed Random Variable Expression: x + y Components: x ~ CUQI Normal. y ~ CUQI Normal. x - y yields: Transformed Random Variable Expression: x - y Components: x ~ CUQI Normal. y ~ CUQI Normal. x * y yields: Transformed Random Variable Expression: x * y Components: x ~ CUQI Normal. y ~ CUQI Normal. x / y yields: Transformed Random Variable Expression: x / y Components: x ~ CUQI Normal. y ~ CUQI Normal. .. GENERATED FROM PYTHON SOURCE LINES 40-44 .. code-block:: Python print("Complex operations: \n") print(f"x**2 + 2*x*y + y**2 yields:\n{x**2 + 2*x*y + y**2}\n") print(f"(x + y)**2 yields\n{(x + y)**2}\n") .. rst-class:: sphx-glr-script-out .. code-block:: none Complex operations: x**2 + 2*x*y + y**2 yields: Transformed Random Variable Expression: x^2 + (x * 2) * y + y^2 Components: x ~ CUQI Normal. y ~ CUQI Normal. (x + y)**2 yields Transformed Random Variable Expression: (x + y)^2 Components: x ~ CUQI Normal. y ~ CUQI Normal. .. GENERATED FROM PYTHON SOURCE LINES 45-48 .. code-block:: Python print("Array operations: \n") print(f"x[0] + y[1] yields:\n{x[0] + y[1]}\n") .. rst-class:: sphx-glr-script-out .. code-block:: none Array operations: x[0] + y[1] yields: Transformed Random Variable Expression: x[0] + y[1] Components: x ~ CUQI Normal. y ~ CUQI Normal. .. GENERATED FROM PYTHON SOURCE LINES 49-53 Utilizing the recorded operations --------------------------------- We can evaluate the recorded operations by calling the random variable object with the desired values for the random variables. .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: Python # Define a new random variable 'z' z = (x + y)**2 # Evaluate the expression (using the __call__ method) print(f"z={z.expression} evaluated at x=1, y=2 yields: {z(x=1, y=2)}") .. rst-class:: sphx-glr-script-out .. code-block:: none z=(x + y)^2 evaluated at x=1, y=2 yields: 9 .. GENERATED FROM PYTHON SOURCE LINES 60-64 Building Bayesian Problems -------------------------- Random variables can be used to define Bayesian problems. In this example we build an example Bayesian problem using the Deconvolution1D test problem. .. GENERATED FROM PYTHON SOURCE LINES 64-112 .. code-block:: Python from cuqi.testproblem import Deconvolution1D from cuqi.distribution import Gaussian, Gamma, GMRF from cuqi.experimental.algebra import RandomVariable from cuqi.problem import BayesianProblem from cuqi.distribution import JointDistribution from cuqi.experimental.mcmc import HybridGibbs, LinearRTO, Conjugate, ConjugateApprox import numpy as np import matplotlib.pyplot as plt # Forward model A, y_obs, info = Deconvolution1D().get_components() # Bayesian Problem (defined using Random Variables) d = Gamma(1, 1e-4).rv s = Gamma(1, 1e-4).rv x = GMRF(np.zeros(A.domain_dim), d).rv y = Gaussian(A @ x, 1/s).rv # Combine into a Bayesian Problem and perform UQ BP = BayesianProblem(y, x, s, d) BP.set_data(y=y_obs) BP.UQ(exact={"x": info.exactSolution}) # Random variables can also be used to define JointDistribution. Here we solve the same # problem above by explictly forming a target distribution and then drawing samples with # the HybridGibbs sampler. target = JointDistribution(y, x, s, d)(y=y_obs) # Sampling strategy sampling_strategy = { "x" : LinearRTO(), "s" : Conjugate(), "d" : Conjugate() } # Gibbs sampler sampler = HybridGibbs(target, sampling_strategy) # Run sampler sampler.warmup(200) sampler.sample(1000) samples = sampler.get_samples() # Plot plt.figure() samples["x"].plot_ci(exact=info.exactSolution) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_001.png :alt: RandomVariablesAndAlgebra :srcset: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_001.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_002.png :alt: s, s :srcset: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_002.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_003.png :alt: d, d :srcset: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_003.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_004.png :alt: RandomVariablesAndAlgebra :srcset: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Computing 1000 samples !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Automatic sampler selection is a work-in-progress. !!! !!! Always validate the computed results. !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Using Gibbs sampler burn-in: 20% Automatically determined sampling strategy: x: LinearRTO s: Conjugate d: Conjugate Warmup 2 / 200 Warmup 4 / 200 Warmup 6 / 200 Warmup 8 / 200 Warmup 10 / 200 Warmup 12 / 200 Warmup 14 / 200 Warmup 16 / 200 Warmup 18 / 200 Warmup 20 / 200 Warmup 22 / 200 Warmup 24 / 200 Warmup 26 / 200 Warmup 28 / 200 Warmup 30 / 200 Warmup 32 / 200 Warmup 34 / 200 Warmup 36 / 200 Warmup 38 / 200 Warmup 40 / 200 Warmup 42 / 200 Warmup 44 / 200 Warmup 46 / 200 Warmup 48 / 200 Warmup 50 / 200 Warmup 52 / 200 Warmup 54 / 200 Warmup 56 / 200 Warmup 58 / 200 Warmup 60 / 200 Warmup 62 / 200 Warmup 64 / 200 Warmup 66 / 200 Warmup 68 / 200 Warmup 70 / 200 Warmup 72 / 200 Warmup 74 / 200 Warmup 76 / 200 Warmup 78 / 200 Warmup 80 / 200 Warmup 82 / 200 Warmup 84 / 200 Warmup 86 / 200 Warmup 88 / 200 Warmup 90 / 200 Warmup 92 / 200 Warmup 94 / 200 Warmup 96 / 200 Warmup 98 / 200 Warmup 100 / 200 Warmup 102 / 200 Warmup 104 / 200 Warmup 106 / 200 Warmup 108 / 200 Warmup 110 / 200 Warmup 112 / 200 Warmup 114 / 200 Warmup 116 / 200 Warmup 118 / 200 Warmup 120 / 200 Warmup 122 / 200 Warmup 124 / 200 Warmup 126 / 200 Warmup 128 / 200 Warmup 130 / 200 Warmup 132 / 200 Warmup 134 / 200 Warmup 136 / 200 Warmup 138 / 200 Warmup 140 / 200 Warmup 142 / 200 Warmup 144 / 200 Warmup 146 / 200 Warmup 148 / 200 Warmup 150 / 200 Warmup 152 / 200 Warmup 154 / 200 Warmup 156 / 200 Warmup 158 / 200 Warmup 160 / 200 Warmup 162 / 200 Warmup 164 / 200 Warmup 166 / 200 Warmup 168 / 200 Warmup 170 / 200 Warmup 172 / 200 Warmup 174 / 200 Warmup 176 / 200 Warmup 178 / 200 Warmup 180 / 200 Warmup 182 / 200 Warmup 184 / 200 Warmup 186 / 200 Warmup 188 / 200 Warmup 190 / 200 Warmup 192 / 200 Warmup 194 / 200 Warmup 196 / 200 Warmup 198 / 200 Warmup 200 / 200 Warmup 200 / 200 Sample 10 / 1000 Sample 20 / 1000 Sample 30 / 1000 Sample 40 / 1000 Sample 50 / 1000 Sample 60 / 1000 Sample 70 / 1000 Sample 80 / 1000 Sample 90 / 1000 Sample 100 / 1000 Sample 110 / 1000 Sample 120 / 1000 Sample 130 / 1000 Sample 140 / 1000 Sample 150 / 1000 Sample 160 / 1000 Sample 170 / 1000 Sample 180 / 1000 Sample 190 / 1000 Sample 200 / 1000 Sample 210 / 1000 Sample 220 / 1000 Sample 230 / 1000 Sample 240 / 1000 Sample 250 / 1000 Sample 260 / 1000 Sample 270 / 1000 Sample 280 / 1000 Sample 290 / 1000 Sample 300 / 1000 Sample 310 / 1000 Sample 320 / 1000 Sample 330 / 1000 Sample 340 / 1000 Sample 350 / 1000 Sample 360 / 1000 Sample 370 / 1000 Sample 380 / 1000 Sample 390 / 1000 Sample 400 / 1000 Sample 410 / 1000 Sample 420 / 1000 Sample 430 / 1000 Sample 440 / 1000 Sample 450 / 1000 Sample 460 / 1000 Sample 470 / 1000 Sample 480 / 1000 Sample 490 / 1000 Sample 500 / 1000 Sample 510 / 1000 Sample 520 / 1000 Sample 530 / 1000 Sample 540 / 1000 Sample 550 / 1000 Sample 560 / 1000 Sample 570 / 1000 Sample 580 / 1000 Sample 590 / 1000 Sample 600 / 1000 Sample 610 / 1000 Sample 620 / 1000 Sample 630 / 1000 Sample 640 / 1000 Sample 650 / 1000 Sample 660 / 1000 Sample 670 / 1000 Sample 680 / 1000 Sample 690 / 1000 Sample 700 / 1000 Sample 710 / 1000 Sample 720 / 1000 Sample 730 / 1000 Sample 740 / 1000 Sample 750 / 1000 Sample 760 / 1000 Sample 770 / 1000 Sample 780 / 1000 Sample 790 / 1000 Sample 800 / 1000 Sample 810 / 1000 Sample 820 / 1000 Sample 830 / 1000 Sample 840 / 1000 Sample 850 / 1000 Sample 860 / 1000 Sample 870 / 1000 Sample 880 / 1000 Sample 890 / 1000 Sample 900 / 1000 Sample 910 / 1000 Sample 920 / 1000 Sample 930 / 1000 Sample 940 / 1000 Sample 950 / 1000 Sample 960 / 1000 Sample 970 / 1000 Sample 980 / 1000 Sample 990 / 1000 Sample 1000 / 1000 Sample 1000 / 1000 Elapsed time: 5.021841526031494 Plotting results Warmup: 0%| | 0/200 [00:00, , ] .. GENERATED FROM PYTHON SOURCE LINES 113-114 Conditioning on random variables (example 1) .. GENERATED FROM PYTHON SOURCE LINES 114-123 .. code-block:: Python s = Gaussian(0, 1).rv x = Gaussian(0, s).rv y = Gaussian(0, lambda d: d).rv z = x+y z.condition(s=1) z.condition(d=2) .. rst-class:: sphx-glr-script-out .. code-block:: none Transformed Random Variable Expression: x + y Components: x ~ CUQI Gaussian. Conditioning variables ['s']. y ~ CUQI Gaussian. .. GENERATED FROM PYTHON SOURCE LINES 124-125 Or conditioning on the variables s, or d .. GENERATED FROM PYTHON SOURCE LINES 125-127 .. code-block:: Python z.condition(s=1) .. rst-class:: sphx-glr-script-out .. code-block:: none Transformed Random Variable Expression: x + y Components: x ~ CUQI Gaussian. y ~ CUQI Gaussian. Conditioning variables ['d']. .. GENERATED FROM PYTHON SOURCE LINES 128-129 Conditioning on random variables (example 2) .. GENERATED FROM PYTHON SOURCE LINES 129-149 .. code-block:: Python from cuqi.testproblem import Deconvolution1D from cuqi.distribution import Gaussian, Gamma, GMRF from cuqi.experimental.algebra import RandomVariable from cuqi.problem import BayesianProblem import numpy as np # Forward model A, y_obs, info = Deconvolution1D(dim=4).get_components() # Bayesian Problem (defined using Random Variables) d = Gamma(1, 1e-4).rv s = Gamma(1, 1e-4).rv x = GMRF(np.zeros(A.domain_dim), d).rv y = Gaussian(A @ x, 1/s).rv z = x+y z.condition(x=np.zeros(A.domain_dim)) .. rst-class:: sphx-glr-script-out .. code-block:: none Transformed Random Variable Expression: [0. 0. 0. 0.] + y Components: y ~ CUQI Gaussian. Conditioning variables ['s']. .. GENERATED FROM PYTHON SOURCE LINES 150-154 Sampling from random variables ------------------------------ Random variables can be sampled using the `sample` method. The method returns a sample from the distribution of the random variable. .. GENERATED FROM PYTHON SOURCE LINES 154-159 .. code-block:: Python x = RandomVariable(Normal(0, 1)) print(f"Sample from x: {x.sample()}") .. rst-class:: sphx-glr-script-out .. code-block:: none Sample from x: 1.1917018637858252 .. GENERATED FROM PYTHON SOURCE LINES 160-162 This can be combined with algebraic operations to sample from more complex random variables. .. GENERATED FROM PYTHON SOURCE LINES 162-167 .. code-block:: Python z = x + x**2 + 25 print(f"Sample from z: {z.sample()}") .. rst-class:: sphx-glr-script-out .. code-block:: none Sample from z: 28.449155125413096 .. GENERATED FROM PYTHON SOURCE LINES 168-176 Constructing a Beta distribution using Gamma random variables ------------------------------------------------------------- Random variables can also be combined to create new distributions. This is primarily useful for sampling at this stage. For example, a Beta distribution can be constructed from two Gamma distributions: If X ~ Gamma(a1, 1) and Y ~ Gamma(a2, 1), then Z = X / (X + Y) ~ Beta(a1, a2). We illustrate this by comparing samples from a Beta distribution to samples constructed using two Gamma distributions. .. GENERATED FROM PYTHON SOURCE LINES 176-197 .. code-block:: Python from cuqi.distribution import Beta, Gamma # Define the shape parameters of the Beta distribution a1, a2 = 3, 2 # Step 1: Directly define a Beta distribution z_ref = RandomVariable(Beta(a1, a2)) # Step 2: Construct the Beta distribution using Gamma random variables x = RandomVariable(Gamma(a1, 1)) # X ~ Gamma(a1, 1) y = RandomVariable(Gamma(a2, 1)) # Y ~ Gamma(a2, 1) z = x / (x + y) # Z ~ Beta(a1, a2) # Step 3: Sample from both distributions z_samples = z.sample(10000) # Samples from constructed Beta distribution z_ref_samples = z_ref.sample(10000) # Samples from direct Beta distribution # Step 4: Plot histograms of the samples for comparison z_samples.hist_chain([0], bins=100) z_ref_samples.hist_chain([0], bins=100) .. image-sg:: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_005.png :alt: RandomVariablesAndAlgebra :srcset: /user/_auto_howtos/images/sphx_glr_RandomVariablesAndAlgebra_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.134 seconds) .. _sphx_glr_download_user__auto_howtos_RandomVariablesAndAlgebra.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: RandomVariablesAndAlgebra.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: RandomVariablesAndAlgebra.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: RandomVariablesAndAlgebra.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_