.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user/_auto_howtos/random_variables_and_algebra.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_random_variables_and_algebra.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.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.algebra import RandomVariable from cuqi.problem import BayesianProblem from cuqi.distribution import JointDistribution from cuqi.sampler 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_random_variables_and_algebra_001.png :alt: random variables and algebra :srcset: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_001.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_002.png :alt: s, s :srcset: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_002.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_003.png :alt: d, d :srcset: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_003.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_004.png :alt: random variables and algebra :srcset: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_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 cuqi.sampler HybridGibbs sampler burn-in: 20% Automatically determined sampling strategy: x: LinearRTO (mcmc.sampler) s: Conjugate (mcmc.sampler) d: Conjugate (mcmc.sampler) 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.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.5313150719729232 .. 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: 25.609017298194665 .. 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_random_variables_and_algebra_005.png :alt: random variables and algebra :srcset: /user/_auto_howtos/images/sphx_glr_random_variables_and_algebra_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 17.850 seconds) .. _sphx_glr_download_user__auto_howtos_random_variables_and_algebra.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: random_variables_and_algebra.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: random_variables_and_algebra.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: random_variables_and_algebra.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_