1 - Linear regression#

Open In Colab


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()
1 linear regression

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")
1 linear regression

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:
{'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"
    ],
    "matrix solvers": [
        "scipy.linalg.lstsq",
        "cofi.simple_newton"
    ],
    "sampling": [
        "emcee"
    ]
}
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: `scipy.linalg.lstsq (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")
1 linear regression

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.961508008942795
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")
1 linear regression

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?

Upload to Jamboard 1

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")
1 linear regression
Inferred curve with n = 2

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\).

sigma = 1.0                                     # common noise standard deviation
Cdinv = np.eye(len(data_y))/(sigma**2)      # 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, initial_state=walkers_start, progress=True)

######## 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]
  1%|1         | 119/10000 [00:00<00:08, 1184.89it/s]
  2%|2         | 241/10000 [00:00<00:08, 1199.70it/s]
  4%|3         | 363/10000 [00:00<00:07, 1206.12it/s]
  5%|4         | 485/10000 [00:00<00:07, 1209.09it/s]
  6%|6         | 607/10000 [00:00<00:07, 1210.57it/s]
  7%|7         | 729/10000 [00:00<00:07, 1211.48it/s]
  9%|8         | 851/10000 [00:00<00:07, 1212.32it/s]
 10%|9         | 973/10000 [00:00<00:07, 1210.65it/s]
 11%|#         | 1095/10000 [00:00<00:07, 1209.34it/s]
 12%|#2        | 1217/10000 [00:01<00:07, 1210.89it/s]
 13%|#3        | 1339/10000 [00:01<00:07, 1212.02it/s]
 15%|#4        | 1461/10000 [00:01<00:07, 1213.07it/s]
 16%|#5        | 1583/10000 [00:01<00:06, 1213.55it/s]
 17%|#7        | 1705/10000 [00:01<00:06, 1214.41it/s]
 18%|#8        | 1827/10000 [00:01<00:06, 1214.63it/s]
 19%|#9        | 1949/10000 [00:01<00:06, 1214.90it/s]
 21%|##        | 2071/10000 [00:01<00:06, 1214.88it/s]
 22%|##1       | 2193/10000 [00:01<00:06, 1215.37it/s]
 23%|##3       | 2315/10000 [00:01<00:06, 1215.12it/s]
 24%|##4       | 2437/10000 [00:02<00:06, 1215.09it/s]
 26%|##5       | 2559/10000 [00:02<00:06, 1214.87it/s]
 27%|##6       | 2681/10000 [00:02<00:06, 1215.42it/s]
 28%|##8       | 2803/10000 [00:02<00:05, 1215.43it/s]
 29%|##9       | 2925/10000 [00:02<00:05, 1215.07it/s]
 30%|###       | 3047/10000 [00:02<00:05, 1215.30it/s]
 32%|###1      | 3169/10000 [00:02<00:05, 1215.72it/s]
 33%|###2      | 3291/10000 [00:02<00:05, 1215.53it/s]
 34%|###4      | 3413/10000 [00:02<00:05, 1215.59it/s]
 35%|###5      | 3535/10000 [00:02<00:05, 1215.03it/s]
 37%|###6      | 3657/10000 [00:03<00:05, 1215.26it/s]
 38%|###7      | 3779/10000 [00:03<00:05, 1215.50it/s]
 39%|###9      | 3901/10000 [00:03<00:05, 1216.04it/s]
 40%|####      | 4023/10000 [00:03<00:04, 1215.59it/s]
 41%|####1     | 4145/10000 [00:03<00:04, 1214.76it/s]
 43%|####2     | 4267/10000 [00:03<00:04, 1214.03it/s]
 44%|####3     | 4389/10000 [00:03<00:04, 1213.55it/s]
 45%|####5     | 4511/10000 [00:03<00:04, 1213.90it/s]
 46%|####6     | 4633/10000 [00:03<00:04, 1214.65it/s]
 48%|####7     | 4755/10000 [00:03<00:04, 1213.12it/s]
 49%|####8     | 4877/10000 [00:04<00:04, 1212.61it/s]
 50%|####9     | 4999/10000 [00:04<00:04, 1211.53it/s]
 51%|#####1    | 5121/10000 [00:04<00:04, 1211.75it/s]
 52%|#####2    | 5243/10000 [00:04<00:03, 1211.44it/s]
 54%|#####3    | 5365/10000 [00:04<00:03, 1212.04it/s]
 55%|#####4    | 5487/10000 [00:04<00:03, 1212.33it/s]
 56%|#####6    | 5609/10000 [00:04<00:03, 1212.29it/s]
 57%|#####7    | 5731/10000 [00:04<00:03, 1206.94it/s]
 59%|#####8    | 5853/10000 [00:04<00:03, 1208.38it/s]
 60%|#####9    | 5974/10000 [00:04<00:03, 1208.35it/s]
 61%|######    | 6096/10000 [00:05<00:03, 1210.29it/s]
 62%|######2   | 6218/10000 [00:05<00:03, 1211.54it/s]
 63%|######3   | 6340/10000 [00:05<00:03, 1212.08it/s]
 65%|######4   | 6462/10000 [00:05<00:02, 1212.19it/s]
 66%|######5   | 6584/10000 [00:05<00:02, 1212.84it/s]
 67%|######7   | 6706/10000 [00:05<00:02, 1212.80it/s]
 68%|######8   | 6828/10000 [00:05<00:02, 1213.00it/s]
 70%|######9   | 6950/10000 [00:05<00:02, 1213.07it/s]
 71%|#######   | 7072/10000 [00:05<00:02, 1213.26it/s]
 72%|#######1  | 7194/10000 [00:05<00:02, 1212.99it/s]
 73%|#######3  | 7316/10000 [00:06<00:02, 1212.93it/s]
 74%|#######4  | 7438/10000 [00:06<00:02, 1213.56it/s]
 76%|#######5  | 7560/10000 [00:06<00:02, 1214.57it/s]
 77%|#######6  | 7682/10000 [00:06<00:01, 1214.10it/s]
 78%|#######8  | 7804/10000 [00:06<00:01, 1213.40it/s]
 79%|#######9  | 7926/10000 [00:06<00:01, 1212.32it/s]
 80%|########  | 8048/10000 [00:06<00:01, 1212.34it/s]
 82%|########1 | 8170/10000 [00:06<00:01, 1212.49it/s]
 83%|########2 | 8292/10000 [00:06<00:01, 1212.80it/s]
 84%|########4 | 8414/10000 [00:06<00:01, 1210.90it/s]
 85%|########5 | 8536/10000 [00:07<00:01, 1211.56it/s]
 87%|########6 | 8658/10000 [00:07<00:01, 1211.48it/s]
 88%|########7 | 8780/10000 [00:07<00:01, 1211.52it/s]
 89%|########9 | 8902/10000 [00:07<00:00, 1211.73it/s]
 90%|######### | 9024/10000 [00:07<00:00, 1212.14it/s]
 91%|#########1| 9146/10000 [00:07<00:00, 1212.21it/s]
 93%|#########2| 9268/10000 [00:07<00:00, 1212.65it/s]
 94%|#########3| 9390/10000 [00:07<00:00, 1212.68it/s]
 95%|#########5| 9512/10000 [00:07<00:00, 1212.73it/s]
 96%|#########6| 9634/10000 [00:07<00:00, 1212.39it/s]
 98%|#########7| 9756/10000 [00:08<00:00, 1212.92it/s]
 99%|#########8| 9878/10000 [00:08<00:00, 1213.09it/s]
100%|##########| 10000/10000 [00:08<00:00, 1213.17it/s]
100%|##########| 10000/10000 [00:08<00:00, 1212.71it/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"]

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.Dataset>
Dimensions:  (chain: 32, draw: 10000)
Coordinates:
  * chain    (chain) int64 0 1 2 3 4 5 6 7 8 9 ... 22 23 24 25 26 27 28 29 30 31
  * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 9993 9994 9995 9996 9997 9998 9999
Data variables:
    m0       (chain, draw) float64 -0.0001787 -0.000345 ... -6.586 -6.586
    m1       (chain, draw) float64 -2.378e-05 0.0002541 ... -5.178 -5.178
    m2       (chain, draw) float64 -9.894e-05 -0.0002379 ... 2.169 2.169
    m3       (chain, draw) float64 -0.0001222 -0.0004012 ... 1.081 1.081
Attributes:
    created_at:                 2022-12-14T02:06:34.443873
    arviz_version:              0.12.1
    inference_library:          emcee
    inference_library_version:  3.1.2


# a standard `trace` plot
axes = az.plot_trace(az_idata, backend_kwargs={"constrained_layout":True});

# add legends
for i, axes_pair in enumerate(axes):
    ax1 = axes_pair[0]
    ax2 = axes_pair[1]
    ax1.axvline(true_model[i], linestyle='dotted', color='red')
    ax1.set_xlabel("parameter value")
    ax1.set_ylabel("density value")
    ax2.set_xlabel("number of iterations")
    ax2.set_ylabel("parameter value")
m0, m0, m1, m1, m2, m2, m3, m3
tau = sampler.get_autocorr_time()
print(f"autocorrelation time: {tau}")
autocorrelation time: [68.29219144 76.96574306 68.65854986 73.0627644 ]
# a Corner plot

fig, axes = plt.subplots(nparams, nparams, figsize=(12,8))

if(False): # if we are plotting the model ensemble use this
    az.plot_pair(
        az_idata.sel(draw=slice(300,None)),
        marginals=True,
        reference_values=dict(zip([f"m{i}" for i in range(4)], true_model.tolist())),
        ax=axes,
    );
else: # if we wish to plot a kernel density plot then use this option
    az.plot_pair(
        az_idata.sel(draw=slice(300,None)),
        marginals=True,
        reference_values=dict(zip([f"m{i}" for i in range(4)], true_model.tolist())),
        kind="kde",
        kde_kwargs={
            "hdi_probs": [0.3, 0.6, 0.9],  # Plot 30%, 60% and 90% HDI contours
            "contourf_kwargs": {"cmap": "Blues"},
        },
        ax=axes,
    );
1 linear regression

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()
plot_models(flat_samples[inds])
plot_model(x,true_y, "True model", color="darkorange")
1 linear regression

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  -5.718 [ -6.430,  -5.019]
m1  -5.114 [ -5.598,  -4.632]
m2   1.826 [  1.457,   2.204]
m3   0.976 [  0.843,   1.109]
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
 [[ 0.18833668  0.05827435 -0.08045926 -0.02521073]
 [ 0.05827435  0.08626787 -0.03301353 -0.01800818]
 [-0.08045926 -0.03301353  0.05142114  0.0168967 ]
 [-0.02521073 -0.01800818  0.0168967   0.00669061]]

 Posterior estimate of model standard deviations in each parameter
    m0  0.4340
    m1  0.2937
    m2  0.2268
    m3  0.0818

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 = (-1.,-10.,-10.,-10.)\), and \({\mathbf u}^T = (2.,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?

Upload to Jamboard 2

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,-10,-10,-10])             # lower bound for uniform prior
m_upper_bound = np.array([2,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])
plot_model(x, true_y, "True model", color="darkorange")
1 linear regression
  0%|          | 0/10000 [00:00<?, ?it/s]
  1%|1         | 118/10000 [00:00<00:08, 1174.97it/s]
  2%|2         | 239/10000 [00:00<00:08, 1190.69it/s]
  4%|3         | 360/10000 [00:00<00:08, 1195.72it/s]
  5%|4         | 480/10000 [00:00<00:07, 1195.07it/s]
  6%|6         | 600/10000 [00:00<00:07, 1193.87it/s]
  7%|7         | 721/10000 [00:00<00:07, 1196.80it/s]
  8%|8         | 841/10000 [00:00<00:07, 1194.99it/s]
 10%|9         | 962/10000 [00:00<00:07, 1197.48it/s]
 11%|#         | 1083/10000 [00:00<00:07, 1199.79it/s]
 12%|#2        | 1203/10000 [00:01<00:07, 1198.03it/s]
 13%|#3        | 1323/10000 [00:01<00:07, 1190.41it/s]
 14%|#4        | 1443/10000 [00:01<00:07, 1186.44it/s]
 16%|#5        | 1563/10000 [00:01<00:07, 1189.83it/s]
 17%|#6        | 1684/10000 [00:01<00:06, 1193.31it/s]
 18%|#8        | 1804/10000 [00:01<00:06, 1191.66it/s]
 19%|#9        | 1925/10000 [00:01<00:06, 1194.84it/s]
 20%|##        | 2046/10000 [00:01<00:06, 1196.87it/s]
 22%|##1       | 2167/10000 [00:01<00:06, 1198.06it/s]
 23%|##2       | 2288/10000 [00:01<00:06, 1198.92it/s]
 24%|##4       | 2408/10000 [00:02<00:06, 1195.40it/s]
 25%|##5       | 2529/10000 [00:02<00:06, 1197.39it/s]
 26%|##6       | 2650/10000 [00:02<00:06, 1199.45it/s]
 28%|##7       | 2771/10000 [00:02<00:06, 1200.41it/s]
 29%|##8       | 2892/10000 [00:02<00:05, 1201.15it/s]
 30%|###       | 3013/10000 [00:02<00:05, 1195.88it/s]
 31%|###1      | 3133/10000 [00:02<00:05, 1196.78it/s]
 33%|###2      | 3254/10000 [00:02<00:05, 1197.91it/s]
 34%|###3      | 3374/10000 [00:02<00:05, 1198.52it/s]
 35%|###4      | 3494/10000 [00:02<00:05, 1198.58it/s]
 36%|###6      | 3614/10000 [00:03<00:05, 1194.91it/s]
 37%|###7      | 3735/10000 [00:03<00:05, 1196.91it/s]
 39%|###8      | 3856/10000 [00:03<00:05, 1198.89it/s]
 40%|###9      | 3977/10000 [00:03<00:05, 1200.58it/s]
 41%|####      | 4098/10000 [00:03<00:04, 1201.41it/s]
 42%|####2     | 4219/10000 [00:03<00:04, 1197.17it/s]
 43%|####3     | 4340/10000 [00:03<00:04, 1198.61it/s]
 45%|####4     | 4461/10000 [00:03<00:04, 1199.33it/s]
 46%|####5     | 4582/10000 [00:03<00:04, 1199.89it/s]
 47%|####7     | 4702/10000 [00:03<00:04, 1198.99it/s]
 48%|####8     | 4822/10000 [00:04<00:04, 1195.48it/s]
 49%|####9     | 4943/10000 [00:04<00:04, 1197.94it/s]
 51%|#####     | 5064/10000 [00:04<00:04, 1199.54it/s]
 52%|#####1    | 5185/10000 [00:04<00:04, 1200.38it/s]
 53%|#####3    | 5306/10000 [00:04<00:03, 1195.00it/s]
 54%|#####4    | 5426/10000 [00:04<00:03, 1192.27it/s]
 55%|#####5    | 5547/10000 [00:04<00:03, 1195.26it/s]
 57%|#####6    | 5668/10000 [00:04<00:03, 1197.34it/s]
 58%|#####7    | 5789/10000 [00:04<00:03, 1198.61it/s]
 59%|#####9    | 5909/10000 [00:04<00:03, 1190.83it/s]
 60%|######    | 6029/10000 [00:05<00:03, 1188.86it/s]
 62%|######1   | 6150/10000 [00:05<00:03, 1192.40it/s]
 63%|######2   | 6271/10000 [00:05<00:03, 1195.64it/s]
 64%|######3   | 6392/10000 [00:05<00:03, 1197.94it/s]
 65%|######5   | 6513/10000 [00:05<00:02, 1199.52it/s]
 66%|######6   | 6633/10000 [00:05<00:02, 1181.72it/s]
 68%|######7   | 6753/10000 [00:05<00:02, 1185.26it/s]
 69%|######8   | 6873/10000 [00:05<00:02, 1188.21it/s]
 70%|######9   | 6994/10000 [00:05<00:02, 1192.32it/s]
 71%|#######1  | 7115/10000 [00:05<00:02, 1195.31it/s]
 72%|#######2  | 7235/10000 [00:06<00:02, 1193.06it/s]
 74%|#######3  | 7356/10000 [00:06<00:02, 1195.75it/s]
 75%|#######4  | 7477/10000 [00:06<00:02, 1197.58it/s]
 76%|#######5  | 7598/10000 [00:06<00:02, 1199.52it/s]
 77%|#######7  | 7719/10000 [00:06<00:01, 1200.36it/s]
 78%|#######8  | 7840/10000 [00:06<00:01, 1185.76it/s]
 80%|#######9  | 7960/10000 [00:06<00:01, 1189.38it/s]
 81%|########  | 8080/10000 [00:06<00:01, 1192.05it/s]
 82%|########2 | 8201/10000 [00:06<00:01, 1194.59it/s]
 83%|########3 | 8322/10000 [00:06<00:01, 1196.45it/s]
 84%|########4 | 8442/10000 [00:07<00:01, 1193.44it/s]
 86%|########5 | 8562/10000 [00:07<00:01, 1195.35it/s]
 87%|########6 | 8683/10000 [00:07<00:01, 1196.96it/s]
 88%|########8 | 8804/10000 [00:07<00:00, 1198.05it/s]
 89%|########9 | 8925/10000 [00:07<00:00, 1198.85it/s]
 90%|######### | 9045/10000 [00:07<00:00, 1185.96it/s]
 92%|#########1| 9165/10000 [00:07<00:00, 1188.66it/s]
 93%|#########2| 9286/10000 [00:07<00:00, 1192.33it/s]
 94%|#########4| 9406/10000 [00:07<00:00, 1194.43it/s]
 95%|#########5| 9527/10000 [00:07<00:00, 1196.13it/s]
 96%|#########6| 9647/10000 [00:08<00:00, 1192.55it/s]
 98%|#########7| 9767/10000 [00:08<00:00, 1194.69it/s]
 99%|#########8| 9887/10000 [00:08<00:00, 1196.13it/s]
100%|##########| 10000/10000 [00:08<00:00, 1195.23it/s]
Resulting samples with prior model lower bounds of [-1,-10,-10,-10], upper bounds of [2,10,10,10]

Why do you think the posterior distribution looks like this?


Challenge: Change the data uncertainty#

To change the data uncertainty we increase sigma and then redefine the log-Likelihood.

Here we increase the assumed data standard deviation by a factor of of 50! So we are telling the inversion that the data are far less accurate than they actually are.

sigma = 50.0                                     # common noise standard deviation
Cdinv = np.eye(len(data_y))/(sigma**2)      # inverse data covariance matrix

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

Lets return the prior to the original bounds.

m_lower_bound = np.ones(4) * (-10.)             # lower bound for uniform prior
m_upper_bound = np.ones(4) * 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)

Your challenge is then to tell CoFI that the Likelihood and prior have changed and then to rerun the sample, and plot results.

Upload to Jamboard 3

Feel free to start from the code below:

######## CoFI BaseProblem - update information
inv_problem.set_log_likelihood(<CHANGE ME>)
inv_problem.set_log_prior(<CHANGE ME>)

######## CoFI Inversion - run it
inv_5 = Inversion(inv_problem, inv_options_3)
inv_result_5 = inv_5.run()

flat_samples = inv_result_5.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 from changed data uncertainty")
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

The answer is in the next cells if you want to run them.

#@title Solution

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

######## CoFI Inversion - run it
inv_5 = Inversion(inv_problem, inv_options_3)
inv_result_5 = inv_5.run()

flat_samples = inv_result_5.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 from changed data uncertainty")
plot_data()
plot_models(flat_samples[inds])
plot_model(x,true_y, "True model", color="darkorange")
1 linear regression
  0%|          | 0/10000 [00:00<?, ?it/s]
  1%|1         | 120/10000 [00:00<00:08, 1194.22it/s]
  2%|2         | 243/10000 [00:00<00:08, 1211.54it/s]
  4%|3         | 366/10000 [00:00<00:07, 1217.31it/s]
  5%|4         | 489/10000 [00:00<00:07, 1220.76it/s]
  6%|6         | 612/10000 [00:00<00:07, 1216.79it/s]
  7%|7         | 735/10000 [00:00<00:07, 1219.95it/s]
  9%|8         | 858/10000 [00:00<00:07, 1221.46it/s]
 10%|9         | 981/10000 [00:00<00:07, 1223.14it/s]
 11%|#1        | 1104/10000 [00:00<00:07, 1223.84it/s]
 12%|#2        | 1227/10000 [00:01<00:07, 1219.29it/s]
 13%|#3        | 1349/10000 [00:01<00:07, 1215.03it/s]
 15%|#4        | 1471/10000 [00:01<00:07, 1214.87it/s]
 16%|#5        | 1594/10000 [00:01<00:06, 1216.90it/s]
 17%|#7        | 1717/10000 [00:01<00:06, 1218.98it/s]
 18%|#8        | 1840/10000 [00:01<00:06, 1219.82it/s]
 20%|#9        | 1963/10000 [00:01<00:06, 1220.81it/s]
 21%|##        | 2086/10000 [00:01<00:06, 1222.30it/s]
 22%|##2       | 2209/10000 [00:01<00:06, 1223.27it/s]
 23%|##3       | 2332/10000 [00:01<00:06, 1224.01it/s]
 25%|##4       | 2455/10000 [00:02<00:06, 1222.83it/s]
 26%|##5       | 2578/10000 [00:02<00:06, 1223.06it/s]
 27%|##7       | 2701/10000 [00:02<00:05, 1223.75it/s]
 28%|##8       | 2824/10000 [00:02<00:05, 1224.62it/s]
 29%|##9       | 2947/10000 [00:02<00:05, 1224.63it/s]
 31%|###       | 3070/10000 [00:02<00:05, 1223.22it/s]
 32%|###1      | 3193/10000 [00:02<00:05, 1223.36it/s]
 33%|###3      | 3316/10000 [00:02<00:05, 1223.18it/s]
 34%|###4      | 3439/10000 [00:02<00:05, 1223.45it/s]
 36%|###5      | 3562/10000 [00:02<00:05, 1223.40it/s]
 37%|###6      | 3685/10000 [00:03<00:05, 1222.03it/s]
 38%|###8      | 3808/10000 [00:03<00:05, 1217.97it/s]
 39%|###9      | 3930/10000 [00:03<00:04, 1217.78it/s]
 41%|####      | 4053/10000 [00:03<00:04, 1218.45it/s]
 42%|####1     | 4175/10000 [00:03<00:04, 1218.65it/s]
 43%|####2     | 4297/10000 [00:03<00:04, 1214.65it/s]
 44%|####4     | 4419/10000 [00:03<00:04, 1215.83it/s]
 45%|####5     | 4542/10000 [00:03<00:04, 1217.34it/s]
 47%|####6     | 4664/10000 [00:03<00:04, 1217.95it/s]
 48%|####7     | 4786/10000 [00:03<00:04, 1218.52it/s]
 49%|####9     | 4908/10000 [00:04<00:04, 1214.57it/s]
 50%|#####     | 5031/10000 [00:04<00:04, 1216.25it/s]
 52%|#####1    | 5154/10000 [00:04<00:03, 1218.19it/s]
 53%|#####2    | 5277/10000 [00:04<00:03, 1218.99it/s]
 54%|#####4    | 5400/10000 [00:04<00:03, 1219.95it/s]
 55%|#####5    | 5522/10000 [00:04<00:03, 1219.22it/s]
 56%|#####6    | 5644/10000 [00:04<00:03, 1218.94it/s]
 58%|#####7    | 5766/10000 [00:04<00:03, 1218.77it/s]
 59%|#####8    | 5889/10000 [00:04<00:03, 1219.29it/s]
 60%|######    | 6011/10000 [00:04<00:03, 1218.06it/s]
 61%|######1   | 6133/10000 [00:05<00:03, 1217.30it/s]
 63%|######2   | 6255/10000 [00:05<00:03, 1209.08it/s]
 64%|######3   | 6376/10000 [00:05<00:03, 1189.13it/s]
 65%|######4   | 6499/10000 [00:05<00:02, 1198.42it/s]
 66%|######6   | 6622/10000 [00:05<00:02, 1205.25it/s]
 67%|######7   | 6744/10000 [00:05<00:02, 1208.88it/s]
 69%|######8   | 6867/10000 [00:05<00:02, 1212.39it/s]
 70%|######9   | 6990/10000 [00:05<00:02, 1215.25it/s]
 71%|#######1  | 7113/10000 [00:05<00:02, 1216.78it/s]
 72%|#######2  | 7236/10000 [00:05<00:02, 1217.88it/s]
 74%|#######3  | 7358/10000 [00:06<00:02, 1217.35it/s]
 75%|#######4  | 7481/10000 [00:06<00:02, 1218.23it/s]
 76%|#######6  | 7603/10000 [00:06<00:01, 1218.55it/s]
 77%|#######7  | 7725/10000 [00:06<00:01, 1218.64it/s]
 78%|#######8  | 7847/10000 [00:06<00:01, 1218.58it/s]
 80%|#######9  | 7969/10000 [00:06<00:01, 1218.17it/s]
 81%|########  | 8091/10000 [00:06<00:01, 1218.44it/s]
 82%|########2 | 8214/10000 [00:06<00:01, 1219.25it/s]
 83%|########3 | 8336/10000 [00:06<00:01, 1219.19it/s]
 85%|########4 | 8458/10000 [00:06<00:01, 1219.31it/s]
 86%|########5 | 8580/10000 [00:07<00:01, 1210.59it/s]
 87%|########7 | 8702/10000 [00:07<00:01, 1210.36it/s]
 88%|########8 | 8824/10000 [00:07<00:00, 1213.02it/s]
 89%|########9 | 8947/10000 [00:07<00:00, 1215.49it/s]
 91%|######### | 9069/10000 [00:07<00:00, 1216.53it/s]
 92%|#########1| 9191/10000 [00:07<00:00, 1216.47it/s]
 93%|#########3| 9314/10000 [00:07<00:00, 1217.79it/s]
 94%|#########4| 9437/10000 [00:07<00:00, 1219.34it/s]
 96%|#########5| 9560/10000 [00:07<00:00, 1220.49it/s]
 97%|#########6| 9683/10000 [00:07<00:00, 1221.41it/s]
 98%|#########8| 9806/10000 [00:08<00:00, 1219.72it/s]
 99%|#########9| 9929/10000 [00:08<00:00, 1219.98it/s]
100%|##########| 10000/10000 [00:08<00:00, 1217.79it/s]
Resulting samples from changed data uncertainty

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

Now lets decrease the number of steps performed by the McMC algorithm. It will be faster but perform less exploration of the model parameters.

We suggest you reduce the number of steps taken by all 32 random walkers and see how it affects the posterior ensemble.

Upload to Jamboard 4

You can start from code template below:

# change number of steps
nsteps = <CHANGE ME>              # instead of 10000

# change number of walkers
nwalkers = <CHANGE ME>            # instead of 32
walkers_start = np.zeros(nparams) + 1e-4 * np.random.randn(nwalkers, ndim)

# let's return to the old uncertainty settings
sigma = 1.0                                     # common noise standard deviation
Cdinv = np.eye(len(data_y))/(sigma**2)      # inverse data covariance matrix

######## CoFI InversionOptions - get a different tool
inv_options_3.set_params(nsteps=nsteps, nwalkers=nwalkers, initial_state=walkers_start)

######## CoFI Inversion - run it
inv_6 = Inversion(inv_problem, inv_options_3)
inv_result_6 = inv_6.run()

######## CoFI InversionResult - plot result
flat_samples = inv_result_6.sampler.get_chain(discard=300, thin=30, flat=True)
inds = np.random.randint(len(flat_samples), size=10) # get a random selection from posterior ensemble

print(f"Inference results from {nsteps} steps and {nwalkers} walkers")
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

# change number of steps
nsteps = 400              # instead of 10000

# change number of walkers
nwalkers = 30             # instead of 32
walkers_start = np.zeros(nparams) + 1e-4 * np.random.randn(nwalkers, ndim)

# let's return to the old uncertainty settings
sigma = 1.0                                     # common noise standard deviation
Cdinv = np.eye(len(data_y))/(sigma**2)      # inverse data covariance matrix

######## CoFI InversionOptions - get a different tool
inv_options_3.set_params(nsteps=nsteps, nwalkers=nwalkers, initial_state=walkers_start)

######## CoFI Inversion - run it
inv_6 = Inversion(inv_problem, inv_options_3)
inv_result_6 = inv_6.run()

######## CoFI InversionResult - plot result
flat_samples = inv_result_6.sampler.get_chain(discard=300, thin=30, flat=True)
inds = np.random.randint(len(flat_samples), size=10) # get a random selection from posterior ensemble

print(f"Inference results from {nsteps} steps and {nwalkers} walkers")
plot_data()
plot_models(flat_samples[inds])
plot_model(x,true_y, "True model", color="darkorange")
1 linear regression
  0%|          | 0/400 [00:00<?, ?it/s]
 31%|###       | 123/400 [00:00<00:00, 1222.74it/s]
 62%|######2   | 250/400 [00:00<00:00, 1245.87it/s]
 94%|#########3| 376/400 [00:00<00:00, 1252.22it/s]
100%|##########| 400/400 [00:00<00:00, 1248.22it/s]
Inference results from 400 steps and 30 walkers

Where to next?#


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.1.2.dev22
numpy 1.21.6
scipy 1.9.1
matplotlib 3.5.3
emcee 3.1.2
arviz 0.12.1

sphinx_gallery_thumbnail_number = -1

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

Gallery generated by Sphinx-Gallery