.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user/_auto_howtos/multiple_likelihoods.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_multiple_likelihoods.py: Setting a Bayesian model with multiple likelihoods ================================================== In this example we build a PDE-based Bayesian inverse problem where the Bayesian model has multiple likelihood functions (two different likelihood functions in this case, but it can be readily extended to more functions) for the same Bayesian parameter `theta`, which represents the conductivity parameters in a 1D Poisson problem. Each likelihood is associated with a different model set up. The models we use here are obtained from the test problem :class:`cuqi.testproblem.Poisson1D`. See the class :class:`cuqi.testproblem.Poisson1D` documentation for more details about the forward model. .. GENERATED FROM PYTHON SOURCE LINES 16-17 First we import the python libraries needed. .. GENERATED FROM PYTHON SOURCE LINES 17-22 .. code-block:: Python import cuqi import numpy as np import matplotlib.pyplot as plt from math import ceil .. GENERATED FROM PYTHON SOURCE LINES 23-31 Choose one of the two cases we study in this demo ------------------------------------------------- We can choose between two cases: Choose `set_up = set_ups[0]` for the case where we have two 1D Poisson models that differ in the observation operator only. And choose `set_up = set_ups[1]` for the case where we have two 1D Poisson models that differ in the source term only. Here we demonstrate the first case, `set_up = set_ups[0]`. .. GENERATED FROM PYTHON SOURCE LINES 31-38 .. code-block:: Python set_ups = ["multi_observation", "multi_source"] set_up = set_ups[0] assert set_up == "multi_observation" or set_up == "multi_source",\ "set_up must be either 'multi_observation' or 'multi_source'" .. GENERATED FROM PYTHON SOURCE LINES 39-41 Set up the parameters used in both models ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 41-56 .. code-block:: Python dim = 50 # Number of the model grid points endpoint = 1 # The model domain is the interval [0, endpoint] field_type = "Step" # The conductivity (or diffusivity) field type. # We choose step function parameterization here. SNR = 400 # Signal-to-noise ratio n_steps = 2 # Number of steps in the conductivity (or diffusivity) step field magnitude = 100 # Magnitude of the source term in the Poisson problem # Exact solution x_exact = np.empty(dim) x_exact[:ceil(dim/2)] = 2 x_exact[ceil(dim/2):] = 3 .. GENERATED FROM PYTHON SOURCE LINES 57-64 Set up the first model ---------------------- We set up the first forward model to have observations at the first half of the domain (or observation everywhere if `set_up = set_ups[1]`). We then plot the true conductivity field (the exact solution), the exact data and the noisy data. .. GENERATED FROM PYTHON SOURCE LINES 64-90 .. code-block:: Python observation_grid_map1 = None if set_up == "multi_observation": # Observe on the first half of the domain observation_grid_map1 = lambda x: x[np.where(x<.5)] # The source term signal source1 = lambda xs: magnitude*np.sin(xs*2*np.pi/endpoint)+magnitude # Obtain the forward model from the test problem model1, data1, problemInfo1 = cuqi.testproblem.Poisson1D(dim=dim, endpoint=endpoint, field_type=field_type, field_params={"n_steps": n_steps}, observation_grid_map=observation_grid_map1, exactSolution=x_exact, source=source1, SNR=SNR).get_components() # Plot data, exact data and exact solution plt.figure() data1.plot(label='data') problemInfo1.exactData.plot(label='exact data') problemInfo1.exactSolution.plot(label='exact solution') plt.legend() .. image-sg:: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_001.png :alt: multiple likelihoods :srcset: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 91-98 Set up the second model ----------------------- We set up the second forward model to have observations at the second half of the domain (or observation everywhere and and different source term if `set_up = set_ups[1]`). We then plot the true conductivity field (the exact solution), the exact data and the noisy data. .. GENERATED FROM PYTHON SOURCE LINES 98-127 .. code-block:: Python observation_grid_map2 = None if set_up == "multi_observation": # Observe on the second half of the domain observation_grid_map2 = lambda x: x[np.where(x>=.5)] # The source term signal if set_up == "multi_source": source2 = lambda xs: magnitude*np.sin(2*xs*2*np.pi/endpoint)+magnitude else: source2 = source1 # Obtain the forward model from the test problem model2, data2, problemInfo2 = cuqi.testproblem.Poisson1D(dim=dim, endpoint=endpoint, field_type=field_type, field_params={"n_steps": n_steps}, observation_grid_map=observation_grid_map2, exactSolution=x_exact, source=source2, SNR=SNR).get_components() # Plot data, exact data and exact solution plt.figure() data2.plot(label='data2') problemInfo2.exactData.plot(label='exact data2') problemInfo2.exactSolution.plot(label='exact solution2') plt.legend() .. image-sg:: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_002.png :alt: multiple likelihoods :srcset: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 128-134 Create the prior ---------------- Create the prior for the Bayesian parameter `theta`, which is the expansion coefficients of the conductivity (or diffusivity) step function. We use a Gaussian prior. .. GENERATED FROM PYTHON SOURCE LINES 134-142 .. code-block:: Python theta = cuqi.distribution.Gaussian( mean=np.zeros(model1.domain_dim), cov=3, geometry=model1.domain_geometry, ) .. GENERATED FROM PYTHON SOURCE LINES 143-145 Create the data distributions using the two forward models ---------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 145-163 .. code-block:: Python # Estimate the data noise standard deviation sigma_noise1 = np.linalg.norm(problemInfo1.exactData)/SNR sigma_noise2 = np.linalg.norm(problemInfo2.exactData)/SNR # Create the data distributions y1 = cuqi.distribution.Gaussian( mean=model1(theta), cov=sigma_noise1**2, geometry=model1.range_geometry, ) y2 = cuqi.distribution.Gaussian( mean=model2(theta), cov=sigma_noise2**2, geometry=model2.range_geometry, ) .. GENERATED FROM PYTHON SOURCE LINES 164-167 Formulate the Bayesian inverse problem using the first data distribution (single likelihood) ---------------------------------------------------------------------------------------------------------- We first formulate the Bayesian inverse problem using the first data distribution and analyze the posterior samples. .. GENERATED FROM PYTHON SOURCE LINES 169-171 Create the posterior ~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 171-174 .. code-block:: Python z1 = cuqi.distribution.JointDistribution(theta,y1)(y1=data1) .. GENERATED FROM PYTHON SOURCE LINES 175-176 We print the joint distribution `z1`: .. GENERATED FROM PYTHON SOURCE LINES 176-178 .. code-block:: Python print(z1) .. rst-class:: sphx-glr-script-out .. code-block:: none Posterior( Equation: p(theta|y1) ∝ L(theta|y1)p(theta) Densities: y1 ~ CUQI Gaussian Likelihood function. Parameters ['theta']. theta ~ CUQI Gaussian. ) .. GENERATED FROM PYTHON SOURCE LINES 179-183 We see that we obtain a :class:`cuqi.distribution.Posterior` object, which represents the posterior distribution of the parameters `theta` given the data `y1`. The posterior distribution in this case is proportional to the product of the likelihood obtained from the first data distribution and the prior. .. GENERATED FROM PYTHON SOURCE LINES 185-187 Sample from the posterior ~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 187-197 .. code-block:: Python # Sample from the posterior sampler = cuqi.sampler.MH(z1) sampler.warmup(8000) samples = sampler.get_samples() # Plot the credible interval and compute the ESS samples.burnthin(1000).plot_ci(95, exact=problemInfo1.exactSolution) samples.compute_ess() .. image-sg:: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_003.png :alt: multiple likelihoods :srcset: /user/_auto_howtos/images/sphx_glr_multiple_likelihoods_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Warmup: 0%| | 0/8000 [00:00` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: multiple_likelihoods.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: multiple_likelihoods.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_