Linear regression with Eustatic Sea-level data

Note

Go to the end to download the full example code.

Linear regression with Eustatic Sea-level data#

Open In Colab

If you are running this notebook locally, make sure you’ve followed steps here to set up the environment. (This environment.yml file specifies a list of packages required to run the notebooks)

What we do in this notebook#

Here we demonstrate use of CoFI on a real dataset linear regression problem, where we fit a polynomial function to Eustatic Sea-level heights.

  • by solution of a linear system of equations,

  • by optimization of a data misfit function

  • by Bayesian sampling of a Likelihood multiplied by a prior.

Data set is from

Sea level and global ice volumes from the Last Glacial Maximum to the Holocene, K. Lambeck, H. Rouby, A. Purcell, Y. Sun, and M. Sambridge, 2014. Proc. Nat. Acad. Sci., 111, no. 43, 15296-15303, doi:10.1073/pnas.1411762111.

# Environment setup (uncomment code below)

# !pip install -U cofi
# !git clone https://github.com/inlab-geo/cofi-examples.git
# %cd cofi-examples/examples/linear_regression

Linear regression#

Lets start with some (x,y) data.

import numpy as np
import matplotlib.pyplot as plt

#
# load x,y data
def load_data_xy(filename):
    x, y, sx, sy = np.loadtxt(filename, skiprows=1).T
    sy /= 2  # input data are 2 sigma and we require sigma
    return x, y, sy

def load_data_ref(filename):

    f = open(filename, 'r')
    lines = f.readlines()
    dx = np.array([]) # Age data
    dy = np.array([]) # ESL height
    dz = np.array([]) # derivative of ESL w.r.t. age
    for line in lines:
        columns = line.split()
        dx = np.append(dx,float(columns[0]))
        dy = np.append(dy,float(columns[1]))
    datavals = np.column_stack((dx,dy))  # Stack data

    return datavals

data_x,data_y,sy = load_data_xy("../../data/eustatic_sea_level/ESL-ff11-sorted.txt")  # Load x,sx,y,sy ESL data (x=time, Y=ESL)

maxtime = 20.
ndata = np.where(data_x>maxtime)[0][0]

data_x,data_y,sy = data_x[:ndata],data_y[:ndata],sy[:ndata]

And now lets plot the data.

def plot_data(x=data_x,y=data_y,sigma=sy,title=None):
    fig, axes = plt.subplots(figsize=(9,6))
    plt.errorbar(x, y, yerr=sy, fmt='.',color="lightcoral",ecolor='lightgrey',ms=2)
    plt.xlabel(' Time before present (ka)')
    plt.ylabel(' ESL height (m)')
    if(title != None): plt.title(title)
plot_data(title='Eustatic sea-level')

# display theory on travel time tomography
from IPython.display import display, Markdown

with open("../../theory/linear_regression.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

We now build the Jacobian/G matrix for this problem and define a forward function which simply multiplies \(\mathbf m\) by \(G\).

nparams = 5 # Number of model parameters to be solved for

def jacobian(x=data_x, n=nparams):
    return np.array([x**i for i in range(n)]).T

def forward(model):
    return jacobian().dot(model)

def Cd_inv(sigma=sy):
    factor= 10                                   # factor to inflate observational errors
    return np.diag(1./sy*1./sy)/(factor**2)

Define a reference model for later.

# Reference model for plotting
ESLref = load_data_ref("../../data/eustatic_sea_level/ESL-f11_yonly.txt") # Load x, y, z reference model and estimated derivative (x=time, Y=ESL, z=dESL/dt)
ndata2 = np.where(ESLref.T[0]>maxtime)[0][0]
ESLref = ESLref[:ndata2]
ref_x,ref_y = ESLref.T[0],ESLref.T[1]

Now lets plot the data with the reference curve

# Some plotting utilities
def plot_model(x,y, label, color=None,lw=0.5):
    plt.plot(x, y, color=color or "green", label=label,lw=lw)
    #plt.xlabel("X")
    #plt.ylabel("ESL")
    plt.legend()

def plot_models(models, label="Posterior samples", color="seagreen", alpha=0.1,lw=0.5):
    G = jacobian(data_x)
    plt.plot(data_x, G.dot(models[0]), color=color, label=label, alpha=alpha,lw=lw)
    for m in models:
        plt.plot(data_x, G.dot(m), color=color, alpha=alpha,lw=lw)
    plt.legend()

plot_data(title="Eustatic sea-level")
plot_model(ref_x,ref_y, "Reference model")

Now we have the data and the forward model we can start to try and estimate the coefficients of the polynomial from the data.

The structure of CoFI#

In the workflow of cofi, there are three main components: BaseProblem, InversionOptions, and Inversion.

  • BaseProblem defines the inverse problem including any user supplied quantities such as data vector, number of model parameters and measure of fit between model predictions and data. python     inv_problem = BaseProblem()     inv_problem.set_objective(some_function_here)     inv_problem.set_jacobian(some_function_here)     inv_problem.set_initial_model(a_starting_point) # if needed, e.g. we are solving a nonlinear problem by optimization

  • InversionOptions describes details about how one wants to run the inversion, including the backend tool and solver-specific parameters. It is based on the concept of a method and tool.

    inv_options = InversionOptions()
inv_options.suggest_solving_methods()
inv_options.set_solving_method("matrix solvers")
inv_options.suggest_tools()
inv_options.set_tool("scipy.linalg.lstsq")
inv_options.summary()

  • Inversion can be seen as an inversion engine that takes in the above two as information, and will produce an InversionResult upon running.

    inv = Inversion(inv_problem, inv_options)
result = inv.run()

Internally CoFI decides the nature of the problem from the quantities set by the user and performs internal checks to ensure it has all that it needs to solve a problem.

1. Linear system solver#

from cofi import BaseProblem, InversionOptions, Inversion

Step 1. Define CoFI BaseProblem#

inv_problem = BaseProblem()
inv_problem.set_data(data_y)
inv_problem.set_jacobian(jacobian())
inv_problem.set_data_covariance_inv(Cd_inv())

Step 2. Define CoFI InversionOptions#

inv_options = InversionOptions()

Using the information supplied, we can ask CoFI to suggest some solving methods.

inv_options.suggest_solving_methods()

The following solving methods are supported:
{'matrix solvers', 'optimization', 'sampling'}

Use `suggest_tools()` to see a full list of backend tools for each method

We can ask CoFI to suggest some specific software tools as well.

inv_options.suggest_tools()

Here's a complete list of inversion tools supported by CoFI (grouped by methods):
{
    "optimization": [
        "scipy.optimize.minimize",
        "scipy.optimize.least_squares",
        "torch.optim",
        "cofi.border_collie_optimization",
        "neighpyI",
        "mealpy.sma",
        "mealpy.slime_mould"
    ],
    "matrix solvers": [
        "scipy.linalg.lstsq",
        "cofi.simple_newton",
        "scipy.sparse.linalg"
    ],
    "sampling": [
        "emcee",
        "bayesbay",
        "neighpy",
        "neighpyII"
    ]
}

inv_options.set_solving_method("matrix solvers") # lets decide to use a matrix solver.
inv_options.summary()

=============================
Summary for inversion options
=============================
Solving method: matrix solvers
Use `suggest_solving_methods()` to check available solving methods.
-----------------------------
Backend tool: `<class 'cofi.tools._scipy_lstsq.ScipyLstSq'> (by default)` - SciPy's wrapper function over LAPACK's linear least-squares solver, using 'gelsd', 'gelsy' (default), or 'gelss' as backend driver
References: ['https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html', 'https://www.netlib.org/lapack/lug/node27.html']
Use `suggest_tools()` to check available backend tools.
-----------------------------
Solver-specific parameters: None set
Use `suggest_solver_params()` to check required/optional solver-specific parameters.

# below is optional, as this has already been the default tool under "linear least square"
inv_options.set_tool("scipy.linalg.lstsq")

Step 3. Define CoFI Inversion and run#

Our choices so far have defined a linear parameter estimation problem (without any regularization) to be solved within a least squares framework. In this case the selection of a matrix solvers method will mean we are calculating the standard least squares solution

\[m = (G^T C_d^{-1} G)^{-1} G^T C_d^{-1} d\]

and our choice of backend tool scipy.linalg.lstsq, means that we will employ scipy’s linalg package to perform the numerics.

Lets run CoFI.

inv = Inversion(inv_problem, inv_options)
inv_result = inv.run()

print(f"The inversion result from `scipy.linalg.lstsq`: {inv_result.model}\n")
inv_result.summary()

The inversion result from `scipy.linalg.lstsq`: [ 1.44051039 -3.11381469  1.412872   -0.20910136  0.00653572]

============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [ 1.44051039 -3.11381469  1.412872   -0.20910136  0.00653572]
sum_of_squared_residuals: []
effective_rank: 5
singular_values: [1.72339368e+09 1.35569994e+06 3.54123758e+03 1.10384400e+02
 7.16903392e+00]
model_covariance: [[ 6.70145209e-02 -6.28173834e-02  1.67266578e-02 -1.60867740e-03
   4.77025592e-05]
 [-6.28173834e-02  7.43365995e-02 -2.23488563e-02  2.32179303e-03
  -7.20214235e-05]
 [ 1.67266578e-02 -2.23488563e-02  7.38371078e-03 -8.20763652e-04
   2.65025504e-05]
 [-1.60867740e-03  2.32179303e-03 -8.20763652e-04  9.62437921e-05
  -3.21300646e-06]
 [ 4.77025592e-05 -7.20214235e-05  2.65025504e-05 -3.21300646e-06
   1.10114955e-07]]

Lets plot the solution.

plot_data(title="Eustatic sea-level")
plot_model(data_x,jacobian(data_x).dot(inv_result.model), "linear system solver", color="seagreen")
plot_model(ref_x,ref_y, "Reference model", color="darkorange")

2. Optimizer#

The same overdetermined linear problem, \(\textbf{d} = G\textbf{m}\), with Gaussian data noise can also be solved by minimising the squares of the residual of the linear equations, e.g. \(\textbf{r}^T \textbf{C}_d^{-1}\textbf{r}\) where \(\textbf{r}=\textbf{d}-G\textbf{m}\). The above matrix solver solution gives us the best data fitting model, but a direct optimisation approach could also be used, say when the number of unknowns is large and we do not wish, or are unable to provide the Jacobian function.

So we use a plain optimizer scipy.optimize.minimize to demonstrate this ability.

######## CoFI BaseProblem - provide additional information
inv_problem.set_initial_model(np.ones(nparams))
#inv_problem.set_initial_model(inv_result.model)
inv_problem.set_forward(forward)
inv_problem.set_data_misfit("squared error")

# inv_problem.set_objective(your_own_misfit_function)    # (optionally) if you'd like to define your own misfit
# inv_problem.set_gradient(your_own_gradient_of_misfit_function)    # (optionally) if you'd like to define your own misfit gradient

######## CoFI InversionOptions - set a different tool
inv_options_2 = InversionOptions()
inv_options_2.set_tool("scipy.optimize.minimize")
inv_options_2.set_params(method="Nelder-Mead")

######## CoFI Inversion - run it
inv_2 = Inversion(inv_problem, inv_options_2)
inv_result_2 = inv_2.run()

######## CoFI InversionResult - check result
print(f"The inversion result from `scipy.optimize.minimize`: {inv_result_2.model}\n")
inv_result_2.summary()

The inversion result from `scipy.optimize.minimize`: [-0.81319113 -0.46081553  0.61740809 -0.12663625  0.00398136]

============================
Summary for inversion result
============================
SUCCESS
----------------------------
fun: 310.0689929995379
nit: 560
nfev: 916
status: 0
message: Optimization terminated successfully.
final_simplex: (array([[-0.81319113, -0.46081553,  0.61740809, -0.12663625,  0.00398136],
       [-0.81328903, -0.46081801,  0.61741429, -0.1266371 ,  0.0039814 ],
       [-0.81317886, -0.46081571,  0.61738786, -0.12663294,  0.00398125],
       [-0.81315239, -0.46081502,  0.61739705, -0.12663469,  0.0039813 ],
       [-0.81320724, -0.46081746,  0.61741117, -0.12663795,  0.00398147],
       [-0.81318091, -0.4608153 ,  0.61740441, -0.12663565,  0.00398136]]), array([310.068993  , 310.06899312, 310.0689937 , 310.06899379,
       310.06899388, 310.06899454]))
model: [-0.81319113 -0.46081553  0.61740809 -0.12663625  0.00398136]

plot_data()
plot_model(data_x,jacobian(data_x).dot(inv_result_2.model), "optimization solution", color="cornflowerblue")
plot_model(ref_x,ref_y, "Reference model", color="darkorange")

The optimization fails to convergence for this problem (with default settings).

Challenge - Change the polynomial degree#

Try and replace the 3rd order polynomial with a 2nd order polynomial (i.e. \(M=2\)) by adding the required commands below. What does the plot looks like?

Start from code below:

inv_problem = BaseProblem()
inv_problem.set_data(data_y)
inv_problem.set_jacobian(jacobian(n=<CHANGE ME>))
inv_problem.set_data_covariance_inv(Cd_inv())
inv_options.set_solving_method("matrix solvers") # lets decide to use a matrix solver.
inv = Inversion(inv_problem, inv_options)
inv_result = inv.run()

print("Inferred curve with n = <CHANGE ME> ")
plot_data()
plot_model(data_x,jacobian(x,n=<CHANGE ME>).dot(inv_result.model), "optimization solution", color="cornflowerblue")
plot_model(ref_x,ref_y, "Reference model", color="darkorange")

# Copy the template above, Replace <CHANGE ME> with your answer

#@title Solution

inv_problem = BaseProblem()
inv_problem.set_data(data_y)
inv_problem.set_jacobian(jacobian(n=3))
inv_problem.set_data_covariance_inv(Cd_inv())
inv_options.set_solving_method("matrix solvers") # lets decide to use a matrix solver.
inv = Inversion(inv_problem, inv_options)
inv_result = inv.run()

print("Inferred curve with n = 3 ")
plot_data()
plot_model(data_x,jacobian(data_x,n=3).dot(inv_result.model), "optimization solution", color="cornflowerblue")
plot_model(ref_x,ref_y, "Reference model", color="darkorange")

Inferred curve with n = 3

Changing to a second order polynomial does converge but gives a poor fit.

3. Bayesian sampling#

Likelihood#

Since data errors follow a Gaussian in this example, we can define a Likelihood function, \(p({\mathbf d}_{obs}| {\mathbf m})\).

\[p({\mathbf d}_{obs} | {\mathbf m}) \propto \exp \left\{- \frac{1}{2} ({\mathbf d}_{obs}-{\mathbf d}_{pred}({\mathbf m}))^T C_D^{-1} ({\mathbf d}_{obs}-{\mathbf d}_{pred}({\mathbf m})) \right\}\]

where \({\mathbf d}_{obs}\) represents the observed y values and \({\mathbf d}_{pred}({\mathbf m})\) are those predicted by the polynomial model \(({\mathbf m})\). The Likelihood is defined as the probability of observing the data actually observed, given a model. In practice we usually only need to evaluate the log of the Likelihood, \(\log p({\mathbf d}_{obs} | {\mathbf m})\). To do so, we require the inverse data covariance matrix describing the statistics of the noise in the data, \(C_D^{-1}\) . For this problem the data errors are independent with identical standard deviation in noise for each datum. Hence \(C_D^{-1} = \frac{1}{\sigma^2}I\) where \(\sigma=1\).

Here we artificially increase the observational errors on the data so that the spread of the posterior samples are visible.

Cdinv = Cd_inv()                 # inverse data covariance matrix

def log_likelihood(model):
    y_synthetics = forward(model)
    residual = data_y - y_synthetics
    return -0.5 * residual @ (Cdinv @ residual).T

Note that the user could specify any appropriate Likelihood function of their choosing here.

Prior#

Bayesian sampling requires a prior probability density function. A common problem with polynomial coefficients as model parameters is that it is not at all obvious what a prior should be. Here we choose a uniform prior with specified bounds

\[\begin{split}\begin{align} p({\mathbf m}) &= \frac{1}{V},\quad l_i \le m_i \le u_i, \quad (i=1,\dots,M)\\ \\ &= 0, \quad {\rm otherwise}, \end{align}\end{split}\]

where \(l_i\) and \(u_i\) are lower and upper bounds on the \(i\)th model coefficient.

Here use the uniform distribution with \({\mathbf l}^T = (-10.,-10.,-10.,-10.)\), and \({\mathbf u}^T = (10.,10.,10.,10.)\).

m_lower_bound = np.ones(nparams) * (-10.)             # lower bound for uniform prior
m_upper_bound = np.ones(nparams) * 10                 # upper bound for uniform prior

def log_prior(model):    # uniform distribution
    for i in range(len(m_lower_bound)):
        if model[i] < m_lower_bound[i] or model[i] > m_upper_bound[i]: return -np.inf
    return 0.0 # model lies within bounds -> return log(1)

Note that the user could specify any appropriate Prior PDF of their choosing here.

Bayesian sampling#

In this aproach we sample a probability distribution rather than find a single best fit solution. Bayes’ theorem tells us the the posterior distribution is proportional to the Likelihood and the prior.

\[p(\mathbf{m}|\mathbf{d}) = K p(\mathbf{d}|\mathbf{m})p(\mathbf{m})\]

where \(K\) is some constant. Under the assumptions specified \(p(\mathbf{m}|\mathbf{d})\) gives a probability density of models that are supported by the data. We seek to draw random samples from \(p(\mathbf{m}|\mathbf{d})\) over model space and then to make inferences from the resulting ensemble of model parameters.

In this example we make use of The Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler Goodman and Weare 2010 to sample the posterior distribution of the model. (See more details about emcee).

Starting points for random walkers#

Now we define some hyperparameters (e.g. the number of walkers and steps), and initialise the starting positions of walkers. We start all walkers in a small ball about a chosen point \((0, 0, 0, 0)\).

nwalkers = 32
ndim = nparams
nsteps = 10000
walkers_start = np.zeros(nparams) + 1e-4 * np.random.randn(nwalkers, ndim)

Add the information and run with CoFI#

######## CoFI BaseProblem - provide additional information
inv_problem.set_log_prior(log_prior)
inv_problem.set_log_likelihood(log_likelihood)
inv_problem.set_model_shape(ndim)

######## CoFI InversionOptions - get a different tool
inv_options_3 = InversionOptions()
inv_options_3.set_tool("emcee")      # Here we use to Affine Invariant McMC sampler from Goodman and Weare (2010).
inv_options_3.set_params(nwalkers=nwalkers, nsteps=nsteps, progress=True, initial_state=walkers_start)

######## CoFI Inversion - run it
inv_3 = Inversion(inv_problem, inv_options_3)
inv_result_3 = inv_3.run()

######## CoFI InversionResult - check result
print(f"The inversion result from `emcee`:")
inv_result_3.summary()

  0%|          | 0/10000 [00:00<?, ?it/s]
  0%|          | 36/10000 [00:00<00:28, 353.83it/s]
  1%|          | 72/10000 [00:00<00:28, 350.64it/s]
  1%|          | 108/10000 [00:00<00:28, 341.30it/s]
  1%|▏         | 143/10000 [00:00<00:37, 263.51it/s]
  2%|▏         | 176/10000 [00:00<00:34, 282.95it/s]
  2%|▏         | 206/10000 [00:00<00:35, 274.04it/s]
  2%|▏         | 235/10000 [00:00<00:39, 245.90it/s]
  3%|▎         | 271/10000 [00:00<00:35, 275.25it/s]
  3%|▎         | 307/10000 [00:01<00:32, 297.44it/s]
  3%|▎         | 343/10000 [00:01<00:30, 314.00it/s]
  4%|▍         | 379/10000 [00:01<00:29, 325.96it/s]
  4%|▍         | 415/10000 [00:01<00:28, 334.70it/s]
  5%|▍         | 451/10000 [00:01<00:28, 340.94it/s]
  5%|▍         | 487/10000 [00:01<00:27, 345.54it/s]
  5%|▌         | 523/10000 [00:01<00:27, 348.65it/s]
  6%|▌         | 559/10000 [00:01<00:26, 351.01it/s]
  6%|▌         | 595/10000 [00:01<00:27, 345.95it/s]
  6%|▋         | 631/10000 [00:01<00:26, 347.86it/s]
  7%|▋         | 667/10000 [00:02<00:26, 350.14it/s]
  7%|▋         | 703/10000 [00:02<00:27, 343.41it/s]
  7%|▋         | 739/10000 [00:02<00:26, 346.78it/s]
  8%|▊         | 775/10000 [00:02<00:26, 349.12it/s]
  8%|▊         | 811/10000 [00:02<00:26, 349.61it/s]
  8%|▊         | 847/10000 [00:02<00:26, 351.43it/s]
  9%|▉         | 883/10000 [00:02<00:25, 351.97it/s]
  9%|▉         | 919/10000 [00:02<00:25, 353.21it/s]
 10%|▉         | 955/10000 [00:02<00:25, 353.79it/s]
 10%|▉         | 991/10000 [00:03<00:25, 353.62it/s]
 10%|█         | 1027/10000 [00:03<00:25, 354.51it/s]
 11%|█         | 1063/10000 [00:03<00:25, 352.65it/s]
 11%|█         | 1099/10000 [00:03<00:25, 353.93it/s]
 11%|█▏        | 1135/10000 [00:03<00:25, 353.27it/s]
 12%|█▏        | 1171/10000 [00:03<00:25, 352.48it/s]
 12%|█▏        | 1207/10000 [00:03<00:24, 353.46it/s]
 12%|█▏        | 1243/10000 [00:03<00:24, 353.97it/s]
 13%|█▎        | 1279/10000 [00:03<00:24, 353.09it/s]
 13%|█▎        | 1315/10000 [00:03<00:24, 349.16it/s]
 14%|█▎        | 1351/10000 [00:04<00:24, 350.78it/s]
 14%|█▍        | 1387/10000 [00:04<00:24, 350.38it/s]
 14%|█▍        | 1423/10000 [00:04<00:24, 351.62it/s]
 15%|█▍        | 1459/10000 [00:04<00:24, 349.22it/s]
 15%|█▍        | 1495/10000 [00:04<00:24, 351.10it/s]
 15%|█▌        | 1531/10000 [00:04<00:24, 351.77it/s]
 16%|█▌        | 1567/10000 [00:04<00:23, 352.73it/s]
 16%|█▌        | 1603/10000 [00:04<00:23, 352.29it/s]
 16%|█▋        | 1639/10000 [00:04<00:23, 353.02it/s]
 17%|█▋        | 1675/10000 [00:04<00:23, 353.41it/s]
 17%|█▋        | 1711/10000 [00:05<00:23, 352.10it/s]
 17%|█▋        | 1747/10000 [00:05<00:23, 351.78it/s]
 18%|█▊        | 1783/10000 [00:05<00:23, 352.33it/s]
 18%|█▊        | 1819/10000 [00:05<00:23, 352.39it/s]
 19%|█▊        | 1855/10000 [00:05<00:25, 322.92it/s]
 19%|█▉        | 1888/10000 [00:05<00:27, 289.72it/s]
 19%|█▉        | 1921/10000 [00:05<00:27, 298.34it/s]
 20%|█▉        | 1957/10000 [00:05<00:25, 313.38it/s]
 20%|█▉        | 1993/10000 [00:05<00:24, 324.64it/s]
 20%|██        | 2029/10000 [00:06<00:23, 332.56it/s]
 21%|██        | 2065/10000 [00:06<00:23, 338.97it/s]
 21%|██        | 2101/10000 [00:06<00:23, 342.74it/s]
 21%|██▏       | 2136/10000 [00:06<00:22, 344.33it/s]
 22%|██▏       | 2172/10000 [00:06<00:22, 347.22it/s]
 22%|██▏       | 2208/10000 [00:06<00:22, 349.31it/s]
 22%|██▏       | 2244/10000 [00:06<00:22, 351.16it/s]
 23%|██▎       | 2280/10000 [00:06<00:21, 352.47it/s]
 23%|██▎       | 2316/10000 [00:06<00:21, 353.24it/s]
 24%|██▎       | 2352/10000 [00:06<00:21, 353.93it/s]
 24%|██▍       | 2388/10000 [00:07<00:21, 353.48it/s]
 24%|██▍       | 2424/10000 [00:07<00:21, 352.58it/s]
 25%|██▍       | 2460/10000 [00:07<00:21, 353.49it/s]
 25%|██▍       | 2496/10000 [00:07<00:21, 353.45it/s]
 25%|██▌       | 2532/10000 [00:07<00:21, 353.85it/s]
 26%|██▌       | 2568/10000 [00:07<00:21, 350.39it/s]
 26%|██▌       | 2604/10000 [00:07<00:21, 352.15it/s]
 26%|██▋       | 2640/10000 [00:07<00:21, 349.84it/s]
 27%|██▋       | 2676/10000 [00:07<00:20, 350.67it/s]
 27%|██▋       | 2712/10000 [00:07<00:20, 348.65it/s]
 27%|██▋       | 2748/10000 [00:08<00:20, 350.07it/s]
 28%|██▊       | 2784/10000 [00:08<00:20, 350.08it/s]
 28%|██▊       | 2820/10000 [00:08<00:20, 350.73it/s]
 29%|██▊       | 2856/10000 [00:08<00:20, 351.26it/s]
 29%|██▉       | 2892/10000 [00:08<00:20, 351.83it/s]
 29%|██▉       | 2928/10000 [00:08<00:20, 351.35it/s]
 30%|██▉       | 2964/10000 [00:08<00:20, 351.67it/s]
 30%|███       | 3000/10000 [00:08<00:19, 351.78it/s]
 30%|███       | 3036/10000 [00:08<00:19, 352.36it/s]
 31%|███       | 3072/10000 [00:09<00:19, 352.36it/s]
 31%|███       | 3108/10000 [00:09<00:19, 352.93it/s]
 31%|███▏      | 3144/10000 [00:09<00:19, 353.08it/s]
 32%|███▏      | 3180/10000 [00:09<00:19, 352.78it/s]
 32%|███▏      | 3216/10000 [00:09<00:19, 352.98it/s]
 33%|███▎      | 3252/10000 [00:09<00:19, 353.48it/s]
 33%|███▎      | 3288/10000 [00:09<00:19, 352.23it/s]
 33%|███▎      | 3324/10000 [00:09<00:18, 352.84it/s]
 34%|███▎      | 3360/10000 [00:09<00:18, 353.39it/s]
 34%|███▍      | 3396/10000 [00:09<00:18, 354.13it/s]
 34%|███▍      | 3432/10000 [00:10<00:18, 352.16it/s]
 35%|███▍      | 3468/10000 [00:10<00:20, 325.87it/s]
 35%|███▌      | 3501/10000 [00:10<00:27, 239.11it/s]
 35%|███▌      | 3533/10000 [00:10<00:25, 256.37it/s]
 36%|███▌      | 3562/10000 [00:10<00:24, 258.12it/s]
 36%|███▌      | 3590/10000 [00:10<00:26, 241.64it/s]
 36%|███▌      | 3623/10000 [00:10<00:24, 261.78it/s]
 37%|███▋      | 3651/10000 [00:10<00:25, 253.08it/s]
 37%|███▋      | 3678/10000 [00:11<00:27, 230.85it/s]
 37%|███▋      | 3712/10000 [00:11<00:24, 257.11it/s]
 37%|███▋      | 3748/10000 [00:11<00:22, 282.97it/s]
 38%|███▊      | 3784/10000 [00:11<00:20, 301.71it/s]
 38%|███▊      | 3820/10000 [00:11<00:19, 316.13it/s]
 39%|███▊      | 3856/10000 [00:11<00:18, 326.36it/s]
 39%|███▉      | 3892/10000 [00:11<00:18, 334.23it/s]
 39%|███▉      | 3928/10000 [00:11<00:17, 340.57it/s]
 40%|███▉      | 3964/10000 [00:11<00:17, 344.44it/s]
 40%|████      | 4000/10000 [00:12<00:17, 347.33it/s]
 40%|████      | 4036/10000 [00:12<00:17, 349.61it/s]
 41%|████      | 4072/10000 [00:12<00:16, 351.09it/s]
 41%|████      | 4108/10000 [00:12<00:16, 351.17it/s]
 41%|████▏     | 4144/10000 [00:12<00:16, 351.73it/s]
 42%|████▏     | 4180/10000 [00:12<00:16, 352.58it/s]
 42%|████▏     | 4216/10000 [00:12<00:16, 353.53it/s]
 43%|████▎     | 4252/10000 [00:12<00:16, 353.93it/s]
 43%|████▎     | 4288/10000 [00:12<00:16, 354.72it/s]
 43%|████▎     | 4324/10000 [00:12<00:15, 355.23it/s]
 44%|████▎     | 4360/10000 [00:13<00:15, 355.25it/s]
 44%|████▍     | 4396/10000 [00:13<00:15, 355.26it/s]
 44%|████▍     | 4432/10000 [00:13<00:15, 354.38it/s]
 45%|████▍     | 4468/10000 [00:13<00:15, 354.62it/s]
 45%|████▌     | 4504/10000 [00:13<00:15, 355.06it/s]
 45%|████▌     | 4540/10000 [00:13<00:15, 355.07it/s]
 46%|████▌     | 4576/10000 [00:13<00:15, 354.43it/s]
 46%|████▌     | 4612/10000 [00:13<00:15, 354.38it/s]
 46%|████▋     | 4648/10000 [00:13<00:15, 354.83it/s]
 47%|████▋     | 4684/10000 [00:13<00:15, 354.15it/s]
 47%|████▋     | 4720/10000 [00:14<00:14, 352.45it/s]
 48%|████▊     | 4756/10000 [00:14<00:14, 353.38it/s]
 48%|████▊     | 4792/10000 [00:14<00:14, 353.81it/s]
 48%|████▊     | 4828/10000 [00:14<00:14, 354.39it/s]
 49%|████▊     | 4864/10000 [00:14<00:14, 354.78it/s]
 49%|████▉     | 4900/10000 [00:14<00:14, 353.77it/s]
 49%|████▉     | 4936/10000 [00:14<00:14, 354.37it/s]
 50%|████▉     | 4972/10000 [00:14<00:14, 354.63it/s]
 50%|█████     | 5008/10000 [00:14<00:14, 355.00it/s]
 50%|█████     | 5044/10000 [00:14<00:13, 355.28it/s]
 51%|█████     | 5080/10000 [00:15<00:13, 355.00it/s]
 51%|█████     | 5116/10000 [00:15<00:13, 349.92it/s]
 52%|█████▏    | 5152/10000 [00:15<00:13, 351.77it/s]
 52%|█████▏    | 5188/10000 [00:15<00:13, 352.86it/s]
 52%|█████▏    | 5224/10000 [00:15<00:13, 353.97it/s]
 53%|█████▎    | 5260/10000 [00:15<00:13, 353.11it/s]
 53%|█████▎    | 5296/10000 [00:15<00:13, 353.44it/s]
 53%|█████▎    | 5332/10000 [00:15<00:13, 353.99it/s]
 54%|█████▎    | 5368/10000 [00:15<00:13, 354.32it/s]
 54%|█████▍    | 5404/10000 [00:15<00:12, 354.51it/s]
 54%|█████▍    | 5440/10000 [00:16<00:12, 354.71it/s]
 55%|█████▍    | 5476/10000 [00:16<00:12, 354.92it/s]
 55%|█████▌    | 5512/10000 [00:16<00:12, 353.19it/s]
 55%|█████▌    | 5548/10000 [00:16<00:12, 353.54it/s]
 56%|█████▌    | 5584/10000 [00:16<00:12, 354.11it/s]
 56%|█████▌    | 5620/10000 [00:16<00:12, 354.43it/s]
 57%|█████▋    | 5656/10000 [00:16<00:12, 354.42it/s]
 57%|█████▋    | 5692/10000 [00:16<00:12, 354.51it/s]
 57%|█████▋    | 5728/10000 [00:16<00:12, 354.56it/s]
 58%|█████▊    | 5764/10000 [00:16<00:11, 354.21it/s]
 58%|█████▊    | 5800/10000 [00:17<00:11, 354.38it/s]
 58%|█████▊    | 5836/10000 [00:17<00:12, 346.01it/s]
 59%|█████▊    | 5872/10000 [00:17<00:11, 348.81it/s]
 59%|█████▉    | 5908/10000 [00:17<00:11, 349.86it/s]
 59%|█████▉    | 5944/10000 [00:17<00:11, 351.79it/s]
 60%|█████▉    | 5980/10000 [00:17<00:11, 352.60it/s]
 60%|██████    | 6016/10000 [00:17<00:11, 353.33it/s]
 61%|██████    | 6052/10000 [00:17<00:11, 353.88it/s]
 61%|██████    | 6088/10000 [00:17<00:11, 354.71it/s]
 61%|██████    | 6124/10000 [00:18<00:10, 354.43it/s]
 62%|██████▏   | 6160/10000 [00:18<00:10, 354.58it/s]
 62%|██████▏   | 6196/10000 [00:18<00:10, 354.80it/s]
 62%|██████▏   | 6232/10000 [00:18<00:10, 354.70it/s]
 63%|██████▎   | 6268/10000 [00:18<00:10, 355.11it/s]
 63%|██████▎   | 6304/10000 [00:18<00:10, 355.37it/s]
 63%|██████▎   | 6340/10000 [00:18<00:10, 355.15it/s]
 64%|██████▍   | 6376/10000 [00:18<00:10, 355.21it/s]
 64%|██████▍   | 6412/10000 [00:18<00:10, 348.83it/s]
 64%|██████▍   | 6447/10000 [00:18<00:10, 343.46it/s]
 65%|██████▍   | 6483/10000 [00:19<00:10, 346.00it/s]
 65%|██████▌   | 6519/10000 [00:19<00:09, 348.45it/s]
 66%|██████▌   | 6555/10000 [00:19<00:09, 350.45it/s]
 66%|██████▌   | 6591/10000 [00:19<00:09, 351.75it/s]
 66%|██████▋   | 6627/10000 [00:19<00:09, 353.11it/s]
 67%|██████▋   | 6663/10000 [00:19<00:09, 351.52it/s]
 67%|██████▋   | 6699/10000 [00:19<00:09, 352.65it/s]
 67%|██████▋   | 6735/10000 [00:19<00:09, 352.09it/s]
 68%|██████▊   | 6771/10000 [00:19<00:09, 353.11it/s]
 68%|██████▊   | 6807/10000 [00:19<00:09, 353.45it/s]
 68%|██████▊   | 6843/10000 [00:20<00:08, 352.55it/s]
 69%|██████▉   | 6879/10000 [00:20<00:08, 353.68it/s]
 69%|██████▉   | 6915/10000 [00:20<00:08, 354.08it/s]
 70%|██████▉   | 6951/10000 [00:20<00:08, 352.94it/s]
 70%|██████▉   | 6987/10000 [00:20<00:08, 352.76it/s]
 70%|███████   | 7023/10000 [00:20<00:08, 352.62it/s]
 71%|███████   | 7059/10000 [00:20<00:08, 353.58it/s]
 71%|███████   | 7095/10000 [00:20<00:08, 348.10it/s]
 71%|███████▏  | 7131/10000 [00:20<00:08, 349.95it/s]
 72%|███████▏  | 7167/10000 [00:20<00:08, 351.72it/s]
 72%|███████▏  | 7203/10000 [00:21<00:07, 352.45it/s]
 72%|███████▏  | 7239/10000 [00:21<00:07, 353.48it/s]
 73%|███████▎  | 7275/10000 [00:21<00:07, 352.27it/s]
 73%|███████▎  | 7311/10000 [00:21<00:07, 353.43it/s]
 73%|███████▎  | 7347/10000 [00:21<00:07, 353.77it/s]
 74%|███████▍  | 7383/10000 [00:21<00:07, 353.44it/s]
 74%|███████▍  | 7419/10000 [00:21<00:07, 352.53it/s]
 75%|███████▍  | 7455/10000 [00:21<00:07, 352.67it/s]
 75%|███████▍  | 7491/10000 [00:21<00:07, 353.67it/s]
 75%|███████▌  | 7527/10000 [00:21<00:06, 354.07it/s]
 76%|███████▌  | 7563/10000 [00:22<00:06, 352.14it/s]
 76%|███████▌  | 7599/10000 [00:22<00:06, 351.56it/s]
 76%|███████▋  | 7635/10000 [00:22<00:06, 349.68it/s]
 77%|███████▋  | 7671/10000 [00:22<00:06, 351.26it/s]
 77%|███████▋  | 7707/10000 [00:22<00:06, 352.60it/s]
 77%|███████▋  | 7743/10000 [00:22<00:06, 350.69it/s]
 78%|███████▊  | 7779/10000 [00:22<00:06, 352.12it/s]
 78%|███████▊  | 7815/10000 [00:22<00:06, 352.18it/s]
 79%|███████▊  | 7851/10000 [00:22<00:06, 353.12it/s]
 79%|███████▉  | 7887/10000 [00:23<00:05, 353.52it/s]
 79%|███████▉  | 7923/10000 [00:23<00:05, 354.18it/s]
 80%|███████▉  | 7959/10000 [00:23<00:05, 353.64it/s]
 80%|███████▉  | 7995/10000 [00:23<00:05, 350.71it/s]
 80%|████████  | 8031/10000 [00:23<00:05, 351.95it/s]
 81%|████████  | 8067/10000 [00:23<00:05, 350.17it/s]
 81%|████████  | 8103/10000 [00:23<00:05, 351.73it/s]
 81%|████████▏ | 8139/10000 [00:23<00:05, 351.59it/s]
 82%|████████▏ | 8175/10000 [00:23<00:05, 352.27it/s]
 82%|████████▏ | 8211/10000 [00:23<00:05, 352.57it/s]
 82%|████████▏ | 8247/10000 [00:24<00:04, 352.58it/s]
 83%|████████▎ | 8283/10000 [00:24<00:04, 352.94it/s]
 83%|████████▎ | 8319/10000 [00:24<00:04, 352.74it/s]
 84%|████████▎ | 8355/10000 [00:24<00:04, 353.46it/s]
 84%|████████▍ | 8391/10000 [00:24<00:04, 353.89it/s]
 84%|████████▍ | 8427/10000 [00:24<00:04, 353.74it/s]
 85%|████████▍ | 8463/10000 [00:24<00:04, 353.96it/s]
 85%|████████▍ | 8499/10000 [00:24<00:04, 353.18it/s]
 85%|████████▌ | 8535/10000 [00:24<00:04, 353.35it/s]
 86%|████████▌ | 8571/10000 [00:24<00:04, 353.89it/s]
 86%|████████▌ | 8607/10000 [00:25<00:03, 354.69it/s]
 86%|████████▋ | 8643/10000 [00:25<00:03, 354.87it/s]
 87%|████████▋ | 8679/10000 [00:25<00:03, 354.67it/s]
 87%|████████▋ | 8715/10000 [00:25<00:03, 353.72it/s]
 88%|████████▊ | 8751/10000 [00:25<00:03, 353.33it/s]
 88%|████████▊ | 8787/10000 [00:25<00:03, 353.42it/s]
 88%|████████▊ | 8823/10000 [00:25<00:03, 353.71it/s]
 89%|████████▊ | 8859/10000 [00:25<00:03, 353.68it/s]
 89%|████████▉ | 8895/10000 [00:25<00:04, 271.60it/s]
 89%|████████▉ | 8926/10000 [00:26<00:04, 267.35it/s]
 90%|████████▉ | 8955/10000 [00:26<00:04, 240.53it/s]
 90%|████████▉ | 8981/10000 [00:26<00:04, 222.45it/s]
 90%|█████████ | 9013/10000 [00:26<00:04, 243.91it/s]
 90%|█████████ | 9048/10000 [00:26<00:03, 270.40it/s]
 91%|█████████ | 9084/10000 [00:26<00:03, 292.50it/s]
 91%|█████████ | 9120/10000 [00:26<00:02, 309.54it/s]
 92%|█████████▏| 9156/10000 [00:26<00:02, 321.93it/s]
 92%|█████████▏| 9192/10000 [00:27<00:02, 330.99it/s]
 92%|█████████▏| 9228/10000 [00:27<00:02, 337.82it/s]
 93%|█████████▎| 9264/10000 [00:27<00:02, 342.76it/s]
 93%|█████████▎| 9300/10000 [00:27<00:02, 345.14it/s]
 93%|█████████▎| 9335/10000 [00:27<00:01, 346.54it/s]
 94%|█████████▎| 9371/10000 [00:27<00:01, 348.66it/s]
 94%|█████████▍| 9407/10000 [00:27<00:01, 350.22it/s]
 94%|█████████▍| 9443/10000 [00:27<00:01, 350.06it/s]
 95%|█████████▍| 9479/10000 [00:27<00:01, 351.05it/s]
 95%|█████████▌| 9515/10000 [00:27<00:01, 352.07it/s]
 96%|█████████▌| 9551/10000 [00:28<00:01, 352.39it/s]
 96%|█████████▌| 9587/10000 [00:28<00:01, 352.74it/s]
 96%|█████████▌| 9623/10000 [00:28<00:01, 352.71it/s]
 97%|█████████▋| 9659/10000 [00:28<00:00, 352.98it/s]
 97%|█████████▋| 9695/10000 [00:28<00:00, 352.83it/s]
 97%|█████████▋| 9731/10000 [00:28<00:00, 353.00it/s]
 98%|█████████▊| 9767/10000 [00:28<00:00, 353.47it/s]
 98%|█████████▊| 9803/10000 [00:28<00:00, 353.35it/s]
 98%|█████████▊| 9839/10000 [00:28<00:00, 353.38it/s]
 99%|█████████▉| 9875/10000 [00:28<00:00, 352.63it/s]
 99%|█████████▉| 9911/10000 [00:29<00:00, 353.11it/s]
 99%|█████████▉| 9947/10000 [00:29<00:00, 353.43it/s]
100%|█████████▉| 9983/10000 [00:29<00:00, 353.21it/s]
100%|██████████| 10000/10000 [00:29<00:00, 341.34it/s]
The inversion result from `emcee`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
sampler: <emcee.ensemble.EnsembleSampler object>
blob_names: ['log_likelihood', 'log_prior']

Post-sampling analysis#

By default the raw sampler resulting object is attached to cofi’s inversion result.

Optionally, you can convert that into an arviz data structure to have access to a range of analysis functions. (See more details in arviz documentation).

import arviz as az

labels = ["m0", "m1", "m2","m3","m4"]

sampler = inv_result_3.sampler
az_idata = az.from_emcee(sampler, var_names=labels)
# az_idata = inv_result_3.to_arviz()      # alternatively

az_idata.get("posterior")

<xarray.DataTree 'posterior'>
Group: /posterior
    Dimensions:  (draw: 10000, chain: 32)
    Coordinates:
      * draw     (draw) int64 80kB 0 1 2 3 4 5 6 ... 9994 9995 9996 9997 9998 9999
      * chain    (chain) int64 256B 0 1 2 3 4 5 6 7 8 ... 23 24 25 26 27 28 29 30 31
    Data variables:
        m0       (draw, chain) float64 3MB 0.0001384 0.0001054 ... 1.483 1.274
        m1       (draw, chain) float64 3MB 2.934e-05 -1.438e-06 ... -3.354 -2.681
        m2       (draw, chain) float64 3MB -7.781e-05 8.234e-05 ... 1.487 1.253
        m3       (draw, chain) float64 3MB 1.866e-05 -8.024e-05 ... -0.2165 -0.191
        m4       (draw, chain) float64 3MB -3.765e-05 -0.0001327 ... 0.005969
    Attributes:
        created_at:                 2026-05-25T05:48:42.583202+00:00
        creation_library:           ArviZ
        creation_library_version:   1.1.0
        creation_library_language:  Python
        inference_library:          emcee
        inference_library_version:  3.1.6
        sample_dims:                ['draw', 'chain']


# a standard `trace` plot
labels = ["m0", "m1", "m2", "m3", "m4"]
pc = az.plot_trace_dist(az_idata, visuals={"xlabel_trace": False, "trace": {"lw": 0.5}, "dist": {"lw": 0.5}});

# add legends
n_vars = len(az_idata.posterior.data_vars)
for i in range(n_vars):
    ax1 = pc.iget_target(i, 0)  # dist (KDE) column
    ax2 = pc.iget_target(i, 1)  # trace column
    ax1.set_title(labels[i])
    ax1.set_xlabel("parameter value")
    ax1.set_ylabel("density value")
    ax2.set_title(labels[i])
    ax2.set_xlabel("number of iterations")
    ax2.set_ylabel("parameter value")

#tau = sampler.get_autocorr_time()
#print(f"autocorrelation time: {tau}")

# a Corner plot

if(True): # if we are plotting the model ensemble use this
    pm = az.plot_pair(
        az_idata.sel(draw=slice(300,None)),
        marginal=True,
        triangle="lower",
        visuals={"scatter": {"s": 2}},
    )
else: # if we wish to plot a kernel density plot then use this option
    pm = az.plot_pair(
        az_idata.sel(draw=slice(300,None)),
        marginal=True,
        triangle="lower",
        visuals={"scatter": False, "dist": {"kind": "kde"}},
    )

Now we plot the predicted curves for the posterior ensemble of solutions.

flat_samples = sampler.get_chain(discard=300, thin=30, flat=True)
inds = np.random.randint(len(flat_samples), size=100) # get a random selection from posterior ensemble

plot_data(title="Eustatic sea-level")
plt.xlim(0,maxtime)
plot_models(flat_samples[inds],color="lightgrey")
plot_model(ref_x,ref_y, "Reference model", color="darkorange")
#plt.xlim(15,20.)
#plt.ylim(-140,-100)

Expected values, credible intervals and model covariance matrix from the ensemble#

print("\n Expected value and 95% credible intervals ")
for i in range(ndim):
    mcmc = np.percentile(flat_samples[:, i], [5, 50, 95])
    print(" {} {:7.3f} [{:7.3f}, {:7.3f}]".format(labels[i],mcmc[1],mcmc[0],mcmc[2]))

Expected value and 95% credible intervals
m0   1.443 [  1.013,   1.865]
m1  -3.110 [ -3.558,  -2.666]
m2   1.412 [  1.273,   1.552]
m3  -0.209 [ -0.225,  -0.193]
m4   0.007 [  0.006,   0.007]

CMpost = np.cov(flat_samples.T)
CM_std= np.std(flat_samples,axis=0)
print('Posterior model covariance matrix\n',CMpost)
print('\n Posterior estimate of model standard deviations in each parameter')
for i in range(ndim):
    print("    {} {:7.4f}".format(labels[i],CM_std[i]))

Posterior model covariance matrix
 [[ 6.80084803e-02 -6.32698260e-02  1.67246604e-02 -1.59738238e-03
   4.72195087e-05]
 [-6.32698260e-02  7.40276118e-02 -2.20756670e-02  2.27920023e-03
  -7.05562797e-05]
 [ 1.67246604e-02 -2.20756670e-02  7.23406771e-03 -7.99827996e-04
   2.57885217e-05]
 [-1.59738238e-03  2.27920023e-03 -7.99827996e-04  9.34068185e-05
  -3.11572352e-06]
 [ 4.72195087e-05 -7.05562797e-05  2.57885217e-05 -3.11572352e-06
   1.06746135e-07]]

 Posterior estimate of model standard deviations in each parameter
    m0  0.2608
    m1  0.2721
    m2  0.0850
    m3  0.0097
    m4  0.0003

Challenge - Change the prior model bounds#

Replace the previous prior bounds to new values

The original uniform bounds had

\({\mathbf l}^T = (-10.,-10.,-10.,-10.)\), and \({\mathbf u}^T = (10.,10.,10.,10.)\).

Lets replace with

\({\mathbf l}^T = (-0.5,-10.,-10.,-10.)\), and \({\mathbf u}^T = (0.5,10.,10.,10.)\).

We have only changed the bounds of the first parameter. However since the true value of constant term was 6, these bounds are now inconsistent with the true model.

What does this do to the posterior distribution?

Start from the code template below:

m_lower_bound = <CHANGE ME>             # lower bound for uniform prior
m_upper_bound = <CHANGE ME>             # upper bound for uniform prior

def log_prior(model):    # uniform distribution
    for i in range(len(m_lower_bound)):
        if model[i] < m_lower_bound[i] or model[i] > m_upper_bound[i]: return -np.inf
    return 0.0 # model lies within bounds -> return log(1)

######## CoFI BaseProblem - update information
inv_problem.set_log_prior(log_prior)

######## CoFI Inversion - run it
inv_4 = Inversion(inv_problem, inv_options_3)
inv_result_4 = inv_4.run()

flat_samples = inv_result_4.sampler.get_chain(discard=300, thin=30, flat=True)
inds = np.random.randint(len(flat_samples), size=100) # get a random selection from posterior ensemble

print("Resulting samples with prior model lower bounds of <CHANGE ME>, upper bounds of <CHANGE ME>")
plot_data()
plot_models(flat_samples[inds])
plot_model(x, true_y, "True model", color="darkorange")

# Copy the template above, Replace <CHANGE ME> with your answer

#@title Solution

m_lower_bound = np.array([-1.0,-10,-10,-10])             # lower bound for uniform prior
m_upper_bound = np.array([1.0,10,10,10])                 # upper bound for uniform prior

def log_prior(model):    # uniform distribution
    for i in range(len(m_lower_bound)):
        if model[i] < m_lower_bound[i] or model[i] > m_upper_bound[i]: return -np.inf
    return 0.0 # model lies within bounds -> return log(1)

######## CoFI BaseProblem - update information
inv_problem.set_log_prior(log_prior)

######## CoFI Inversion - run it
inv_4 = Inversion(inv_problem, inv_options_3)
inv_result_4 = inv_4.run()

flat_samples = inv_result_4.sampler.get_chain(discard=300, thin=30, flat=True)
inds = np.random.randint(len(flat_samples), size=100) # get a random selection from posterior ensemble

print("Resulting samples with prior model lower bounds of [-1,-10,-10,-10], upper bounds of [2,10,10,10]")
plot_data()
plot_models(flat_samples[inds],color="lightgrey")
plot_model(ref_x, ref_y, "Reference model", color="darkorange")

  0%|          | 0/10000 [00:00<?, ?it/s]
  0%|          | 35/10000 [00:00<00:28, 349.02it/s]
  1%|          | 70/10000 [00:00<00:29, 336.96it/s]
  1%|          | 106/10000 [00:00<00:28, 345.99it/s]
  1%|▏         | 142/10000 [00:00<00:28, 349.84it/s]
  2%|▏         | 178/10000 [00:00<00:27, 352.02it/s]
  2%|▏         | 214/10000 [00:00<00:27, 351.75it/s]
  3%|▎         | 251/10000 [00:00<00:27, 356.77it/s]
  3%|▎         | 290/10000 [00:00<00:26, 364.33it/s]
  3%|▎         | 329/10000 [00:00<00:26, 369.71it/s]
  4%|▎         | 368/10000 [00:01<00:25, 374.44it/s]
  4%|▍         | 407/10000 [00:01<00:25, 378.87it/s]
  4%|▍         | 447/10000 [00:01<00:24, 384.26it/s]
  5%|▍         | 486/10000 [00:01<00:25, 374.29it/s]
  5%|▌         | 524/10000 [00:01<00:25, 371.04it/s]
  6%|▌         | 564/10000 [00:01<00:24, 378.66it/s]
  6%|▌         | 603/10000 [00:01<00:24, 380.00it/s]
  6%|▋         | 642/10000 [00:01<00:24, 382.36it/s]
  7%|▋         | 681/10000 [00:01<00:24, 383.32it/s]
  7%|▋         | 720/10000 [00:01<00:24, 383.26it/s]
  8%|▊         | 759/10000 [00:02<00:23, 385.14it/s]
  8%|▊         | 798/10000 [00:02<00:23, 386.13it/s]
  8%|▊         | 837/10000 [00:02<00:23, 384.33it/s]
  9%|▉         | 876/10000 [00:02<00:23, 384.86it/s]
  9%|▉         | 916/10000 [00:02<00:23, 387.09it/s]
 10%|▉         | 956/10000 [00:02<00:23, 389.74it/s]
 10%|▉         | 996/10000 [00:02<00:22, 392.64it/s]
 10%|█         | 1036/10000 [00:02<00:22, 391.27it/s]
 11%|█         | 1076/10000 [00:02<00:22, 391.61it/s]
 11%|█         | 1116/10000 [00:02<00:22, 391.61it/s]
 12%|█▏        | 1156/10000 [00:03<00:22, 388.61it/s]
 12%|█▏        | 1195/10000 [00:03<00:22, 387.58it/s]
 12%|█▏        | 1234/10000 [00:03<00:22, 387.27it/s]
 13%|█▎        | 1273/10000 [00:03<00:22, 387.98it/s]
 13%|█▎        | 1312/10000 [00:03<00:22, 385.16it/s]
 14%|█▎        | 1352/10000 [00:03<00:22, 387.45it/s]
 14%|█▍        | 1391/10000 [00:03<00:22, 387.86it/s]
 14%|█▍        | 1430/10000 [00:03<00:22, 385.00it/s]
 15%|█▍        | 1469/10000 [00:03<00:22, 383.59it/s]
 15%|█▌        | 1508/10000 [00:03<00:22, 383.14it/s]
 15%|█▌        | 1547/10000 [00:04<00:23, 363.85it/s]
 16%|█▌        | 1584/10000 [00:04<00:23, 356.84it/s]
 16%|█▌        | 1620/10000 [00:04<00:28, 292.35it/s]
 17%|█▋        | 1652/10000 [00:04<00:31, 263.31it/s]
 17%|█▋        | 1688/10000 [00:04<00:29, 286.31it/s]
 17%|█▋        | 1727/10000 [00:04<00:26, 311.08it/s]
 18%|█▊        | 1767/10000 [00:04<00:24, 332.93it/s]
 18%|█▊        | 1806/10000 [00:04<00:23, 347.27it/s]
 18%|█▊        | 1846/10000 [00:05<00:22, 361.62it/s]
 19%|█▉        | 1886/10000 [00:05<00:21, 370.01it/s]
 19%|█▉        | 1926/10000 [00:05<00:21, 377.18it/s]
 20%|█▉        | 1966/10000 [00:05<00:21, 382.40it/s]
 20%|██        | 2005/10000 [00:05<00:21, 366.45it/s]
 20%|██        | 2045/10000 [00:05<00:21, 375.10it/s]
 21%|██        | 2085/10000 [00:05<00:20, 379.99it/s]
 21%|██        | 2124/10000 [00:05<00:20, 382.57it/s]
 22%|██▏       | 2163/10000 [00:05<00:20, 384.74it/s]
 22%|██▏       | 2203/10000 [00:05<00:20, 389.03it/s]
 22%|██▏       | 2242/10000 [00:06<00:20, 387.60it/s]
 23%|██▎       | 2282/10000 [00:06<00:19, 389.72it/s]
 23%|██▎       | 2322/10000 [00:06<00:19, 391.02it/s]
 24%|██▎       | 2362/10000 [00:06<00:19, 391.22it/s]
 24%|██▍       | 2402/10000 [00:06<00:19, 391.15it/s]
 24%|██▍       | 2442/10000 [00:06<00:19, 389.28it/s]
 25%|██▍       | 2481/10000 [00:06<00:19, 386.28it/s]
 25%|██▌       | 2520/10000 [00:06<00:19, 386.68it/s]
 26%|██▌       | 2559/10000 [00:06<00:19, 385.05it/s]
 26%|██▌       | 2599/10000 [00:06<00:19, 387.06it/s]
 26%|██▋       | 2638/10000 [00:07<00:19, 381.24it/s]
 27%|██▋       | 2677/10000 [00:07<00:19, 383.16it/s]
 27%|██▋       | 2716/10000 [00:07<00:19, 381.51it/s]
 28%|██▊       | 2755/10000 [00:07<00:18, 381.71it/s]
 28%|██▊       | 2795/10000 [00:07<00:18, 386.25it/s]
 28%|██▊       | 2835/10000 [00:07<00:18, 387.82it/s]
 29%|██▊       | 2874/10000 [00:07<00:18, 387.74it/s]
 29%|██▉       | 2914/10000 [00:07<00:18, 389.49it/s]
 30%|██▉       | 2953/10000 [00:07<00:18, 389.17it/s]
 30%|██▉       | 2993/10000 [00:08<00:17, 390.57it/s]
 30%|███       | 3033/10000 [00:08<00:17, 387.69it/s]
 31%|███       | 3073/10000 [00:08<00:17, 388.71it/s]
 31%|███       | 3112/10000 [00:08<00:17, 388.04it/s]
 32%|███▏      | 3151/10000 [00:08<00:17, 388.23it/s]
 32%|███▏      | 3191/10000 [00:08<00:17, 389.74it/s]
 32%|███▏      | 3230/10000 [00:08<00:17, 388.91it/s]
 33%|███▎      | 3270/10000 [00:08<00:17, 389.32it/s]
 33%|███▎      | 3309/10000 [00:08<00:17, 389.45it/s]
 33%|███▎      | 3349/10000 [00:08<00:17, 390.25it/s]
 34%|███▍      | 3389/10000 [00:09<00:17, 388.16it/s]
 34%|███▍      | 3428/10000 [00:09<00:17, 384.43it/s]
 35%|███▍      | 3467/10000 [00:09<00:20, 315.24it/s]
 35%|███▌      | 3503/10000 [00:09<00:20, 324.55it/s]
 35%|███▌      | 3538/10000 [00:09<00:23, 272.21it/s]
 36%|███▌      | 3568/10000 [00:09<00:23, 270.93it/s]
 36%|███▌      | 3607/10000 [00:09<00:21, 300.34it/s]
 36%|███▋      | 3647/10000 [00:09<00:19, 325.47it/s]
 37%|███▋      | 3687/10000 [00:10<00:18, 343.92it/s]
 37%|███▋      | 3727/10000 [00:10<00:17, 357.56it/s]
 38%|███▊      | 3766/10000 [00:10<00:17, 366.53it/s]
 38%|███▊      | 3806/10000 [00:10<00:16, 373.96it/s]
 38%|███▊      | 3845/10000 [00:10<00:16, 377.11it/s]
 39%|███▉      | 3884/10000 [00:10<00:16, 379.94it/s]
 39%|███▉      | 3923/10000 [00:10<00:15, 382.26it/s]
 40%|███▉      | 3962/10000 [00:10<00:15, 382.04it/s]
 40%|████      | 4001/10000 [00:10<00:15, 382.76it/s]
 40%|████      | 4041/10000 [00:10<00:15, 386.42it/s]
 41%|████      | 4081/10000 [00:11<00:15, 389.70it/s]
 41%|████      | 4121/10000 [00:11<00:15, 391.20it/s]
 42%|████▏     | 4161/10000 [00:11<00:14, 391.60it/s]
 42%|████▏     | 4201/10000 [00:11<00:14, 391.27it/s]
 42%|████▏     | 4241/10000 [00:11<00:14, 393.48it/s]
 43%|████▎     | 4281/10000 [00:11<00:14, 390.06it/s]
 43%|████▎     | 4321/10000 [00:11<00:14, 390.55it/s]
 44%|████▎     | 4361/10000 [00:11<00:14, 385.88it/s]
 44%|████▍     | 4400/10000 [00:11<00:14, 386.20it/s]
 44%|████▍     | 4440/10000 [00:11<00:14, 387.78it/s]
 45%|████▍     | 4479/10000 [00:12<00:14, 386.61it/s]
 45%|████▌     | 4519/10000 [00:12<00:14, 389.66it/s]
 46%|████▌     | 4559/10000 [00:12<00:13, 389.76it/s]
 46%|████▌     | 4598/10000 [00:12<00:13, 389.01it/s]
 46%|████▋     | 4637/10000 [00:12<00:13, 386.89it/s]
 47%|████▋     | 4676/10000 [00:12<00:13, 386.85it/s]
 47%|████▋     | 4715/10000 [00:12<00:13, 386.51it/s]
 48%|████▊     | 4754/10000 [00:12<00:13, 385.52it/s]
 48%|████▊     | 4794/10000 [00:12<00:13, 387.40it/s]
 48%|████▊     | 4834/10000 [00:12<00:13, 389.47it/s]
 49%|████▊     | 4873/10000 [00:13<00:13, 388.28it/s]
 49%|████▉     | 4912/10000 [00:13<00:13, 387.08it/s]
 50%|████▉     | 4951/10000 [00:13<00:13, 385.00it/s]
 50%|████▉     | 4990/10000 [00:13<00:13, 358.42it/s]
 50%|█████     | 5029/10000 [00:13<00:13, 365.49it/s]
 51%|█████     | 5069/10000 [00:13<00:13, 372.79it/s]
 51%|█████     | 5109/10000 [00:13<00:12, 378.59it/s]
 51%|█████▏    | 5148/10000 [00:13<00:12, 381.00it/s]
 52%|█████▏    | 5187/10000 [00:13<00:12, 381.37it/s]
 52%|█████▏    | 5226/10000 [00:13<00:12, 383.52it/s]
 53%|█████▎    | 5265/10000 [00:14<00:12, 385.35it/s]
 53%|█████▎    | 5304/10000 [00:14<00:12, 384.74it/s]
 53%|█████▎    | 5343/10000 [00:14<00:13, 351.40it/s]
 54%|█████▍    | 5379/10000 [00:14<00:14, 310.11it/s]
 54%|█████▍    | 5418/10000 [00:14<00:13, 329.94it/s]
 55%|█████▍    | 5457/10000 [00:14<00:13, 343.71it/s]
 55%|█████▍    | 5496/10000 [00:14<00:12, 354.57it/s]
 55%|█████▌    | 5535/10000 [00:14<00:12, 361.96it/s]
 56%|█████▌    | 5574/10000 [00:14<00:11, 369.79it/s]
 56%|█████▌    | 5613/10000 [00:15<00:11, 375.63it/s]
 57%|█████▋    | 5652/10000 [00:15<00:11, 378.18it/s]
 57%|█████▋    | 5691/10000 [00:15<00:11, 381.44it/s]
 57%|█████▋    | 5731/10000 [00:15<00:11, 385.48it/s]
 58%|█████▊    | 5770/10000 [00:15<00:10, 386.70it/s]
 58%|█████▊    | 5809/10000 [00:15<00:10, 387.47it/s]
 58%|█████▊    | 5849/10000 [00:15<00:10, 388.80it/s]
 59%|█████▉    | 5888/10000 [00:15<00:10, 386.91it/s]
 59%|█████▉    | 5927/10000 [00:15<00:10, 387.78it/s]
 60%|█████▉    | 5966/10000 [00:16<00:10, 388.09it/s]
 60%|██████    | 6005/10000 [00:16<00:10, 386.13it/s]
 60%|██████    | 6044/10000 [00:16<00:10, 386.17it/s]
 61%|██████    | 6083/10000 [00:16<00:10, 381.88it/s]
 61%|██████    | 6123/10000 [00:16<00:10, 386.73it/s]
 62%|██████▏   | 6163/10000 [00:16<00:09, 390.06it/s]
 62%|██████▏   | 6203/10000 [00:16<00:09, 389.82it/s]
 62%|██████▏   | 6243/10000 [00:16<00:09, 391.90it/s]
 63%|██████▎   | 6283/10000 [00:16<00:09, 393.26it/s]
 63%|██████▎   | 6323/10000 [00:16<00:09, 394.55it/s]
 64%|██████▎   | 6363/10000 [00:17<00:09, 393.58it/s]
 64%|██████▍   | 6403/10000 [00:17<00:09, 390.59it/s]
 64%|██████▍   | 6443/10000 [00:17<00:09, 390.93it/s]
 65%|██████▍   | 6483/10000 [00:17<00:08, 390.80it/s]
 65%|██████▌   | 6523/10000 [00:17<00:08, 390.28it/s]
 66%|██████▌   | 6563/10000 [00:17<00:08, 392.47it/s]
 66%|██████▌   | 6603/10000 [00:17<00:08, 390.96it/s]
 66%|██████▋   | 6643/10000 [00:17<00:08, 390.68it/s]
 67%|██████▋   | 6683/10000 [00:17<00:08, 390.77it/s]
 67%|██████▋   | 6723/10000 [00:17<00:08, 389.65it/s]
 68%|██████▊   | 6762/10000 [00:18<00:08, 388.85it/s]
 68%|██████▊   | 6802/10000 [00:18<00:08, 389.86it/s]
 68%|██████▊   | 6841/10000 [00:18<00:08, 386.76it/s]
 69%|██████▉   | 6880/10000 [00:18<00:08, 385.68it/s]
 69%|██████▉   | 6919/10000 [00:18<00:07, 386.27it/s]
 70%|██████▉   | 6958/10000 [00:18<00:07, 383.81it/s]
 70%|██████▉   | 6997/10000 [00:18<00:07, 383.81it/s]
 70%|███████   | 7036/10000 [00:18<00:07, 384.36it/s]
 71%|███████   | 7075/10000 [00:18<00:07, 383.55it/s]
 71%|███████   | 7114/10000 [00:18<00:07, 383.86it/s]
 72%|███████▏  | 7153/10000 [00:19<00:07, 382.11it/s]
 72%|███████▏  | 7192/10000 [00:19<00:07, 383.96it/s]
 72%|███████▏  | 7231/10000 [00:19<00:07, 384.60it/s]
 73%|███████▎  | 7270/10000 [00:19<00:07, 371.66it/s]
 73%|███████▎  | 7308/10000 [00:19<00:07, 353.02it/s]
 73%|███████▎  | 7344/10000 [00:19<00:08, 311.26it/s]
 74%|███████▍  | 7383/10000 [00:19<00:07, 331.23it/s]
 74%|███████▍  | 7418/10000 [00:19<00:10, 246.41it/s]
 75%|███████▍  | 7452/10000 [00:20<00:09, 266.34it/s]
 75%|███████▍  | 7490/10000 [00:20<00:08, 293.50it/s]
 75%|███████▌  | 7529/10000 [00:20<00:07, 317.32it/s]
 76%|███████▌  | 7568/10000 [00:20<00:07, 335.96it/s]
 76%|███████▌  | 7606/10000 [00:20<00:06, 347.52it/s]
 76%|███████▋  | 7645/10000 [00:20<00:06, 357.11it/s]
 77%|███████▋  | 7685/10000 [00:20<00:06, 367.30it/s]
 77%|███████▋  | 7725/10000 [00:20<00:06, 375.23it/s]
 78%|███████▊  | 7765/10000 [00:20<00:05, 380.95it/s]
 78%|███████▊  | 7805/10000 [00:20<00:05, 383.92it/s]
 78%|███████▊  | 7844/10000 [00:21<00:05, 385.62it/s]
 79%|███████▉  | 7884/10000 [00:21<00:05, 387.08it/s]
 79%|███████▉  | 7923/10000 [00:21<00:05, 387.18it/s]
 80%|███████▉  | 7962/10000 [00:21<00:05, 385.88it/s]
 80%|████████  | 8001/10000 [00:21<00:05, 387.04it/s]
 80%|████████  | 8040/10000 [00:21<00:05, 385.71it/s]
 81%|████████  | 8079/10000 [00:21<00:04, 384.83it/s]
 81%|████████  | 8119/10000 [00:21<00:04, 387.72it/s]
 82%|████████▏ | 8158/10000 [00:21<00:04, 388.11it/s]
 82%|████████▏ | 8198/10000 [00:22<00:04, 389.11it/s]
 82%|████████▏ | 8238/10000 [00:22<00:04, 391.51it/s]
 83%|████████▎ | 8278/10000 [00:22<00:04, 391.59it/s]
 83%|████████▎ | 8318/10000 [00:22<00:04, 392.90it/s]
 84%|████████▎ | 8358/10000 [00:22<00:04, 393.89it/s]
 84%|████████▍ | 8398/10000 [00:22<00:04, 393.80it/s]
 84%|████████▍ | 8438/10000 [00:22<00:03, 393.07it/s]
 85%|████████▍ | 8478/10000 [00:22<00:03, 392.04it/s]
 85%|████████▌ | 8518/10000 [00:22<00:03, 393.20it/s]
 86%|████████▌ | 8558/10000 [00:22<00:03, 394.21it/s]
 86%|████████▌ | 8598/10000 [00:23<00:03, 392.89it/s]
 86%|████████▋ | 8638/10000 [00:23<00:03, 390.52it/s]
 87%|████████▋ | 8678/10000 [00:23<00:03, 390.78it/s]
 87%|████████▋ | 8718/10000 [00:23<00:03, 389.38it/s]
 88%|████████▊ | 8758/10000 [00:23<00:03, 390.98it/s]
 88%|████████▊ | 8798/10000 [00:23<00:03, 392.17it/s]
 88%|████████▊ | 8838/10000 [00:23<00:02, 390.62it/s]
 89%|████████▉ | 8878/10000 [00:23<00:02, 391.36it/s]
 89%|████████▉ | 8918/10000 [00:23<00:02, 389.76it/s]
 90%|████████▉ | 8957/10000 [00:23<00:02, 388.91it/s]
 90%|████████▉ | 8996/10000 [00:24<00:02, 387.83it/s]
 90%|█████████ | 9035/10000 [00:24<00:02, 386.77it/s]
 91%|█████████ | 9075/10000 [00:24<00:02, 388.35it/s]
 91%|█████████ | 9115/10000 [00:24<00:02, 391.04it/s]
 92%|█████████▏| 9155/10000 [00:24<00:02, 388.56it/s]
 92%|█████████▏| 9194/10000 [00:24<00:02, 364.12it/s]
 92%|█████████▏| 9231/10000 [00:24<00:02, 317.51it/s]
 93%|█████████▎| 9264/10000 [00:24<00:02, 310.68it/s]
 93%|█████████▎| 9296/10000 [00:24<00:02, 273.49it/s]
 93%|█████████▎| 9327/10000 [00:25<00:02, 282.41it/s]
 94%|█████████▎| 9366/10000 [00:25<00:02, 309.32it/s]
 94%|█████████▍| 9405/10000 [00:25<00:01, 330.32it/s]
 94%|█████████▍| 9444/10000 [00:25<00:01, 345.71it/s]
 95%|█████████▍| 9483/10000 [00:25<00:01, 358.24it/s]
 95%|█████████▌| 9523/10000 [00:25<00:01, 369.86it/s]
 96%|█████████▌| 9563/10000 [00:25<00:01, 376.99it/s]
 96%|█████████▌| 9602/10000 [00:25<00:01, 378.59it/s]
 96%|█████████▋| 9641/10000 [00:25<00:00, 379.95it/s]
 97%|█████████▋| 9681/10000 [00:26<00:00, 383.27it/s]
 97%|█████████▋| 9722/10000 [00:26<00:00, 388.23it/s]
 98%|█████████▊| 9762/10000 [00:26<00:00, 391.05it/s]
 98%|█████████▊| 9802/10000 [00:26<00:00, 393.53it/s]
 98%|█████████▊| 9842/10000 [00:26<00:00, 389.95it/s]
 99%|█████████▉| 9882/10000 [00:26<00:00, 388.90it/s]
 99%|█████████▉| 9922/10000 [00:26<00:00, 390.88it/s]
100%|█████████▉| 9962/10000 [00:26<00:00, 391.14it/s]
100%|██████████| 10000/10000 [00:26<00:00, 372.83it/s]
Resulting samples with prior model lower bounds of [-1,-10,-10,-10], upper bounds of [2,10,10,10]

Is there much change to the posterior distribution?

Trans-dimensional partition modelling#

A separate CoFI example notebook that performs trans-dimensional partition modelling on the same dataset with is also available.

This is called Partition_modelling_sealevel_bayesbay.ipynb and is available in the directory examples/partition_modelling

Watermark#

watermark_list = ["cofi", "numpy", "scipy", "matplotlib", "emcee", "arviz"]
for pkg in watermark_list:
    pkg_var = __import__(pkg)
    print(pkg, getattr(pkg_var, "__version__"))

cofi 0.2.11+71.gb28b5b0
numpy 2.2.6
scipy 1.17.1
matplotlib 3.10.9
emcee 3.1.6
arviz 1.1.0

sphinx_gallery_thumbnail_number = -1

Total running time of the script: (1 minutes 5.615 seconds)

Gallery generated by Sphinx-Gallery

Contents