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.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)
Computing 1000 samples
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Automatic sampler selection is a work-in-progress. !!!
!!! Always validate the computed results. !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Using samplers from cuqi.sampler !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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<?, ?it/s]
Warmup: 7%|▋ | 14/200 [00:00<00:01, 132.41it/s]
Warmup: 14%|█▍ | 28/200 [00:00<00:01, 132.46it/s]
Warmup: 21%|██ | 42/200 [00:00<00:01, 133.09it/s]
Warmup: 28%|██▊ | 56/200 [00:00<00:01, 133.06it/s]
Warmup: 35%|███▌ | 70/200 [00:00<00:00, 133.17it/s]
Warmup: 42%|████▏ | 84/200 [00:00<00:00, 133.19it/s]
Warmup: 49%|████▉ | 98/200 [00:00<00:00, 133.59it/s]
Warmup: 56%|█████▌ | 112/200 [00:00<00:00, 133.69it/s]
Warmup: 63%|██████▎ | 126/200 [00:00<00:00, 133.82it/s]
Warmup: 70%|███████ | 140/200 [00:01<00:00, 133.85it/s]
Warmup: 77%|███████▋ | 154/200 [00:01<00:00, 133.70it/s]
Warmup: 84%|████████▍ | 168/200 [00:01<00:00, 133.72it/s]
Warmup: 91%|█████████ | 182/200 [00:01<00:00, 133.76it/s]
Warmup: 98%|█████████▊| 196/200 [00:01<00:00, 133.87it/s]
Warmup: 100%|██████████| 200/200 [00:01<00:00, 133.53it/s]
Sample: 0%| | 0/1000 [00:00<?, ?it/s]
Sample: 1%|▏ | 14/1000 [00:00<00:07, 134.01it/s]
Sample: 3%|▎ | 28/1000 [00:00<00:07, 134.02it/s]
Sample: 4%|▍ | 42/1000 [00:00<00:07, 133.43it/s]
Sample: 6%|▌ | 56/1000 [00:00<00:07, 133.71it/s]
Sample: 7%|▋ | 70/1000 [00:00<00:06, 133.97it/s]
Sample: 8%|▊ | 84/1000 [00:00<00:06, 133.73it/s]
Sample: 10%|▉ | 98/1000 [00:00<00:06, 133.93it/s]
Sample: 11%|█ | 112/1000 [00:00<00:06, 134.11it/s]
Sample: 13%|█▎ | 126/1000 [00:00<00:06, 133.89it/s]
Sample: 14%|█▍ | 140/1000 [00:01<00:06, 133.77it/s]
Sample: 15%|█▌ | 154/1000 [00:01<00:06, 133.89it/s]
Sample: 17%|█▋ | 168/1000 [00:01<00:06, 134.11it/s]
Sample: 18%|█▊ | 182/1000 [00:01<00:06, 133.92it/s]
Sample: 20%|█▉ | 196/1000 [00:01<00:05, 134.01it/s]
Sample: 21%|██ | 210/1000 [00:01<00:05, 134.14it/s]
Sample: 22%|██▏ | 224/1000 [00:01<00:05, 133.96it/s]
Sample: 24%|██▍ | 238/1000 [00:01<00:05, 133.99it/s]
Sample: 25%|██▌ | 252/1000 [00:01<00:05, 134.14it/s]
Sample: 27%|██▋ | 266/1000 [00:01<00:05, 134.20it/s]
Sample: 28%|██▊ | 280/1000 [00:02<00:05, 134.14it/s]
Sample: 29%|██▉ | 294/1000 [00:02<00:05, 134.18it/s]
Sample: 31%|███ | 308/1000 [00:02<00:05, 134.14it/s]
Sample: 32%|███▏ | 322/1000 [00:02<00:05, 134.42it/s]
Sample: 34%|███▎ | 336/1000 [00:02<00:04, 134.60it/s]
Sample: 35%|███▌ | 350/1000 [00:02<00:04, 134.49it/s]
Sample: 36%|███▋ | 364/1000 [00:02<00:04, 134.65it/s]
Sample: 38%|███▊ | 378/1000 [00:02<00:04, 134.57it/s]
Sample: 39%|███▉ | 392/1000 [00:02<00:04, 134.38it/s]
Sample: 41%|████ | 406/1000 [00:03<00:04, 134.32it/s]
Sample: 42%|████▏ | 420/1000 [00:03<00:04, 134.39it/s]
Sample: 43%|████▎ | 434/1000 [00:03<00:04, 134.60it/s]
Sample: 45%|████▍ | 448/1000 [00:03<00:04, 134.67it/s]
Sample: 46%|████▌ | 462/1000 [00:03<00:03, 134.71it/s]
Sample: 48%|████▊ | 476/1000 [00:03<00:03, 134.87it/s]
Sample: 49%|████▉ | 490/1000 [00:03<00:03, 134.75it/s]
Sample: 50%|█████ | 504/1000 [00:03<00:03, 134.53it/s]
Sample: 52%|█████▏ | 518/1000 [00:03<00:03, 134.70it/s]
Sample: 53%|█████▎ | 532/1000 [00:03<00:03, 134.85it/s]
Sample: 55%|█████▍ | 546/1000 [00:04<00:03, 134.89it/s]
Sample: 56%|█████▌ | 560/1000 [00:04<00:03, 134.66it/s]
Sample: 57%|█████▋ | 574/1000 [00:04<00:03, 134.65it/s]
Sample: 59%|█████▉ | 588/1000 [00:04<00:03, 134.67it/s]
Sample: 60%|██████ | 602/1000 [00:04<00:02, 134.77it/s]
Sample: 62%|██████▏ | 616/1000 [00:04<00:02, 134.92it/s]
Sample: 63%|██████▎ | 630/1000 [00:04<00:02, 134.81it/s]
Sample: 64%|██████▍ | 644/1000 [00:04<00:02, 134.85it/s]
Sample: 66%|██████▌ | 658/1000 [00:04<00:02, 134.87it/s]
Sample: 67%|██████▋ | 672/1000 [00:05<00:02, 134.77it/s]
Sample: 69%|██████▊ | 686/1000 [00:05<00:02, 134.74it/s]
Sample: 70%|███████ | 700/1000 [00:05<00:02, 134.81it/s]
Sample: 71%|███████▏ | 714/1000 [00:05<00:02, 134.83it/s]
Sample: 73%|███████▎ | 728/1000 [00:05<00:02, 134.85it/s]
Sample: 74%|███████▍ | 742/1000 [00:05<00:01, 134.75it/s]
Sample: 76%|███████▌ | 756/1000 [00:05<00:01, 134.53it/s]
Sample: 77%|███████▋ | 770/1000 [00:05<00:01, 134.63it/s]
Sample: 78%|███████▊ | 784/1000 [00:05<00:01, 134.62it/s]
Sample: 80%|███████▉ | 798/1000 [00:05<00:01, 134.71it/s]
Sample: 81%|████████ | 812/1000 [00:06<00:01, 134.52it/s]
Sample: 83%|████████▎ | 826/1000 [00:06<00:01, 134.44it/s]
Sample: 84%|████████▍ | 840/1000 [00:06<00:01, 134.32it/s]
Sample: 85%|████████▌ | 854/1000 [00:06<00:01, 134.30it/s]
Sample: 87%|████████▋ | 868/1000 [00:06<00:00, 134.33it/s]
Sample: 88%|████████▊ | 882/1000 [00:06<00:00, 134.39it/s]
Sample: 90%|████████▉ | 896/1000 [00:06<00:00, 134.34it/s]
Sample: 91%|█████████ | 910/1000 [00:06<00:00, 134.26it/s]
Sample: 92%|█████████▏| 924/1000 [00:06<00:00, 134.30it/s]
Sample: 94%|█████████▍| 938/1000 [00:06<00:00, 133.64it/s]
Sample: 95%|█████████▌| 952/1000 [00:07<00:00, 133.61it/s]
Sample: 97%|█████████▋| 966/1000 [00:07<00:00, 133.76it/s]
Sample: 98%|█████████▊| 980/1000 [00:07<00:00, 133.77it/s]
Sample: 99%|█████████▉| 994/1000 [00:07<00:00, 133.73it/s]
Sample: 100%|██████████| 1000/1000 [00:07<00:00, 134.34it/s]
Elapsed time: 8.950470447540283
Plotting results
Warmup: 0%| | 0/200 [00:00<?, ?it/s]
Warmup: 7%|▋ | 14/200 [00:00<00:01, 132.05it/s]
Warmup: 14%|█▍ | 28/200 [00:00<00:01, 132.78it/s]
Warmup: 21%|██ | 42/200 [00:00<00:01, 133.05it/s]
Warmup: 28%|██▊ | 56/200 [00:00<00:01, 133.43it/s]
Warmup: 35%|███▌ | 70/200 [00:00<00:00, 133.51it/s]
Warmup: 42%|████▏ | 84/200 [00:00<00:00, 133.42it/s]
Warmup: 49%|████▉ | 98/200 [00:00<00:00, 133.47it/s]
Warmup: 56%|█████▌ | 112/200 [00:00<00:00, 133.50it/s]
Warmup: 63%|██████▎ | 126/200 [00:00<00:00, 133.10it/s]
Warmup: 70%|███████ | 140/200 [00:01<00:00, 133.26it/s]
Warmup: 77%|███████▋ | 154/200 [00:01<00:00, 133.22it/s]
Warmup: 84%|████████▍ | 168/200 [00:01<00:00, 133.33it/s]
Warmup: 91%|█████████ | 182/200 [00:01<00:00, 133.52it/s]
Warmup: 98%|█████████▊| 196/200 [00:01<00:00, 133.68it/s]
Warmup: 100%|██████████| 200/200 [00:01<00:00, 133.37it/s]
Sample: 0%| | 0/1000 [00:00<?, ?it/s]
Sample: 1%|▏ | 14/1000 [00:00<00:07, 133.89it/s]
Sample: 3%|▎ | 28/1000 [00:00<00:07, 133.79it/s]
Sample: 4%|▍ | 42/1000 [00:00<00:07, 134.11it/s]
Sample: 6%|▌ | 56/1000 [00:00<00:07, 134.35it/s]
Sample: 7%|▋ | 70/1000 [00:00<00:06, 134.43it/s]
Sample: 8%|▊ | 84/1000 [00:00<00:06, 134.13it/s]
Sample: 10%|▉ | 98/1000 [00:00<00:06, 134.14it/s]
Sample: 11%|█ | 112/1000 [00:00<00:06, 134.04it/s]
Sample: 13%|█▎ | 126/1000 [00:00<00:06, 133.94it/s]
Sample: 14%|█▍ | 140/1000 [00:01<00:06, 133.94it/s]
Sample: 15%|█▌ | 154/1000 [00:01<00:06, 133.72it/s]
Sample: 17%|█▋ | 168/1000 [00:01<00:06, 133.88it/s]
Sample: 18%|█▊ | 182/1000 [00:01<00:06, 134.02it/s]
Sample: 20%|█▉ | 196/1000 [00:01<00:05, 134.02it/s]
Sample: 21%|██ | 210/1000 [00:01<00:05, 133.86it/s]
Sample: 22%|██▏ | 224/1000 [00:01<00:05, 133.94it/s]
Sample: 24%|██▍ | 238/1000 [00:01<00:05, 134.13it/s]
Sample: 25%|██▌ | 252/1000 [00:01<00:05, 134.11it/s]
Sample: 27%|██▋ | 266/1000 [00:01<00:05, 134.05it/s]
Sample: 28%|██▊ | 280/1000 [00:02<00:05, 134.14it/s]
Sample: 29%|██▉ | 294/1000 [00:02<00:05, 134.22it/s]
Sample: 31%|███ | 308/1000 [00:02<00:05, 134.20it/s]
Sample: 32%|███▏ | 322/1000 [00:02<00:05, 134.05it/s]
Sample: 34%|███▎ | 336/1000 [00:02<00:04, 134.13it/s]
Sample: 35%|███▌ | 350/1000 [00:02<00:04, 133.96it/s]
Sample: 36%|███▋ | 364/1000 [00:02<00:04, 134.09it/s]
Sample: 38%|███▊ | 378/1000 [00:02<00:04, 134.17it/s]
Sample: 39%|███▉ | 392/1000 [00:02<00:04, 134.32it/s]
Sample: 41%|████ | 406/1000 [00:03<00:04, 134.37it/s]
Sample: 42%|████▏ | 420/1000 [00:03<00:04, 134.46it/s]
Sample: 43%|████▎ | 434/1000 [00:03<00:04, 134.26it/s]
Sample: 45%|████▍ | 448/1000 [00:03<00:04, 134.20it/s]
Sample: 46%|████▌ | 462/1000 [00:03<00:04, 134.03it/s]
Sample: 48%|████▊ | 476/1000 [00:03<00:03, 133.88it/s]
Sample: 49%|████▉ | 490/1000 [00:03<00:03, 133.61it/s]
Sample: 50%|█████ | 504/1000 [00:03<00:03, 133.88it/s]
Sample: 52%|█████▏ | 518/1000 [00:03<00:03, 133.79it/s]
Sample: 53%|█████▎ | 532/1000 [00:03<00:03, 133.69it/s]
Sample: 55%|█████▍ | 546/1000 [00:04<00:03, 133.50it/s]
Sample: 56%|█████▌ | 560/1000 [00:04<00:03, 133.72it/s]
Sample: 57%|█████▋ | 574/1000 [00:04<00:03, 133.69it/s]
Sample: 59%|█████▉ | 588/1000 [00:04<00:03, 133.75it/s]
Sample: 60%|██████ | 602/1000 [00:04<00:02, 133.62it/s]
Sample: 62%|██████▏ | 616/1000 [00:04<00:02, 133.52it/s]
Sample: 63%|██████▎ | 630/1000 [00:04<00:02, 133.40it/s]
Sample: 64%|██████▍ | 644/1000 [00:04<00:02, 133.48it/s]
Sample: 66%|██████▌ | 658/1000 [00:04<00:02, 133.75it/s]
Sample: 67%|██████▋ | 672/1000 [00:05<00:02, 133.73it/s]
Sample: 69%|██████▊ | 686/1000 [00:05<00:02, 133.73it/s]
Sample: 70%|███████ | 700/1000 [00:05<00:02, 132.80it/s]
Sample: 71%|███████▏ | 714/1000 [00:05<00:02, 132.72it/s]
Sample: 73%|███████▎ | 728/1000 [00:05<00:02, 133.26it/s]
Sample: 74%|███████▍ | 742/1000 [00:05<00:01, 132.97it/s]
Sample: 76%|███████▌ | 756/1000 [00:05<00:01, 133.33it/s]
Sample: 77%|███████▋ | 770/1000 [00:05<00:01, 133.69it/s]
Sample: 78%|███████▊ | 784/1000 [00:05<00:01, 133.92it/s]
Sample: 80%|███████▉ | 798/1000 [00:05<00:01, 134.13it/s]
Sample: 81%|████████ | 812/1000 [00:06<00:01, 134.12it/s]
Sample: 83%|████████▎ | 826/1000 [00:06<00:01, 134.14it/s]
Sample: 84%|████████▍ | 840/1000 [00:06<00:01, 134.20it/s]
Sample: 85%|████████▌ | 854/1000 [00:06<00:01, 134.33it/s]
Sample: 87%|████████▋ | 868/1000 [00:06<00:00, 134.31it/s]
Sample: 88%|████████▊ | 882/1000 [00:06<00:00, 134.22it/s]
Sample: 90%|████████▉ | 896/1000 [00:06<00:00, 134.44it/s]
Sample: 91%|█████████ | 910/1000 [00:06<00:00, 134.45it/s]
Sample: 92%|█████████▏| 924/1000 [00:06<00:00, 134.47it/s]
Sample: 94%|█████████▍| 938/1000 [00:07<00:00, 134.34it/s]
Sample: 95%|█████████▌| 952/1000 [00:07<00:00, 134.45it/s]
Sample: 97%|█████████▋| 966/1000 [00:07<00:00, 134.30it/s]
Sample: 98%|█████████▊| 980/1000 [00:07<00:00, 134.18it/s]
Sample: 99%|█████████▉| 994/1000 [00:07<00:00, 134.12it/s]
Sample: 100%|██████████| 1000/1000 [00:07<00:00, 133.95it/s]
[<matplotlib.lines.Line2D object at 0x7f58a63f75f0>, <matplotlib.lines.Line2D object at 0x7f58a63f5b80>, <matplotlib.collections.FillBetweenPolyCollection object at 0x7f58a2f93a10>]
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: 0.07100675251932247
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: 25.07235042596335
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 19.801 seconds)