Note

Go to the end to download the full example code

Linear regression#

What we do in this notebook#

Here we demonstrate use of CoFI on a simple linear regression problem, where we fit a polynomial function to data, using three different algorithms:

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.

Learning outcomes#

A demonstration of running CoFI for a class of parameter fitting problem. Example of a CoFI template.
A demonstration of how CoFI may be used to experiment with different inference approaches under a common interface.
A demonstration of CoFI’s expandability in that it may be used with pre-set, or user defined, misfits, likelihood or priors.

# Environment setup (uncomment code below)

# !pip install -U cofi

Linear regression#

Lets start with some (x,y) data.

import numpy as np
import matplotlib.pyplot as plt

# here is some (x,y) data
data_x = np.array([1.1530612244897958, -0.07142857142857162, -1.7857142857142858,
                1.6428571428571423, -2.642857142857143, -1.0510204081632653,
                1.1530612244897958, -1.295918367346939, -0.806122448979592,
                -2.2755102040816326, -2.2755102040816326, -0.6836734693877551,
                0.7857142857142856, 1.2755102040816322, -0.6836734693877551,
                -3.2551020408163267, -0.9285714285714288, -3.377551020408163,
                -0.6836734693877551, 1.7653061224489797])

data_y = np.array([-7.550931153863841, -6.060810406314714, 3.080063056254076,
                -4.499764131508964, 2.9462042659962333, -0.4645899453212615,
                -7.43068837808917, 1.6273774547833582, -0.05922697815443567,
                3.8462283231266903, 3.425851020301113, -0.05359797104829345,
                -10.235538857712598, -5.929113775071286, -1.1871766078924957,
                -4.124258811692425, 0.6969191559961637, -4.454022624935177,
                -2.352842192972056, -4.25145590011172])
sigma = 1   # estimation on the data noise

And now lets plot the data.

def plot_data(sigma=None):
    if(sigma is None):
        plt.scatter(data_x, data_y, color="lightcoral", label="observed data")
    else:
        plt.errorbar(data_x, data_y, yerr=sigma, fmt='.',color="lightcoral",ecolor='lightgrey',ms=10)
plot_data()

Problem description#

To begin with, we will work with polynomial curves,

\[y(x) = \sum_{j=0}^M m_j x^j\,.\]

Here, \(M\) is the ‘order’ of the polynomial: if \(M=1\) we have a straight line with 2 parameters, if \(M=2\) it will be a quadratic with 3 parameters, and so on. The \(m_j, (j=0,\dots M)\) are the ‘model coefficients’ that we seek to constrain from the data.

For this class of problem the forward operator takes the following form:

\[\begin{split}\left(\begin{array}{c}y_0\\y_1\\\vdots\\y_N\end{array}\right) = \left(\begin{array}{ccc}1&x_0&x_0^2&x_0^3\\1&x_1&x_1^2&x_1^3\\\vdots&\vdots&\vdots\\1&x_N&x_N^2&x_N^3\end{array}\right)\left(\begin{array}{c}m_0\\m_1\\m_2\\m_3\end{array}\right)\end{split}\]

This clearly has the required general form, \(\mathbf{d} =G{\mathbf m}\).

where:

\(\textbf{d}\) is the vector of data values, (\(y_0,y_1,\dots,y_N\));
\(\textbf{m}\) is the vector of model parameters, (\(m_0,m_1,m_2\));
\(G\) is the basis matrix (or design matrix) of this linear regression problem (also called the Jacobian matrix for this linear problem).

We have a set of noisy data values, \(y_i (i=0,\dots,N)\), measured at known locations, \(x_i (i=0,\dots,N)\), and wish to find the best fit degree 3 polynomial.

The function that generated our data is : \(y=-6-5x+2x^2+x^3\), and we have added Gaussian random noise, \({\cal N}(0,\sigma^2)\), with \(\sigma=1.0\).

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

nparams = 4 # 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=sigma, ndata=len(data_x)):
    return 1/sigma**2 * np.identity(ndata)

Define the true model for later.

# True model for plotting
x = np.linspace(-3.5,2.5)              # x values to plot
true_model = np.array([-6, -5, 2, 1])  # we know it for this case which will be useful later for comparison.

true_y = jacobian(x,4).dot(true_model) # y values for true curve

Now lets plot the data with the curve from the true polynomial coefficients.

# Some plotting utilities
def plot_model(x,y, label, color=None):
    #x = np.linspace(-3.5,2.5)
    #y = jacobian(x).dot(model)
    plt.plot(x, y, color=color or "green", label=label)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.legend()

def plot_models(models, label="Posterior samples", color="seagreen", alpha=0.1):
    x = np.linspace(-3.5,2.5)
    G = jacobian(x)
    plt.plot(x, G.dot(models[0]), color=color, label=label, alpha=alpha)
    for m in models:
        plt.plot(x, G.dot(m), color=color, alpha=alpha)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.legend()

plot_data(sigma=sigma)
plot_model(x,true_y, "true 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.

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:
{'sampling', 'optimization', 'matrix solvers'}

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"
    ],
    "matrix solvers": [
        "scipy.linalg.lstsq",
        "cofi.simple_newton"
    ],
    "sampling": [
        "emcee",
        "bayesbay",
        "neighpy"
    ]
}

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`: [-5.71964359 -5.10903808  1.82553662  0.97472374]

============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [-5.71964359 -5.10903808  1.82553662  0.97472374]
sum_of_squared_residuals: []
effective_rank: 4
singular_values: [3765.51775745   69.19268194   16.27124488    3.85437889]
model_covariance: [[ 0.19027447  0.05812534 -0.08168411 -0.02550866]
 [ 0.05812534  0.08673796 -0.03312809 -0.01812686]
 [-0.08168411 -0.03312809  0.05184851  0.01704165]
 [-0.02550866 -0.01812686  0.01704165  0.00676031]]

Lets plot the solution.

plot_data()
plot_model(x,jacobian(x).dot(inv_result.model), "linear system solver", color="seagreen")
plot_model(x,true_y, "true 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_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`: [-5.71967431 -5.10913992  1.82556456  0.9747426 ]

============================
Summary for inversion result
============================
SUCCESS
----------------------------
fun: 14.961508008942793
nit: 193
nfev: 330
status: 0
message: Optimization terminated successfully.
final_simplex: (array([[-5.71967431, -5.10913992,  1.82556456,  0.9747426 ],
       [-5.71958302, -5.10907158,  1.8255083 ,  0.97472628],
       [-5.71969118, -5.10911404,  1.82556388,  0.97474474],
       [-5.7197282 , -5.10917942,  1.82554925,  0.97474097],
       [-5.71960767, -5.10913354,  1.82551338,  0.97473478]]), array([14.96150801, 14.96150804, 14.96150805, 14.9615082 , 14.96150821]))
model: [-5.71967431 -5.10913992  1.82556456  0.9747426 ]

plot_data()
plot_model(x,jacobian(x).dot(inv_result_2.model), "optimization solution", color="cornflowerblue")
plot_model(x,true_y, "true model", color="darkorange")

Challenge: Change the polynomial degree#

Try and replace the 3rd order polynomial with a 1st order polynomial (i.e. \(M=1\)) 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(x,jacobian(x,n=<CHANGE ME>).dot(inv_result.model), "optimization solution", color="cornflowerblue")
plot_model(x,true_y, "true 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=2))
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 = 2 ")
plot_data()
plot_model(x,jacobian(x,n=2).dot(inv_result.model), "optimization solution", color="cornflowerblue")
plot_model(x,true_y, "true model", color="darkorange")

Inferred curve with n = 2

Where to next?#

Linear regression with Eustatic Sea-level data - link to notebook

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.7
numpy 1.24.4
scipy 1.12.0
matplotlib 3.8.3
emcee 3.1.4
arviz 0.17.0

sphinx_gallery_thumbnail_number = -1

Total running time of the script: (0 minutes 29.038 seconds)

Gallery generated by Sphinx-Gallery

Linear regression

Contents

Linear regression#

What we do in this notebook#

Learning outcomes#

Linear regression#

Problem description#

The structure of CoFI#

1. Linear system solver#

Step 1. Define CoFI `BaseProblem`#

Step 2. Define CoFI `InversionOptions`#

Step 3. Define CoFI `Inversion` and run#

2. Optimizer#

Challenge: Change the polynomial degree#

3. Bayesian sampling#

Likelihood#

Prior#

Bayesian sampling#

Starting points for random walkers#

Add the information and run with CoFI#

Post-sampling analysis#

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

Challenge: Change the prior model bounds#

Challenge: Change the data uncertainty#

Challenge: Change the number of walkers / steps in the McMC algorithm (optional)#

Where to next?#

Watermark#

Linear regression

Contents

Linear regression#

What we do in this notebook#

Learning outcomes#

Linear regression#

Problem description#

The structure of CoFI#

1. Linear system solver#

Step 1. Define CoFI BaseProblem#

Step 2. Define CoFI InversionOptions#

Step 3. Define CoFI Inversion and run#

2. Optimizer#

Challenge: Change the polynomial degree#

3. Bayesian sampling#

Likelihood#

Prior#

Bayesian sampling#

Starting points for random walkers#

Add the information and run with CoFI#

Post-sampling analysis#

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

Challenge: Change the prior model bounds#

Challenge: Change the data uncertainty#

Challenge: Change the number of walkers / steps in the McMC algorithm (optional)#

Where to next?#

Watermark#

Step 1. Define CoFI `BaseProblem`#

Step 2. Define CoFI `InversionOptions`#

Step 3. Define CoFI `Inversion` and run#