Linear regression with 1D-correlation#

The n data points (y1, y2, …, yn) generated for this example are sampled from an n-variate normal distribution with mean values given by yi = a * xi + b with a, b being the model parameters and x1, x2, …, xi, …, xn being predefined spatial x-coordinates ranging from 0 to 1. The data points (y1, y2, …, yn) are not independent but correlated in x. The corresponding covariance matrix is defined based on an exponential correlation function parameterized by the constant standard deviation sigma of the n-variate normal distribution and a correlation length l_corr. Hence, the full model has four parameters a, b, sigma, l_corr, all of which are inferred in this example using maximum likelihood estimation as well as sampling via emcee.

First, let’s import the required functions and classes for this example.

# third party imports
import numpy as np
import matplotlib.pyplot as plt
from tripy.utils import correlation_function
from tripy.utils import correlation_matrix

# local imports (problem definition)
from probeye.definition.inverse_problem import InverseProblem
from probeye.definition.forward_model import ForwardModelBase
from probeye.definition.distribution import Normal, Uniform
from probeye.definition.sensor import Sensor
from probeye.definition.likelihood_model import GaussianLikelihoodModel
from probeye.definition.correlation_model import ExpModel

# local imports (problem solving)
from probeye.inference.scipy.solver import MaxLikelihoodSolver
from probeye.inference.emcee.solver import EmceeSolver

# local imports (inference data post-processing)
from probeye.postprocessing.sampling_plots import create_pair_plot
from probeye.postprocessing.sampling_plots import create_posterior_plot
from probeye.postprocessing.sampling_plots import create_trace_plot

We start by generating a synthetic data set from a known linear model to which we will add correlated noise. Afterwards, we will pretend to have forgotten the parameters of this ground-truth model and will instead try to recover them just from the data. The slope (a) and intercept (b) of the ground truth model are set to be:

# ground truth
a_true = 2.5
b_true = 1.7

Now, let’s generate a few data points that we contaminate with a Gaussian error:

# settings for the data generation
n_experiments = 3
n_points = 50
seed = 1
std_noise = 0.5
l_corr = 0.05

# first create the true values without an error model; these 'true' values will be the
# mean values for sampling from a multivariate normal distribution that accounts for the
# intended correlation
np.random.seed(seed)
x_test = np.linspace(0.0, 1.0, n_points)
y_true = a_true * x_test + b_true

# assemble the spatial covariance matrix
x_test_as_column_matrix = x_test.reshape((n_points, -1))
f_corr = lambda a: correlation_function(d=a, correlation_length=l_corr)
cov = std_noise**2 * correlation_matrix(x_test_as_column_matrix, f_corr)

# now generate the noisy test data including correlations; we assume here that
# there are n_experiments test series
data_dict = {}
for i in range(n_experiments):
    y_test_i = np.random.multivariate_normal(mean=y_true, cov=cov)
    data_dict[f"Test_{i}"] = y_test_i
    plt.scatter(x_test, y_test_i, label=f"measured data (test {i+1})", s=10, zorder=10)
# finish the plot
plt.plot(x_test, y_true, label="true model", c="black")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
plot correlation 1D

Until this point, we didn’t use probeye at all, since we just generate some data. In a normal use case, we wouldn’t have to generate our data of course. Instead, it would be provided to us, for example as the result of some test series. As the first step in any calibration problem, one needs to have a parameterized model (in probeye such a model is called ‘forward model’) of which one assumes that it is able to describe the data at hand. In this case, if we took a look at the blue data points in the plot above without knowing the orange line, we might expect a simple linear model. It is now our job to describe this model within the probeye-framework. This is done by defining our own specific model class:

class LinearModel(ForwardModelBase):
    def interface(self):
        self.parameters = ["a", "b"]
        self.input_sensors = Sensor("x")
        self.output_sensors = Sensor("y", std_model="sigma")

    def response(self, inp: dict) -> dict:
        x = inp["x"]
        m = inp["a"]
        b = inp["b"]
        return {"y": m * x + b}

First, note that this model class is based on the probeye class ‘ForwardModelBase’. While this is a requirement, the name of the class can be chosen freely. As you can see, this class has a ‘interface’ and a ‘response’ method. In the ‘interface’ method we define that our model has two parameters, ‘a’ and ‘b’, next to one input and one output sensors, called ‘x’ and ‘y’ respectively. Keeping this interface in mind, let’s now take a look at the ‘response’ method. This method describes the actual forward model evaluation. The method takes one dictionary as an input and returns one dictionary as its output. The input dictionary ‘inp’ will have the keys ‘a’, ‘b’ and ‘x’ because of the definitions given in self.interface. Analogously, the returned dictionary must have the key ‘y’, because we defined an output sensor with the name ‘y’. Note that the entire interface of the ‘response’ method is described by the ‘interface’ method. Parameters and input sensors will be contained in the ‘inp’ dictionary, while the output sensors must be contained in the returned dictionary.

After we now have defined our forward model, we can set up the inverse problem itself. This always begins by initializing an object form the InverseProblem-class, and adding all of the problem’s parameters with priors that reflect our current best guesses of what the parameter’s values might look like. Please check out the ‘Components’-part of this documentation to get more information on the arguments seen below. However, most of the code should be self-explanatory.

# initialize the inverse problem with a useful name
problem = InverseProblem("Linear regression with 1D correlation", print_header=False)

# add the problem's parameters
problem.add_parameter(
    "a",
    tex="$a$",
    info="Slope of the graph",
    prior=Normal(mean=2.0, std=1.0),
)
problem.add_parameter(
    "b",
    info="Intersection of graph with y-axis",
    tex="$b$",
    prior=Normal(mean=1.0, std=1.0),
)
problem.add_parameter(
    "sigma",
    domain="(0, +oo)",
    tex=r"$\sigma$",
    info="Standard deviation, of zero-mean Gaussian noise model",
    prior=Uniform(low=0.0, high=0.8),
)
problem.add_parameter(
    "l_corr",
    domain="(0, +oo)",
    tex=r"$l_\mathrm{corr}$",
    info="Correlation length of correlation model",
    prior=Uniform(low=0.0, high=0.2),
)

As the next step, we need to add our experimental data the forward model and the likelihood model. Note that the order is important and should not be changed.

# experimental data
for exp_name, y_test_i in data_dict.items():
    problem.add_experiment(
        name=exp_name,
        sensor_data={"x": x_test, "y": y_test_i},
    )

# add the forward model to the problem
problem.add_forward_model(LinearModel("LinearModel"), experiments=[*data_dict.keys()])

# likelihood model
for exp_name in data_dict:
    likelihood_model = GaussianLikelihoodModel(
        experiment_name=exp_name,
        model_error="additive",
        correlation=ExpModel(x="l_corr"),
    )
    problem.add_likelihood_model(likelihood_model)

Now, our problem definition is complete, and we can take a look at its summary:

# print problem summary
problem.info(print_header=True)
2023-08-02 14:21:47.461 | INFO     | # ================================================================================================ # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.461 | INFO     | #                                                                                                  # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                                            dP                                                    # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                                            88                                                    # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  88d888b. 88d888b..d8888b. 88d888b. .d8888b. dP    dP .d8888b.                   # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  88'  `88 88'     88'  `88 88'  `88 88ooood8 88    88 88ooood8                   # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  88.  .88 88      88.  .88 88.  .88 88.      88.  .88 88.                        # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  88Y888P' dP      `88888P' 88Y8888' `88888P' `8888P88 `88888P'                   # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  88                                                .88                           # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                  dP                                            d8888P                            # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | #                                                                                                  # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.462 | INFO     | # ================================================================================================ # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.463 | INFO     | #                                                                                                  # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.463 | INFO     | #        Version 3.0.4 - A general framework for setting up parameter estimation problems.         # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.463 | INFO     | #                                                                                                  # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.463 | INFO     | # ================================================================================================ # | probeye.subroutines:print_probeye_header:620
2023-08-02 14:21:47.465 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.465 | INFO     | Problem summary: Linear regression with 1D correlation                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.465 | INFO     | ======================================================                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     | Forward models                                                                                       | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     | ---------------------------------------------------------                                            | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     |  Model name   | Global parameters   | Local parameters                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     | --------------+---------------------+--------------------                                            | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     |  LinearModel  | a, b                | a, b                                                           | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     | Priors                                                                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.466 | INFO     | ------------------------------------------------------------------------------------                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  Prior name     | Global parameters               | Local parameters                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     | ----------------+---------------------------------+---------------------------------                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  a_normal       | a, mean_a, std_a                | a, mean_a, std_a                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  b_normal       | b, mean_b, std_b                | b, mean_b, std_b                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  sigma_uniform  | sigma, low_sigma, high_sigma    | sigma, low_sigma, high_sigma                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  l_corr_uniform | l_corr, low_l_corr, high_l_corr | l_corr, low_l_corr, high_l_corr                  | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     | Parameter overview                                                                                   | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     | ---------------------------------------------------------------------------------------------------------------- | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.467 | INFO     |  Parameter type/role   | Parameter names                                                              |   Count | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     | -----------------------+------------------------------------------------------------------------------+--------- | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Model parameters      | a, b                                                                         |       2 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Prior parameters      | mean_a, std_a, mean_b, std_b, low_sigma, high_sigma, low_l_corr, high_l_corr |       8 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Likelihood parameters | sigma, l_corr                                                                |       2 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Const parameters      | mean_a, std_a, mean_b, std_b, low_sigma, high_sigma, low_l_corr, high_l_corr |       8 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Latent parameters     | a, b, sigma, l_corr                                                          |       4 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     | Parameter explanations                                                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     | -----------------------------------------------------------------------                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.468 | INFO     |  Name        | Short explanation                                                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     | -------------+---------------------------------------------------------                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  mean_a      | Normal prior's parameter for latent parameter 'a'                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  std_a       | Normal prior's parameter for latent parameter 'a'                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  a           | Slope of the graph                                                                    | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  mean_b      | Normal prior's parameter for latent parameter 'b'                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  std_b       | Normal prior's parameter for latent parameter 'b'                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  b           | Intersection of graph with y-axis                                                     | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  low_sigma   | Uniform prior's parameter for latent parameter 'sigma'                                | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  high_sigma  | Uniform prior's parameter for latent parameter 'sigma'                                | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.469 | INFO     |  sigma       | Standard deviation, of zero-mean Gaussian noise model                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  low_l_corr  | Uniform prior's parameter for latent parameter 'l_corr'                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  high_l_corr | Uniform prior's parameter for latent parameter 'l_corr'                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  l_corr      | Correlation length of correlation model                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     | Constant parameters                                                                                  | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     | -----------------------                                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  Name        |   Value                                                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     | -------------+---------                                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  mean_a      |     2                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.470 | INFO     |  std_a       |     1                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  mean_b      |     1                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  std_b       |     1                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  low_sigma   |     0                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  high_sigma  |     0.8                                                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  low_l_corr  |     0                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |  high_l_corr |     0.2                                                                               | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.471 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | Theta interpretation                                                                                 | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | +---------------------------+                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |  Theta  |    Parameter    |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |  index  |      name       |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |---------------------------|                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |      0 --> a              |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |      1 --> b              |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |      2 --> sigma          |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | |      3 --> l_corr         |                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.472 | INFO     | +---------------------------+                                                                        | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     | Added experiments                                                                                    | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     | --------------------------------------------                                                         | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |  Name   | Sensor values   | Forward model                                                            | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     | --------+-----------------+-----------------                                                         | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |  Test_0 | x (50 elements) | LinearModel                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |         | y (50 elements) |                                                                          | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |  Test_1 | x (50 elements) | LinearModel                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |         | y (50 elements) |                                                                          | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.473 | INFO     |  Test_2 | x (50 elements) | LinearModel                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |         | y (50 elements) |                                                                          | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     | Added likelihood models                                                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     | ----------------------------------------------------------                                           | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |  Name   | Parameters    | Target sensors   | Experiment                                              | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     | --------+---------------+------------------+--------------                                           | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |  Test_0 | l_corr, sigma | y                | Test_0                                                  | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |  Test_1 | l_corr, sigma | y                | Test_1                                                  | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |  Test_2 | l_corr, sigma | y                | Test_2                                                  | probeye.definition.inverse_problem:info:978
2023-08-02 14:21:47.474 | INFO     |                                                                                                      | probeye.definition.inverse_problem:info:978

After the problem definition comes the problem solution. There are different solver one can use, but we will just demonstrate how to use two of them: the scipy-solver, which merely provides a point estimate based on a maximum likelihood optimization, and the emcee solver, which is a MCMC-sampling solver. Let’s begin with the scipy-solver:

# this is for using the scipy-solver (maximum likelihood estimation)
scipy_solver = MaxLikelihoodSolver(problem, show_progress=False)
max_like_data = scipy_solver.run()
2023-08-02 14:21:47.478 | INFO     | Solving problem via maximum likelihood estimation                                                    | probeye.inference.scipy.solver:_run_ml_or_map:452
2023-08-02 14:21:47.479 | INFO     | Using start values:                                                                                  | probeye.inference.scipy.solver:_run_ml_or_map:460
2023-08-02 14:21:47.479 | INFO     | a       = 2.0                                                                                        | probeye.subroutines:print_dict_in_rows:728
2023-08-02 14:21:47.479 | INFO     | b       = 1.0                                                                                        | probeye.subroutines:print_dict_in_rows:728
2023-08-02 14:21:47.480 | INFO     | sigma   = 0.4                                                                                        | probeye.subroutines:print_dict_in_rows:728
2023-08-02 14:21:47.480 | INFO     | l_corr  = 0.1                                                                                        | probeye.subroutines:print_dict_in_rows:728
2023-08-02 14:21:47.480 | INFO     | Starting optimizer (using Nelder-Mead)                                                               | probeye.inference.scipy.solver:_run_ml_or_map:464
2023-08-02 14:21:47.480 | INFO     | No solver options specified                                                                          | probeye.inference.scipy.solver:_run_ml_or_map:469
2023-08-02 14:21:47.568 | INFO     |                                                                                                      | probeye.inference.scipy.solver:summarize_point_estimate_results:356
2023-08-02 14:21:47.568 | INFO     | Results of maximum likelihood estimation                                                             | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.568 | INFO     | =====================================                                                                | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.568 | INFO     | Optimization terminated successfully.                                                                | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.568 | INFO     | -------------------------------------                                                                | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.568 | INFO     | Number of iterations:           131                                                                  | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.568 | INFO     | Number of function evaluations: 225                                                                  | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.569 | INFO     | -------------------------------------                                                                | probeye.inference.scipy.solver:summarize_point_estimate_results:359
2023-08-02 14:21:47.569 | INFO     | a_opt      = [2.50521064] (start = 2.0)                                                              | probeye.inference.scipy.solver:summarize_point_estimate_results:382
2023-08-02 14:21:47.569 | INFO     | b_opt      = [1.57413752] (start = 1.0)                                                              | probeye.inference.scipy.solver:summarize_point_estimate_results:382
2023-08-02 14:21:47.569 | INFO     | sigma_opt  = [0.45072398] (start = 0.4)                                                              | probeye.inference.scipy.solver:summarize_point_estimate_results:382
2023-08-02 14:21:47.570 | INFO     | l_corr_opt = [0.05469275] (start = 0.1)                                                              | probeye.inference.scipy.solver:summarize_point_estimate_results:382
2023-08-02 14:21:47.570 | INFO     |                                                                                                      | probeye.inference.scipy.solver:summarize_point_estimate_results:383

All solver have in common that they are first initialized, and then execute a run-method, which returns its result data in the format of an arviz inference-data object (except for the scipy-solver). Let’s now take a look at the emcee-solver.

emcee_solver = EmceeSolver(problem, show_progress=False)
inference_data = emcee_solver.run(n_steps=2000, n_initial_steps=200)
2023-08-02 14:21:47.574 | INFO     | Solving problem using emcee sampler with 200 + 2000 samples and 20 walkers                           | probeye.inference.emcee.solver:run:178
2023-08-02 14:21:47.574 | INFO     | No additional options specified                                                                      | probeye.inference.emcee.solver:run:186
2023-08-02 14:22:31.892 | INFO     | Sampling of the posterior distribution completed: 2000 steps and 20 walkers.                         | probeye.inference.emcee.solver:run:255
2023-08-02 14:22:31.892 | INFO     | Total run-time (including initial sampling): 44s.                                                    | probeye.inference.emcee.solver:run:259
2023-08-02 14:22:31.892 | INFO     |                                                                                                      | probeye.inference.emcee.solver:run:260
2023-08-02 14:22:31.893 | INFO     | Summary of sampling results (emcee)                                                                  | probeye.inference.emcee.solver:run:261
2023-08-02 14:22:31.900 | INFO     |           mean    median    sd    5%    95%                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     | ------  ------  --------  ----  ----  -----                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     | a         2.48      2.48  0.30  2.01   2.99                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     | b         1.57      1.58  0.18  1.26   1.86                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     | sigma     0.51      0.50  0.07  0.42   0.65                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     | l_corr    0.07      0.07  0.03  0.04   0.12                                                          | probeye.inference.emcee.solver:emcee_summary:137
2023-08-02 14:22:31.901 | INFO     |                                                                                                      | probeye.inference.emcee.solver:run:267

Finally, we want to plot the results we obtained. To that end, probeye provides some post-processing routines, which are mostly based on the arviz-plotting routines.

# this is optional, since in most cases we don't know the ground truth
true_values = {"a": a_true, "b": b_true, "sigma": std_noise, "l_corr": l_corr}

# this is an overview plot that allows to visualize correlations
pair_plot_array = create_pair_plot(
    inference_data,
    emcee_solver.problem,
    true_values=true_values,
    focus_on_posterior=True,
    show_legends=False,
    title="Sampling results from emcee-Solver (pair plot)",
)
Sampling results from emcee-Solver (pair plot)
# this is a posterior plot, without including priors
post_plot_array = create_posterior_plot(
    inference_data,
    emcee_solver.problem,
    true_values=true_values,
    title="Sampling results from emcee-Solver (posterior plot)",
)
Sampling results from emcee-Solver (posterior plot), $a$, $b$, $\sigma$, $l_\mathrm{corr}$
# trace plots are used to check for "healthy" sampling
trace_plot_array = create_trace_plot(
    inference_data,
    emcee_solver.problem,
    title="Sampling results from emcee-Solver (trace plot)",
)
Sampling results from emcee-Solver (trace plot), $a$, $a$, $b$, $b$, $\sigma$, $\sigma$, $l_\mathrm{corr}$, $l_\mathrm{corr}$

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