.. _bayesian_reconstruction:

Bayesian parameter reconstruction using Bayesian optimization
==========================================================================================


:Driver: :ref:`BayesianReconstruction`


:Download script: :download:`bayesian_reconstruction.py`


The target of the study is to showcase the solution of a Bayesian parameter reconstruction problembased on a set ot measurements :math:`\mathbf{t} \in \mathbb{R}^K`.

We assume that the measurement process can be accurately modeled by a
function :math:`\mathbf{f}(\mathbf{b}) \in \mathbb{R}^K`.
As an example we consider the `MGH17 problem <https://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml>`_ which consists of fitting a vectorial model :math:`\mathbf{f}(\mathbf{b}) \in \mathbb{R}^{33}` with

.. math::
   f_i(\mathbf{b}) = b_1 + b_2 \exp(-i \cdot b_4) + b_3 \exp(-i \cdot b_5),\, i = 0, \dots, 32

to a measurement vector with :math:`K = 33` entries.

Due to random measurement noise :math:`\mathbf{w} \sim \mathcal{N}(0,
{\rm diag}{[\eta_1^2,...,\eta_K^2]})`, the
vector of measurements :math:`\mathbf{t} = \mathbf{f}(\mathbf{b}) + \mathbf{w}`
is a random vector with probability density

.. math::
   P(\mathbf{t}) = \prod_{i=1}^K
   \frac{1}{\sqrt{2\pi}\eta_i}\exp\left[-\frac{1}{2}\left(\frac{t_i -
   f_i(\mathbf{b})}{\eta_i}\right)^2\right].

The noise variances :math:`\eta_i^2` shall be unknown and estimated
from the measurement data. To this end we assume that the variance is
composed of a background term :math:`c_1^2` and a noise contribution which
scales linearly with :math:`f_i(\mathbf{b})`, i.e.

.. math::

   \eta_i^2(\mathbf{c}) = c_1^2 + \left[c_2 f_i(\mathbf{b})\right]^2.
	   
Taking non-uniform *prior distributions* for the design parameter
vector :math:`P_\text{prior}(\mathbf{b})` and the error model
parameters :math:`P_\text{prior}(\mathbf{c})` into account the *posterior
distribution* is then proportional to :math:`P(\mathbf{b}, \mathbf{c} |
\mathbf{t}) \propto P(\mathbf{t} | \mathbf{b}, \mathbf{c})
P_\text{prior}(\mathbf{b}) P_\text{prior}(\mathbf{c})`.
	   
Alltogether, the target of finding the parameters with maximum
posterior probability density is equivalent of minimizing the value of
the negative log-probability

.. math::
   \begin{split}
   -\log\left(P(\mathbf{b}, \mathbf{c}| \mathbf{t})\right) = &
   \frac{1}{2} K\log(2\pi)
   +\sum_{i=1}^K\log\left(\eta_i(\mathbf{c})\right) 
   +\frac{1}{2}\sum_{i=1}^K \left(
      \frac{t_i - f_i(\mathbf{b})}{\eta_i(\mathbf{c})}
   \right)^2 \\
   &-\log\left(P_\text{prior}(\mathbf{d})\right)
   -\log\left(P_\text{prior}(\mathbf{c})\right).
           \end{split}

In the following, we show how this negative log-probability can be minimized using
the :ref:`BayesianReconstruction` driver. We compare the results of the dirver's
MCMC sampling of the posterior probability to an MCMC sampling of the analytic value.



.. literalinclude:: ./bayesian_reconstruction.py
   :language: python
   :linenos:


 

.. figure:: images/bayesian_reconstruction/corner_analytic.svg
   :alt: MCMC sampling based on surrogate

   Markov-Chain Monte-Carlo (MCMC) sampling of the probability density
   of the model parameters :math:`b_1,\dots,b_5` and the log value of
   the error model parameters :math:`c_1,\dots,c_2` based on the
   analytic model function.

.. figure:: images/bayesian_reconstruction/corner_surrogate.svg
   :alt: MCMC sampling based on surrogate

   Markov-Chain Monte-Carlo (MCMC) sampling of the probability density
   of the model parameters :math:`b_1,\dots,b_5` and the log value of
   the error model parameters :math:`c_1,\dots,c_2` based on the
   trained surrogate of the study. A comparison between the analytic
   and the surrogate model function shows a good quantitative
   agreement.