.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user/_auto_tutorials/uq_in_deconvolution_1d.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_tutorials_uq_in_deconvolution_1d.py: Uncertainty Quantification in one-dimensional deconvolution =========================================================== This tutorial walks through the process of solving a simple 1D deconvolution problem in a Bayesian setting. It also shows how to define such a convolution model in CUQIpy. .. GENERATED FROM PYTHON SOURCE LINES 11-14 Setup ----- We start by importing the necessary modules .. GENERATED FROM PYTHON SOURCE LINES 14-19 .. code-block:: Python import cuqi import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 20-37 Setting up the forward model ---------------------------- We start by defining the forward model. In this case, we will use a simple convolution model. The forward model is defined by the following equation: .. math:: \mathbf{y} = \mathbf{A} \mathbf{x} where :math:`\mathbf{y}` is the data, :math:`\mathbf{A}` is the convolution (forward model) operator, and :math:`\mathbf{x}` is the solution. The easiest way to define the forward model is to use the testproblem module. This module contains a number of pre-defined test problems that contain the forward model and synthetic data. In this case, we will use the :class:`cuqi.testproblem.Deconvolution1D` test problem. We extract the forward model and synthetic data from the test problem by calling the :func:`get_components` method. .. GENERATED FROM PYTHON SOURCE LINES 37-41 .. code-block:: Python # Forward model and data A, y_data, info = cuqi.testproblem.Deconvolution1D().get_components() .. GENERATED FROM PYTHON SOURCE LINES 42-49 There are many parameters that can be set when creating the test problem. For more details see the :class:`cuqi.testproblem.Deconvolution1D` documentation. In this case, we will use the default parameters. The :func:`get_components` method returns the forward model, synthetic data, and a :class:`~ProblemInfo` object that contains information about the test problem. Let's take a look at the forward model .. GENERATED FROM PYTHON SOURCE LINES 49-52 .. code-block:: Python print(A) .. rst-class:: sphx-glr-script-out .. code-block:: none CUQI LinearModel: Continuous1D[128] -> Continuous1D[128]. Forward parameters: ['x']. .. GENERATED FROM PYTHON SOURCE LINES 53-57 We see that the forward model is a a :class:`~cuqi.model.LinearModel` object. This object contains the forward model and the adjoint model. We also see that the domain and range of the forward model are both continuous 1D spaces. Finally, we see that the default forward parameters are set to :math:`\mathbf{x}`. .. GENERATED FROM PYTHON SOURCE LINES 59-61 Let's take a look at the synthetic data and compare with the exact solution that we can find in the :class:`~ProblemInfo` object. .. GENERATED FROM PYTHON SOURCE LINES 61-67 .. code-block:: Python y_data.plot(label="Synthetic data") info.exactSolution.plot(label="Exact solution") plt.title("Deconvolution 1D problem") plt.legend() .. image-sg:: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_001.png :alt: Deconvolution 1D problem :srcset: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 68-74 Setting up the prior -------------------- We now need to define the prior distribution for the solution. In this case, we will use a Gaussian Markov Random Field (GMRF) prior. For more details on the GMRF prior, see the :class:`cuqi.distribution.GMRF` documentation. .. GENERATED FROM PYTHON SOURCE LINES 74-78 .. code-block:: Python x = cuqi.distribution.GMRF(np.zeros(A.domain_dim), 200) .. GENERATED FROM PYTHON SOURCE LINES 79-84 Setting up the likelihood ------------------------- We now need to define the likelihood. First let us take a look at the information provided by the test problem. .. GENERATED FROM PYTHON SOURCE LINES 84-87 .. code-block:: Python print(info.infoString) .. rst-class:: sphx-glr-script-out .. code-block:: none Noise type: Additive Gaussian with std: 0.01 .. GENERATED FROM PYTHON SOURCE LINES 88-91 We see that the noise level is known and that the noise is Gaussian. We can use this information to define the likelihood. In this case, we will use a :class:`~cuqi.distribution.Gaussian` distribution. .. GENERATED FROM PYTHON SOURCE LINES 91-94 .. code-block:: Python y = cuqi.distribution.Gaussian(A @ x, 0.01**2) .. GENERATED FROM PYTHON SOURCE LINES 95-101 Bayesian problem (Joint distribution) ------------------------------------- After defining the prior and likelihood, we can now define the Bayesian problem. The Bayesian problem is defined by the joint distribution of the solution and the data. This can be seen when we print the Bayesian problem. .. GENERATED FROM PYTHON SOURCE LINES 101-106 .. code-block:: Python BP = cuqi.problem.BayesianProblem(y, x) print(BP) .. rst-class:: sphx-glr-script-out .. code-block:: none BayesianProblem with target: JointDistribution( Equation: p(y,x) = p(y|x)p(x) Densities: y ~ CUQI Gaussian. Conditioning variables ['x']. x ~ CUQI GMRF. ) .. GENERATED FROM PYTHON SOURCE LINES 107-111 Setting the data (posterior) ---------------------------- Now to set the data, we need to call the :func:`~cuqi.problem.BayesianProblem.set_data` .. GENERATED FROM PYTHON SOURCE LINES 111-116 .. code-block:: Python BP.set_data(y=y_data) print(BP) .. rst-class:: sphx-glr-script-out .. code-block:: none BayesianProblem with target: Posterior( Equation: p(x|y) ∝ L(x|y)p(x) Densities: y ~ CUQI Gaussian Likelihood function. Parameters ['x']. x ~ CUQI GMRF. ) .. GENERATED FROM PYTHON SOURCE LINES 117-121 Sampling from the posterior --------------------------- We can then use the automatic sampling method to sample from the posterior distribution. .. GENERATED FROM PYTHON SOURCE LINES 121-124 .. code-block:: Python samples = BP.sample_posterior(1000) .. rst-class:: sphx-glr-script-out .. code-block:: none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Automatic sampler selection is a work-in-progress. !!! !!! Always validate the computed results. !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Using cuqi.sampler LinearRTO sampler. burn-in: 20% Warmup: 0%| | 0/200 [00:00, , ] .. GENERATED FROM PYTHON SOURCE LINES 131-136 Unknown noise level ------------------- In the previous example, we assumed that we knew the noise level of the data. In many cases, this is not the case. If we do not know the noise level, we can use a :class:`~cuqi.distribution.Gamma` distribution to model the noise level. .. GENERATED FROM PYTHON SOURCE LINES 136-139 .. code-block:: Python s = cuqi.distribution.Gamma(1, 1e-4) .. GENERATED FROM PYTHON SOURCE LINES 140-142 Update likelihood with unknown noise level ------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 142-145 .. code-block:: Python y = cuqi.distribution.Gaussian(A @ x, prec=lambda s: s) .. GENERATED FROM PYTHON SOURCE LINES 146-148 Bayesian problem (Joint distribution) ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 148-153 .. code-block:: Python BP = cuqi.problem.BayesianProblem(y, x, s) print(BP) .. rst-class:: sphx-glr-script-out .. code-block:: none BayesianProblem with target: JointDistribution( Equation: p(y,x,s) = p(y|x,s)p(x)p(s) Densities: y ~ CUQI Gaussian. Conditioning variables ['x', 's']. x ~ CUQI GMRF. s ~ CUQI Gamma. ) .. GENERATED FROM PYTHON SOURCE LINES 154-157 Setting the data (posterior) ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 157-162 .. code-block:: Python BP.set_data(y=y_data) print(BP) .. rst-class:: sphx-glr-script-out .. code-block:: none BayesianProblem with target: JointDistribution( Equation: p(x,s|y) ∝ L(x,s|y)p(x)p(s) Densities: y ~ CUQI Gaussian Likelihood function. Parameters ['x', 's']. x ~ CUQI GMRF. s ~ CUQI Gamma. ) .. GENERATED FROM PYTHON SOURCE LINES 163-165 Sampling from the posterior --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 165-169 .. code-block:: Python samples = BP.sample_posterior(1000) .. rst-class:: sphx-glr-script-out .. code-block:: none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! 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) Warmup: 0%| | 0/200 [00:00, ]], dtype=object) .. GENERATED FROM PYTHON SOURCE LINES 179-181 We see that the estimated noise level is close to the true noise level. Let's now look at the estimated solution .. GENERATED FROM PYTHON SOURCE LINES 181-186 .. code-block:: Python samples["x"].plot_ci(exact=info.exactSolution) .. image-sg:: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_004.png :alt: uq in deconvolution 1d :srcset: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [, , ] .. GENERATED FROM PYTHON SOURCE LINES 187-188 We can even plot traces of "x" for a few cases and compare .. GENERATED FROM PYTHON SOURCE LINES 188-190 .. code-block:: Python samples["x"].plot_trace(exact=info.exactSolution) .. image-sg:: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_005.png :alt: x6, x6, x73, x73, x80, x80, x81, x81, x105, x105 :srcset: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Selecting 5 randomly chosen variables array([[, ], [, ], [, ], [, ], [, ]], dtype=object) .. GENERATED FROM PYTHON SOURCE LINES 191-192 And finally we note that the UQ method does this analysis automatically and shows a selected number of plots .. GENERATED FROM PYTHON SOURCE LINES 192-194 .. code-block:: Python BP.UQ(exact={"x": info.exactSolution, "s": 1/0.01**2}) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_006.png :alt: uq in deconvolution 1d :srcset: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_006.png :class: sphx-glr-multi-img * .. image-sg:: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_007.png :alt: s, s :srcset: /user/_auto_tutorials/images/sphx_glr_uq_in_deconvolution_1d_007.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) Warmup: 0%| | 0/200 [00:00` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: uq_in_deconvolution_1d.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: uq_in_deconvolution_1d.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_