Note
Go to the end to download the full example code.
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.
Defining Random Variables#
# 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
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.
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")
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.
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")
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.
print("Array operations: \n")
print(f"x[0] + y[1] yields:\n{x[0] + y[1]}\n")
Array operations:
x[0] + y[1] yields:
Transformed Random Variable
Expression: x[0] + y[1]
Components:
x ~ CUQI Normal.
y ~ CUQI Normal.
Utilizing the recorded operations#
We can evaluate the recorded operations by calling the random variable object with the desired values for the random variables.
# 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)}")
z=(x + y)^2 evaluated at x=1, y=2 yields: 9
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.
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)
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.594318628311157
Plotting results
Warmup: 0%| | 0/200 [00:00<?, ?it/s]
Warmup: 10%|█ | 20/200 [00:00<00:00, 199.40it/s]
Warmup: 20%|██ | 40/200 [00:00<00:00, 198.44it/s]
Warmup: 30%|███ | 60/200 [00:00<00:00, 198.09it/s]
Warmup: 40%|████ | 81/200 [00:00<00:00, 199.31it/s]
Warmup: 50%|█████ | 101/200 [00:00<00:00, 199.55it/s]
Warmup: 61%|██████ | 122/200 [00:00<00:00, 200.34it/s]
Warmup: 72%|███████▏ | 143/200 [00:00<00:00, 200.85it/s]
Warmup: 82%|████████▏ | 164/200 [00:00<00:00, 201.43it/s]
Warmup: 92%|█████████▎| 185/200 [00:00<00:00, 201.35it/s]
Warmup: 100%|██████████| 200/200 [00:00<00:00, 200.51it/s]
Sample: 0%| | 0/1000 [00:00<?, ?it/s]
Sample: 2%|▏ | 21/1000 [00:00<00:04, 201.78it/s]
Sample: 4%|▍ | 42/1000 [00:00<00:04, 202.72it/s]
Sample: 6%|▋ | 63/1000 [00:00<00:04, 202.78it/s]
Sample: 8%|▊ | 84/1000 [00:00<00:04, 201.19it/s]
Sample: 10%|█ | 105/1000 [00:00<00:04, 201.81it/s]
Sample: 13%|█▎ | 126/1000 [00:00<00:04, 201.35it/s]
Sample: 15%|█▍ | 147/1000 [00:00<00:04, 201.48it/s]
Sample: 17%|█▋ | 168/1000 [00:00<00:04, 201.98it/s]
Sample: 19%|█▉ | 189/1000 [00:00<00:04, 202.13it/s]
Sample: 21%|██ | 210/1000 [00:01<00:03, 202.37it/s]
Sample: 23%|██▎ | 231/1000 [00:01<00:03, 202.27it/s]
Sample: 25%|██▌ | 252/1000 [00:01<00:03, 200.34it/s]
Sample: 27%|██▋ | 273/1000 [00:01<00:03, 200.75it/s]
Sample: 29%|██▉ | 294/1000 [00:01<00:03, 201.25it/s]
Sample: 32%|███▏ | 315/1000 [00:01<00:03, 200.66it/s]
Sample: 34%|███▎ | 336/1000 [00:01<00:03, 201.14it/s]
Sample: 36%|███▌ | 357/1000 [00:01<00:03, 201.25it/s]
Sample: 38%|███▊ | 378/1000 [00:01<00:03, 201.58it/s]
Sample: 40%|███▉ | 399/1000 [00:01<00:02, 201.61it/s]
Sample: 42%|████▏ | 420/1000 [00:02<00:02, 201.65it/s]
Sample: 44%|████▍ | 441/1000 [00:02<00:02, 201.24it/s]
Sample: 46%|████▌ | 462/1000 [00:02<00:02, 201.22it/s]
Sample: 48%|████▊ | 483/1000 [00:02<00:02, 201.57it/s]
Sample: 50%|█████ | 504/1000 [00:02<00:02, 201.67it/s]
Sample: 52%|█████▎ | 525/1000 [00:02<00:02, 201.63it/s]
Sample: 55%|█████▍ | 546/1000 [00:02<00:02, 200.64it/s]
Sample: 57%|█████▋ | 567/1000 [00:02<00:02, 200.46it/s]
Sample: 59%|█████▉ | 588/1000 [00:02<00:02, 200.58it/s]
Sample: 61%|██████ | 609/1000 [00:03<00:01, 201.53it/s]
Sample: 63%|██████▎ | 630/1000 [00:03<00:01, 201.50it/s]
Sample: 65%|██████▌ | 651/1000 [00:03<00:01, 199.99it/s]
Sample: 67%|██████▋ | 672/1000 [00:03<00:01, 201.02it/s]
Sample: 69%|██████▉ | 693/1000 [00:03<00:01, 201.67it/s]
Sample: 71%|███████▏ | 714/1000 [00:03<00:01, 202.33it/s]
Sample: 74%|███████▎ | 735/1000 [00:03<00:01, 202.57it/s]
Sample: 76%|███████▌ | 756/1000 [00:03<00:01, 202.39it/s]
Sample: 78%|███████▊ | 777/1000 [00:03<00:01, 201.51it/s]
Sample: 80%|███████▉ | 798/1000 [00:03<00:01, 201.15it/s]
Sample: 82%|████████▏ | 819/1000 [00:04<00:00, 201.01it/s]
Sample: 84%|████████▍ | 840/1000 [00:04<00:00, 200.51it/s]
Sample: 86%|████████▌ | 861/1000 [00:04<00:00, 200.02it/s]
Sample: 88%|████████▊ | 882/1000 [00:04<00:00, 200.17it/s]
Sample: 90%|█████████ | 903/1000 [00:04<00:00, 200.68it/s]
Sample: 92%|█████████▏| 924/1000 [00:04<00:00, 200.96it/s]
Sample: 94%|█████████▍| 945/1000 [00:04<00:00, 201.24it/s]
Sample: 97%|█████████▋| 966/1000 [00:04<00:00, 201.04it/s]
Sample: 99%|█████████▊| 987/1000 [00:04<00:00, 201.17it/s]
Sample: 100%|██████████| 1000/1000 [00:04<00:00, 201.32it/s]
[<matplotlib.lines.Line2D object at 0x7f2f13e7bf20>, <matplotlib.lines.Line2D object at 0x7f2f13e78920>, <matplotlib.collections.FillBetweenPolyCollection object at 0x7f2f137fd070>]
Conditioning on random variables (example 1)
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)
Transformed Random Variable
Expression: x + y
Components:
x ~ CUQI Gaussian. Conditioning variables ['s'].
y ~ CUQI Gaussian.
Or conditioning on the variables s, or d
z.condition(s=1)
Transformed Random Variable
Expression: x + y
Components:
x ~ CUQI Gaussian.
y ~ CUQI Gaussian. Conditioning variables ['d'].
Conditioning on random variables (example 2)
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))
Transformed Random Variable
Expression: [0. 0. 0. 0.] + y
Components:
y ~ CUQI Gaussian. Conditioning variables ['s'].
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.
x = RandomVariable(Normal(0, 1))
print(f"Sample from x: {x.sample()}")
Sample from x: 1.1917018637858252
This can be combined with algebraic operations to sample from more complex random variables.
z = x + x**2 + 25
print(f"Sample from z: {z.sample()}")
Sample from z: 28.449155125413096
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.
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)
<BarContainer object of 100 artists>
Total running time of the script: (0 minutes 13.774 seconds)