Under the hood
==============
What is the code doing while it's running? Let's take a peek under the hood.
Parameterization
----------------
Models in **Bayesian1DEM**, are parameterized as a stack of layers, each of which is assigned a constant value of
bulk resistivity. The layers themselves are defined as sequence of layer interface depths. Due to the diffusive
nature of most EM methods, as well as the wide range of bulk resistivities in Earth materials, it is often best to
invert for the base-10 logarithm of both the layer resistivities and the layer interface depths. **Bayesian1DEM**
omits the interface at z=0 and assumes that the layer underneath the final interface is a halfspace.
Following [Blatter19]_, **Bayesian1DEM** samples the subsurface using two model regions, each with their own inversion
parameters. This was designed for inversion of marine CSEM data where the properties of the water column are vastly
different than the properties of the subsurface (yet getting the water column resistivity right is crucial). Thus, by
default model region 1 is above model region 2 in depth and represents the water column, while model region 2 represents
the subsurface. **Bayesian1DEM** samples these two regions separately (the relative frequency determined by an input
parameter), but the forward modeling of course uses the full, composite model.
Markov chains
-------------
**Bayesian1DEM** uses Markov chains, which are sequences of models that grow iteratively by adding models to the end of the
chain. At each step in the algorithm, **Bayesian1DEM** generates a new model (called the "proposal") by modifying
the current model in one of four ways:
* Adding a new layer interface at a random location
* Deleting a current layer interface (chosen at random)
* Moving a randomly chosen existing layer interface to a new location (the new location is chosen by a Gaussian
distribution whose mean is set to the interface's current location)
* Changing the resistivity of each layer using a Gaussian distribution with mean equal to each layer's
current resistivity value
**Bayesian1DEM** then accepts or rejects the proposal model on the basis of how well it fits the data and
our prior assumptions (its posterior probability) compared to the current model, using the standard Metropolis-Hastings
accept-reject criterion. If the proposal is accepted, it becomes the next model in the chain. If it is rejected,
the current model is copied and becomes the next model in the chain as well. The algorithm then repeats, generating
a new proposal model from the new current model.
Output: the model ensemble
--------------------------
After the chain is sufficiently long such that the statistics of the model parameters become stationary (they aren't
changing with time) and the initial stage as the chain finds models that fit the data (the "burn-in") is removed,
the models in the Markov chain represent samples drawn from the Bayesian postrior distribution. In other words, the
statistics of the model parameters across all the models in the chain (the "model ensemble") represent the Bayesian
posterior uncertainty. Thus, once **Bayesian1DEM** has finished running, the output is a collection of models (saved
in the output .mat files under the variable *s_ll*) that can be mined for all sorts of statistical information: mean,
median, variance, interquartile range, 90% credible interval, etc.
Flexibility
-----------
**Bayesian1DEM** is both physical property and data-type invariant. The core of the algorithm knows only that the
"model" is a stack of layers defined by a series of layer interfaces with one or more constant properties ascribed
to each layer. So long as you have a forward modeling code that you can call to generate modeled data given the current
model as defined by **Bayesian1DEM**, you are good to go. This means you can rather easily modify **Bayesian1DEM** to
invert whatever data you want, so long as a trans-dimensional, piecewise-constant stack of layers is an appropriate
model parameterization.