BISIP

DOCUMENTATION IS UNDER CONSTRUCTION

BISIP is an MIT licensed Python package to perform Bayesian Inversion of Spectral Induced Polarization data with Markov-chain Monte Carlo simulation. See our original paper in Computers & Geosciences for more details. It should be noted that the original version of BISIP was built on top of the PyMC framework when the code was developed in 2015. Now, BISIP uses Goodman & Weare’s Affine Invariant Ensemble sampler as implemented in the emcee library to explore the SIP models’ parameter spaces with multiple walkers. You can read the original paper by Foreman-Mackey et al. explaining the emcee algorithm in detail here.

BISIP is being developed on GitHub.

Basic Usage

To perform a Debye decomposition inversion of a SIP data file using all default parameters, you would use the following code:

from bisip import PolynomialDecomposition

# Use one of the example data files provided with BISIP
filepath = '/path/to/bisip/data/SIP-K389175.dat'

model = PolynomialDecomposition(filepath)

# Fit the model to this data file
model.fit()

#   Out:  100%|██████████| 1000/1000 [00:01<00:00, 563.92it/s]

# Plot the fit on top of the input data, discarding the burn-in period
model.plot_fit(discard=200)

A more detailed example is available in the Quickstart tutorial.

Installation

BISIP is compatible with Python 3.6+.

Requirements

The following packages are required and should be installed automatically on setup:

• numpy
• cython
• matplotlib
• emcee

These optional packages are used for progress bars and corner plots:

• tqdm
• corner

Package managers

TODO: Add BISIP to conda-forge.

From source

BISIP is developed on GitHub. Clone the repository to your computer. Then navigate to the bisip directory. Finally run the setup.py script with Python.

git clone https://github.com/clberube/bisip2
cd bisip2
python setup.py install -f

Testing

To test if everything was installed correctly, do the following:

# Last tested on Python 3.7.3 (default, Mar 27 2019, 16:54:48)
import bisip
bisip.run_test()

If everything is OK the code will load a data file and perform inversion of a data file using various models. Then it will print the best parameters and plot traces and fit quality for the last model used. At the end, you should see the following line:

All tests passed. Press ctrl+C or close figure windows to exit.

The Inversion Class

Standard usage of bisip involves instantiating an Inversion object. This is normally done by invoking one of the inversion model classes described in the SIP models section below.

class bisip.models.Inversion(filepath, nwalkers=32, nsteps=5000, headers=1, ph_units='mrad')

Bases: bisip.plotlib.plotlib, bisip.utils.utils

An abstract class to perform inversion of SIP data.

This class is the base constructor for the PolynomialDecomposition and ColeCole classes.

Parameters:

• filepath (str) – The path to the file to perform inversion on.
• nwalkers (int) – Number of walkers to use to explore the parameter space. Defaults to 32.
• nsteps (int) – Number of steps to perform in the MCMC simulation. Defaults to 5000.
• headers (int) – The number of header lines in the file. Defaults to 1.
• ph_units (str) – The units of the phase shift measurements. Choices: ‘mrad’, ‘rad’, ‘deg’. Defaults to ‘mrad’.

data

The input data dictionary.

Type:dict

fit(p0=None, pool=None, moves=None)

Samples the posterior distribution to fit the model to the data.

Parameters:

• p0 (ndarray) – Starting parameter values. Should be a 2D array with shape (nwalkers, ndim). If None, random values will be uniformly drawn from the parameter bounds. Defaults to None.
• pool (pool, optional) – A pool object from the Python multiprocessing library. See https://emcee.readthedocs.io/en/stable/tutorials/parallel/. Defaults to None.
• moves (moves, optional) – A emcee Moves class (see https://emcee.readthedocs.io/en/stable/user/moves/). If None, the emcee algorithm StretchMove is used. Defaults to None.

fitted

Whether the model has been fitted or not.

Type:bool

get_chain(**kwargs)

Gets the MCMC chains from a fitted model.

Keyword Arguments:  

• discard (int) – Number of steps to discard (burn-in period).
• thin (int) – Thinning factor.
• flat (bool) – Whether or not to flatten the walkers. If flat is False, the output chain will have shape (nsteps, nwalkers, ndim). If flat is True, the output chain will have shape (nsteps*nwalkers, ndim).

Returns:

The MCMC chain(s).

Return type:

ndarray

p0

Starting parameter values. Should be a 2D array with shape (nwalkers, ndim).

Type:ndarray

param_bounds

Ordered bounds of the parameters.

Type:list of float

param_names

Ordered names of the parameters.

Type:list of str

params

Parameter names and their bounds.

Type:dict

sampler

A emcee sampler object (see https://emcee.readthedocs.io/en/stable/user/sampler/).

Type:EnsembleSampler

SIP models

Pelton Cole-Cole

class bisip.models.PeltonColeCole(*args, n_modes=1, **kwargs)

Bases: bisip.models.Inversion

A generalized ColeCole inversion scheme for SIP data.

Parameters:

• *args – Arguments passed to the Inversion class.
• n_modes (int) – The number of ColeCole modes to use for the inversion. Defaults to 1.
• **kwargs – Additional keyword arguments passed to the Inversion class.

forward(theta, w)

Returns a ColeCole impedance.

Parameters:

• theta (ndarray) – Ordered array of R0, m_{1}, …, m_{n_modes}, log_tau_{1}, …, log_tau_{n_modes}, c_{1}, …, c_{n_modes}. See https://doi.org/10.1016/j.cageo.2017.05.001.
• w (ndarray) – Array of angular frequencies to compute the impedance for (w = 2*pi*f).

Polynomial decomposition

class bisip.models.PolynomialDecomposition(*args, poly_deg=5, c_exp=1.0, **kwargs)

Bases: bisip.models.Inversion

A polynomial decomposition inversion scheme for SIP data.

Parameters:

• *args – Arguments passed to the Inversion class.
• poly_deg (int) – The polynomial degree to use for the decomposition. Defaults to 5.
• c_exp (float) – The c-exponent to use for the decomposition scheme. 0.5 -> Warburg, 1.0 -> Debye. Defaults to 1.0.
• **kwargs – Additional keyword arguments passed to the Inversion class.

forward(theta, w)

Returns a Polynomial Decomposition impedance.

Parameters:

• theta (ndarray) – Ordered array of R0, a_{poly_deg}, a_{poly_deg-1}, …, a_{0}. See https://doi.org/10.1016/j.cageo.2017.05.001.
• w (ndarray) – Array of angular frequencies to compute the impedance for (w = 2*pi*f).

Dias (2000)

class bisip.models.Dias2000(*args, **kwargs)

Bases: bisip.models.Inversion

A generalized ColeCole inversion scheme for SIP data.

Parameters:

• *args – Arguments passed to the Inversion class.
• **kwargs – Additional keyword arguments passed to the Inversion class.

forward(theta, w)

Returns a Dias (2000) impedance.

Parameters:

• theta (ndarray) – Ordered array of R0, m, log_tau, eta, delta. See https://doi.org/10.1016/j.cageo.2017.05.001.
• w (ndarray) – Array of angular frequencies to compute the impedance for (w = 2*pi*f).

Shin (2015)

class bisip.models.Shin2015(*args, **kwargs)

Bases: bisip.models.Inversion

A Shin (2015) inversion scheme for SIP data.

Parameters:

• *args – Arguments passed to the Inversion class.
• **kwargs – Additional keyword arguments passed to the Inversion class.

Warning

The Shin model implementation is yielding unexpected results and needs to be reviewed.

forward(theta, w)

Returns a Shin (2015) impedance.

Parameters:

• theta (ndarray) – Ordered array of R1, R2, log_Q1, log_Q2, n1, n2.
• w (ndarray) – Array of angular frequencies to compute the impedance for (w = 2*pi*f).

Plotting methods

These functions may be called as methods of the Inversion class after fitting the model to a dataset.

class bisip.plotlib.plotlib

plot_corner(chain=None, **kwargs)

Plots the corner plot of the MCMC simulation.

Parameters:

• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should have a shape (nsteps, nwalkers, ndim) or (nsteps*nwalkers, ndim). If None, the full, unflattened chain will be used and all walkers will be plotted. Defaults to None.
• **kwargs – Additional keyword arguments for the get_chain function (see below). Use these arguments only if not explicitly passing the chain argument.

Keyword Arguments:  

• discard (int) – The number of steps to discard.
• thin (int) – The thinning factor (keep every thin step).
• flat (bool) – Whether to flatten the walkers into a single chain or not.

Returns:

A matplotlib figure.

Return type:

Figure

plot_data(feature='phase', **kwargs)

Plots the input data.

Parameters:

• feature (str) – Which data feature to plot. Choices are ‘phase’, ‘amplitude’, ‘real’, ‘imaginary’. Defaults to ‘phase’.
• **kwargs – Additional keyword arguments passed to the matplotlib errorbar function.

Returns:

A matplotlib figure.

Return type:

Figure

plot_fit(chain=None, p=[2.5, 50, 97.5], **kwargs)

Plots the input data, best fit and confidence interval of a model. Shows the real and imaginary parts.

Parameters:

• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should have a shape (nsteps, nwalkers, ndim) or (nsteps*nwalkers, ndim). If None, the full, unflattened chain will be used and all walkers will be plotted. Defaults to None.
• p (list of int) – Percentile values for lower confidence interval, best fit curve, and upper confidence interval, in that order. Defaults to [2.5, 50, 97.5] for the median and 95% HPD.
• **kwargs – Additional keyword arguments for the get_chain function (see below). Use these arguments only if not explicitly passing the chain argument.

Keyword Arguments:  

• discard (int) – The number of steps to discard.
• thin (int) – The thinning factor (keep every thin step).
• flat (bool) – Whether to flatten the walkers into a single chain or not.

Returns:

A matplotlib figure.

Return type:

Figure

plot_fit_pa(chain=None, p=[2.5, 50, 97.5], **kwargs)

Plots the input data, best fit and confidence interval of a model. Shows the amplitude and phase spectra.

Parameters:

• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should have a shape (nsteps, nwalkers, ndim) or (nsteps*nwalkers, ndim). If None, the full, unflattened chain will be used and all walkers will be plotted. Defaults to None.
• p (list of int) – Percentile values for lower confidence interval, best fit curve, and upper confidence interval, in that order. Defaults to [2.5, 50, 97.5] for the median and 95% HPD.
• **kwargs – Additional keyword arguments for the get_chain function (see below). Use these arguments only if not explicitly passing the chain argument.

Keyword Arguments:  

• discard (int) – The number of steps to discard.
• thin (int) – The thinning factor (keep every thin step).
• flat (bool) – Whether to flatten the walkers into a single chain or not.

Returns:

A matplotlib figure.

Return type:

Figure

plot_histograms(chain=None, bins=25, **kwargs)

Plots histograms of the MCMC simulation chains.

Parameters:

• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should have a shape (nsteps, nwalkers, ndim) or (nsteps*nwalkers, ndim). If None, the full, unflattened chain will be used and all walkers will be plotted. Defaults to None.
• bins (int) – The number of bins to use in the histograms.
• **kwargs – Additional keyword arguments for the get_chain function (see below). Use these arguments only if not explicitly passing the chain argument.

Keyword Arguments:  

• discard (int) – The number of steps to discard.
• thin (int) – The thinning factor (keep every thin step).
• flat (bool) – Whether to flatten the walkers into a single chain or not.

Returns:

A matplotlib figure.

Return type:

Figure

plot_traces(chain=None, **kwargs)

Plots the traces of the MCMC simulation.

Parameters:

• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should have a shape (nsteps, nwalkers, ndim) or (nsteps*nwalkers, ndim). If None, the full, unflattened chain will be used and all walkers will be plotted. Defaults to None.
• **kwargs – Additional keyword arguments for the get_chain function (see below). Use these arguments only if not explicitly passing the chain argument.

Keyword Arguments:  

• discard (int) – The number of steps to discard.
• thin (int) – The thinning factor (keep every thin step).
• flat (bool) – Whether to flatten the walkers into a single chain or not.

Returns:

A matplotlib figure.

Return type:

Figure

Utility methods

These utility functions may be called as methods of the Inversion class.

class bisip.utils.utils

get_model_percentile(p=[2.5, 50, 97.5], chain=None, **kwargs)

Gets percentiles of the model values for a MCMC chain.

Parameters:

• p (float or list of float) – percentiles values to compute. Defaults to [2.5, 50, 97.5].
• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should be a 2D array (nsteps, ndim). If None and no kwargs are passed to discard iterations, will raise a warning and the full chain will be used. Defaults to None.

Keyword Arguments:  

**kwargs – See kwargs of the get_chain method.

get_param_mean(chain=None, **kwargs)

Gets the mean of the model parameters for a MCMC chain.

Parameters:chain (ndarray) – A numpy array containing the MCMC chain to plot. Should be a 2D array (nsteps, ndim). If None and no kwargs are passed to discard iterations, will raise a warning and the full chain will be used. Defaults to None. Keyword Arguments:  **kwargs – See kwargs of the get_chain method.

get_param_percentile(p=[2.5, 50, 97.5], chain=None, **kwargs)

Gets percentiles of the parameter values for a MCMC chain.

Parameters:

• p (float or list of float) – percentiles values to compute. Defaults to [2.5, 50, 97.5].
• chain (ndarray) – A numpy array containing the MCMC chain to plot. Should be a 2D array (nsteps, ndim). If None and no kwargs are passed to discard iterations, will raise a warning and the full chain will be used. Defaults to None.

Keyword Arguments:  

**kwargs – See kwargs of the get_chain method.

get_param_std(chain=None, **kwargs)

Gets the standard deviation of the model parameters.

Parameters:chain (ndarray) – A numpy array containing the MCMC chain to plot. Should be a 2D array (nsteps, ndim). If None and no kwargs are passed to discard iterations, will raise a warning and the full chain will be used. Defaults to None. Keyword Arguments:  **kwargs – See kwargs of the get_chain method.

load_data(filename, headers=1, ph_units='mrad')

Imports a data file and prepares it for inversion.

Parameters:

• filepath (str) – The path to the data file.
• headers (int) – The number of header lines in the file. Defaults to 1.
• ph_units (str) – The units of the phase shift measurements. Choices: ‘mrad’, ‘rad’, ‘deg’. Defaults to ‘mrad’.

print_latex_parameters(names, values, uncertainties, decimals=3)

Prints Pelton parameters and their uncertainties with LaTeX.

Requires the code to be run in interactive mode either from a Jupyter notebook or IPython console to use display and Math functions.

Parameters:

• names (list) – The list of parameter names. Obtain it from Inversion.param_names.
• values (list) – Ordered list of mean parameter values. Obtain it from Inversion.get_param_mean().
• uncertainties (list) – Ordered list of parameter uncertainties. Obtain it from Inversion.get_param_std().
• decimals (int) – The number of decimals to display.

Data format

It is important to respect this format when preparing your data files. We use the format output of the SIP-Fuchs-III instrument and software.

Frequency, Amplitude,   Phase shift,  Amplit error, Phase error
6.000e+03, 1.17152e+05, -2.36226e+02, 1.171527e+01, 9.948376e-02
3.000e+03, 1.22177e+05, -1.46221e+02, 1.392825e+01, 1.134464e-01
1.500e+03, 1.25553e+05, -9.51099e+01, 2.762214e+01, 2.199114e-01
........., ..........., ............, ............, ............
........., ..........., ............, ............, ............
........., ..........., ............, ............, ............
4.575e-02, 1.66153e+05, -1.21143e+01, 1.947314e+02, 1.171115e+00
2.288e-02, 1.67988e+05, -9.36718e+00, 3.306003e+02, 1.9669860+00
1.144e-02, 1.70107e+05, -7.25533e+00, 5.630541e+02, 3.310889e+00

Note

• Save your data in .csv, .txt, .dat or any other format. The extension is not important as long as it is a ASCII file.
• Comma separation between columns is mandatory.
• The order of the columns is very important (Frequency, Amplitude, Phase shift, Amplitude error, Phase error).
• Phase units may be milliradians, radians or degrees.
• Units are specified as an initialization argument (e.g. ph_units=’mrad’).
• Amplitude units may be Ohm-m or Ohm, the data will be normalized.
• A number of header lines may be skipped function argument (e.g. in this case headers=1).
• Scientific or standard notation is OK.

Example data files

Several real-life data files are included with the BISIP package. These files contain the complex resistivity spectra of mineralized rock samples. Once you have installed BISIP in your Python environment, you may get the paths to these data files using the DataFiles class.

class bisip.data.DataFiles(*args, **kwargs)

Bases: dict

The following example shows how to get the absolute file paths to the BISIP data files in your Python installation.

from bisip import DataFiles

files = DataFiles()  # initialize the dictionary
print(files.keys())  # see the various data file names

filepath = files['SIP-K389172']  # extract one path
print(filepath)

Quickstart

In this first tutorial we load a data file, perform single-mode Pelton Cole-Cole inversion on it, visualize the fit quality and posterior distribution and save results to a csv file. A more in-depth look at the Pelton Cole-Cole model and its parameters is available in the Pelton tutorial.

Running your first inversion

To invert a SIP data file, you would use the following approach:

• Import the base PeltonColeCole model.
• Pass filepath to instantiate the model with a specific data file.
• Set the simulation to explore the Pelton model parameter space with nwalkers=32.
• Set the simulation to run for 2000 steps by passing nsteps=2000.
• Skip the first 9 lines in the data file with headers=9.
• Fit the model to the data by calling the fit() method.

import os

import numpy as np
from bisip import PeltonColeCole
from bisip import DataFiles

# This will get one of the example data files in the BISIP package
data_files = DataFiles()
filepath = data_files['SIP-K389172']

model = PeltonColeCole(filepath=filepath,
                       nwalkers=32,  # number of MCMC walkers
                       nsteps=2000,  # number of MCMC steps
                       headers=9,  # number of lines to skip in the data file
                       )

# Fit the model to this data file
model.fit()

100%|██████████| 2000/2000 [00:03<00:00, 593.58it/s]

Note

In this example we use headers=9 to skip the header line and the 8 highest measurement frequencies. As a result we keep only a single polarization mode in the lower frequency range.

Visualizing the parameter traces

Visual inspection of the traces of each MCMC walker is a good way to quickly assess how well the inversion went. They also define what we consider samples of the posterior distribution. The walkers are initialized at random positions in the parameter space, let’s inspect the parameter traces to see if they have moved towards a unique solution to the inverse problem. Plotting the traces is one of the first things the user should do after completing the fit() method. Simply call the plot_traces() method of the Inversion object to visualize the parameter chains.

# Plot the parameter traces
fig = model.plot_traces()

The walkers have correctly moved to similar solutions. As you can see for the first few iterations the walkers are spread out across the entire parameter space. The ensemble has then reached a stationary state after at least 50 steps. To be safe, let’s discard all the values before the 500th step and use the following steps to estimate the best values for our parameters.

You can also store the chain as a numpy ndarray using the get_chain() method. Passing the argument flat=True will flatten the walkers into a single chain, whereas passing flat=False will preserve the nwalkers dimension. The thin argument will pick one sample every thin steps, this can be used to reduce autocorrelation.

# Get chains of all walkers, discarding first 500 steps
chain = model.get_chain(discard=500)
print(chain.shape)  # (nsteps, nwalkers, ndim)

(1500, 32, 4)

# Get chains of all walkers, discarding first 500 steps,
# thinning by a factor of 2 and flattening the walkers
chain = model.get_chain(discard=500, thin=2, flat=True)
print(chain.shape)  # (nsteps*nwalkers/thin, ndim)

(24000, 4)

Note

Here ndim refers to the number of dimensions of the model (the number of parameters).

Plotting models over data

Let’s now plot the best model obtained with this chain using the plot_fit() method, which will show the real and imaginary parts of the complex resistivity as a function of measurement frequency. The plot will also display the median model as a red line and the 95% highest probability density interval as dotted lines.

# Use the chain argument to plot the model
# for a specific flattened chain
fig = model.plot_fit(chain)

# You can then use the `fig` matplotlib object to save the figure or make
# adjustments according to your personal preference. For example:
# fig.savefig('fit_figure_sample_K389172.png', dpi=144, bbox_inches='tight')

You may also directly pass the arguments of the get_chain() method to most plotting functions. By default this will flatten the walkers into a single chain. For example here we use the plot_fit_pa() method, which will display the amplitude and phase shift as a function of measurement frequency. We also adjust the confidence interval to display the 72% HPD with the p argument (p for percentile).

# Use the discard argument to discard samples
# note that we do not need to pass the chain argument
fig = model.plot_fit_pa(discard=500, p=[14, 50, 86])  # 14=lower, 50=median, 86=higher

Printing best parameters

Here we will use the get_param_mean() and get_param_std() utility functions to extract the mean and standard deviations of each parameter in our chain. We will also access the name of each parameter with the param_names property.

values = model.get_param_mean(chain=chain)
uncertainties = model.get_param_std(chain=chain)

for n, v, u in zip(model.param_names, values, uncertainties):
  print(f'{n}: {v:.5f} +/- {u:.5f}')

r0: 1.02413 +/- 0.00484
m1: 0.36238 +/- 0.02619
log_tau1: -2.13531 +/- 0.25427
c1: 0.50193 +/- 0.03324

Looks good! The relative error on the recovered parameters is quite small considering our large error bars on the data.

Note

It is important to note that for every inversion scheme the amplitude values (and the r0 parameter) have been normalized by the maximum amplitude. You may access this normalization factor with model.data['norm_factor']. Therefore the real \(\rho_0\) value is r0 * model.data['norm_factor'].

Inspecting the posterior

Let’s now visualize the posterior distribution of all parameters using a corner plot (from the corner Python package).

fig = model.plot_corner(chain=chain)

This is interesting, we note that the \(\log \tau\) values are anticorrelated with the \(m\) values. It seems like the model can obtain similar fits of the data by adjusting these parameters against each other.

Saving results to CSV file

Finally let’s save the best parameters and some of their statistics as a csv file.

# Get the lower, median and higher percentiles
results = model.get_param_percentile(chain=chain, p=[2.5, 50, 97.5])
# Join the list of parameter names into a comma separated string
headers = ','.join(model.param_names)
# Save to csv with numpy
# The first row is the 2.5th percentile, 2nd the 50th (median), 3rd the 97.5th.
# Parameter names will be listed in the csv file header.

print(headers)

r0,m1,log_tau1,c1

print(results)

[[ 1.01547398  0.31429044 -2.66165185  0.43494749]
 [ 1.02389413  0.3611087  -2.12615936  0.50230639]
 [ 1.03427943  0.41689845 -1.66573259  0.56613694]]

np.savetxt('quickstart_results.csv', results, header=headers,
           delimiter=',', comments='')

The Pelton Cole-Cole model

Introduction

Pelton (1978) noted that the complex resistivity spectra of rocks \(\rho^*\) could be interpreted using Cole-Cole relaxation models:

\[\rho^* = \rho_0 \left[ 1-m\left(1-\frac{1}{1+(i\omega\tau)^c} \right) \right],\]

where \(\omega\) is the measurement angular frequency (\(\omega=2\pi f\)) and \(i\) is the imaginary unit.

BISIP implements the generalized form of the Pelton Cole-Cole model given by Chen et al., (2008):

\[\rho^* = \rho_0 \left[ 1 - \sum_{k=1}^{K} m_k\left(1-\frac{1}{1+(i\omega\tau_k)^{c_k}} \right) \right],\]

where the subscript \(k\) refers to one of \(K\) superimposed Cole-Cole relaxation modes. Here, \(\rho^*\) depends on \(1 + 3K\) parameters:

• \(\rho_0 \in [0, \infty)\), the direct current resistivity \(\rho_0 = \rho^*(\omega\to 0)\).
• \(m_k \in [0, 1]\), each mode’s chargeability \(m=(\rho_0 - \rho_\infty)/\rho_0\).
• \(\tau_k \in [0, \infty)\), each mode’s relaxation time, related to average polarizable particle size.
• \(c_k \in [0, 1]\), the slope of the measured phase shift for each mode.

In this tutorial we will show that the parameter space of the generalized Cole-Cole model can be quite complex for some data sets.

Defining the number of modes

First import the required packages.

import os

import numpy as np
from bisip import PeltonColeCole
from bisip import DataFiles

np.random.seed(42)

To invert a SIP data file with the Pelton model, start by instantiating a PeltonColeCole object and call the fit() method to estimate the parameters.

# This will get one of the example data files in the BISIP package
data_files = DataFiles()
filepath = data_files['SIP-K389174']

# Define MCMC parameters and ColeCole model
model = PeltonColeCole(filepath=filepath,
                       nwalkers=32,
                       nsteps=1000)
# Fit model to data
model.fit()

100%|██████████| 1000/1000 [00:01<00:00, 731.41it/s]

Let’s look at the fit quality first, discarding the first 500 steps. You can inspect the parameter traces to convince yourself that this is an appropriate burn-in period.

# Plot the fit by discarding the first 500 steps
fig = model.plot_fit(discard=500)

Not good. clearly, our data is not well represented by the default ColeCole n_modes=1. In fact, we see two bumps in the imaginary part of the resistivity spectra (the first at 1 Hz and the second at > 10 kHz. Let’s redefine our model with n_modes = 2 and nwalkers = 64. We run the simulation with more walkers because the number of parameters to solve has increased.

# Add the n_modes = 2 argument.
model = PeltonColeCole(filepath=filepath,
                       nwalkers=64,
                       nsteps=1000,
                       n_modes=2)
# Fit model to data
model.fit()

100%|██████████| 1000/1000 [00:02<00:00, 435.49it/s]

Let’s check the traces again. We expect two values for \(m\), \(\log \tau\) and \(c\) because we have two Cole-Cole modes.

fig = model.plot_traces()

This is quite interesting. It’s not even worth yet to look at the fit quality, because the parameter traces tell us that there are at least two possible solutions to the inverse problem. This is because both Cole-Cole modes are interchangeable – meaning that mode 1 can take on the parameters of mode 2 and vice versa. We will now see how to fix this problem by adjusting the parameter boundaries.

First we note that \(\log \tau_1\) has two solutions: one around 0 and another around -15. The same is true for \(\log \tau_2\). We will fix \(\log \tau_1\) so that its values is constrained between -5 and 0, and fix the \(\log \tau_2\) to be constrained between -10 and -15.

# Adjust the boundaries
model.params.update(log_tau1=[-5, 5], log_tau2=[-15, -10])

Next we fit the model using these new boundaries and visualize the traces.

model.fit()
fig = model.plot_traces()

100%|██████████| 1000/1000 [00:02<00:00, 356.54it/s]

Amazing! The stricter priors have allowed all 64 walkers to find a unique solution. With these improved parameter chains the fit quality should be improved. Let’s plot the fit quality while discarding the first 1000 steps.

# Extract the chain
chain = model.get_chain(discard=500, thin=10, flat=True)

# Plot the fit by discarding the first 1000 steps
fig = model.plot_fit(chain)

The adjustment is satisfying, and the 95% HPD reasonable if we consider the relatively small measurement error bars. Let’s inspect the mean parameters and their uncertainties.

# Use this utility function to nicely print parameter values
# for the Pelton model with LaTeX

# Print the mean and std of the parameters after discarding burn-in samples
values = model.get_param_mean(chain)
uncertainties = model.get_param_std(chain)

model.print_latex_parameters(model.param_names, values, uncertainties)

$\displaystyle \rho_0: 1.010 \pm 0.002$

$\displaystyle m_1: 0.140 \pm 0.004$

$\displaystyle m_2: 0.935 \pm 0.055$

$\displaystyle \log\tau_1: -1.574 \pm 0.089$

$\displaystyle \log\tau_2: -12.838 \pm 0.148$

$\displaystyle c_1: 0.451 \pm 0.015$

$\displaystyle c_2: 0.608 \pm 0.015$

It looks like the parameters are very well resolved, and due to the small data error bars the uncertainties on the Pelton parameters are also very small (but you should still not neglect them). Feel free to visualize the posterior distribution of the Pelton model with the plot_corner method, which shows that the \(m_2\) parameter has accumulated at \(m_2 \to 1\). This tells us that although the fit quality is good, the Pelton model does not explain well our data set.

Large uncertainties

Let’s now consider a different rock sample where the measurement uncertainty is relatively lage (K389172), we apply the same method we did previously and put it all together in the following cell:

# Define MCMC parameters and ColeCole model with
# one of the example data files in the BISIP package
model = PeltonColeCole(filepath=data_files['SIP-K389172'],
                       nwalkers=64,
                       nsteps=1000,
                       n_modes=2)

# Adjust the boundaries
model.params.update(log_tau1=[-5, 5], log_tau2=[-20, -10])

# Fit model to data
model.fit()
fig = model.plot_traces()

100%|██████████| 1000/1000 [00:02<00:00, 366.72it/s]

The inversion has reached a unique solution, but due to the large error bars on the data we can see from the traces that the parameters of mode 2 (\(m_2\), \(\log \tau_2\), \(c_2\)) are quite ill-defined.

Weak polarization

Finally, let’s consider rock sample K389176, where the Cole-Cole modes are not well defined in the data (i.e. where the phase shift spectrum shows a relatively flat response.

# Define MCMC parameters and ColeCole model with
# one of the example data files in the BISIP package
model = PeltonColeCole(filepath=data_files['SIP-K389176'],
                       nwalkers=64,
                       nsteps=1000,
                       n_modes=2)

# Adjust the boundaries
model.params.update(log_tau1=[-5, 5], log_tau2=[-20, -10])

# Fit model to data
model.fit()

# Plot fit
fig = model.plot_fit(discard=500)

# Print the mean and std of the parameters after discarding burn-in samples
values = model.get_param_mean(discard=500)
uncertainties = model.get_param_std(discard=500)
model.print_latex_parameters(model.param_names, values, uncertainties, decimals=2)

100%|██████████| 1000/1000 [00:02<00:00, 393.54it/s]

$\displaystyle \rho_0: 1.02 \pm 0.01$

$\displaystyle m_1: 0.11 \pm 0.03$

$\displaystyle m_2: 0.73 \pm 0.16$

$\displaystyle \log\tau_1: -1.49 \pm 1.51$

$\displaystyle \log\tau_2: -13.02 \pm 0.91$

$\displaystyle c_1: 0.24 \pm 0.15$

$\displaystyle c_2: 0.73 \pm 0.13$

Inspect the traces and fit figures, you will notice that the data is well fitted by the double Pelton model. However, we note that due to the weak polarization (flat response) in the low-frequency range, the \(\log \tau_1\) parameter is extremely ill-defined (the relative error is more than 100%!). This is because the absence of a clear polarization peak (low values of \(c\)) can’t be explained with Cole-Cole type models, which expect a well-defined local maximum in the phase shift spectrum. This is discussed in more detail in Bérubé et al. (2017).

Conclusions

From this experiment, we learned that the Pelton Cole-Cole model modes are interchangeable, and that to avoid this problem we must fix the \(\log \tau_k\) parameters to take on boundaries that do not overlap. It also appears that the \(m_2\) parameter is often not well resolved because it has accumulated at its upper limit. This implies that even higher values of \(m_2\) may explain our data better. However, because \(m\) is defined in the interval \([0, 1]\) our results suggest that the Pelton Cole-Cole model may not explain well our data set, especially its high-frequency response. Finally, we saw that data error bars have a high impact of the quality of resolved Cole-Cole parameters, and that weak polarization responses will most certaintly yield ill-defined inversion parameters.

The Dias (2000) model

Introduction

The model proposed by Dias (1972) describes the petrophysical properties of rocks from measurements of their electrical polarization in a frequency range typically from 1 mHz to 100 kHz.

We refer to Dias (2000) for the implementation of the complex resistivity formula in BISIP. This model predicts that the complex resistivity spectra \(\rho^*\) of a polarizable rock sample can be described by

\begin{equation} \rho^* = \rho_0 \left[ 1-m\left(1-\frac{1}{1+i\omega\tau'\left(1+\frac{1}{\mu}\right)} \right) \right], \end{equation}

where \(\omega\) is the measurement angular frequencies (\(\omega=2\pi f\)) and \(i\) is the imaginary unit. Additionally,

\begin{equation} \mu = i\omega\tau + \left(i\omega\tau''\right)^{1/2}, \end{equation}

\begin{equation} \tau' = (\tau/\delta)(1 - \delta)/(1 - m), \end{equation}

\begin{equation} \tau'' = \tau^2 \eta^2. \end{equation}

Here, \(\rho^*\) depends on 5 parameters:

• \(\rho_0 \in [0, \infty)\), the direct current resistivity \(\rho_0 = \rho^*(\omega\to 0)\).
• \(m \in [0, 1)\), the chargeability \(m=(\rho_0 - \rho_\infty)/\rho_0\).
• \(\tau \in [0, \infty)\), the relaxation time, related to average polarizable particle size.
• \(\eta \in [0, 150]\), characteristic of the electrochemical environment producing polarization.
• \(\delta \in [0, 1)\), the pore length fraction of the electrical double layer zone in the material.

In this tutorial we will invert a SIP data with the Dias model to explore the parameter space of this semi-empirical SIP model. We will also quantify the relationship between the various parameters.

Exploring the parameter space

First import the required packages.

import os

import numpy as np
from IPython.display import display, Math
from bisip import Dias2000
from bisip import DataFiles

np.random.seed(42)

To invert a SIP data file with the Dias model, start by instantiating a Dias2000 object and call the fit() method to estimate the parameters.

# This will get one of the example data files in the BISIP package
data_files = DataFiles()
filepath = data_files['SIP-K389172']

# Define MCMC parameters and Dias model
nwalkers = 32
nsteps = 1000
model = Dias2000(filepath=filepath, nwalkers=nwalkers, nsteps=nsteps)
# Fit model to data
model.fit()

100%|██████████| 1000/1000 [00:01<00:00, 511.36it/s]

fig = model.plot_traces()

Some walkers get stuck in local minima because the priors are really wide. Nevertheless, we can see that the median solution of all these chains gives a satisfying result.

# Plot the fit by discarding the first 500 steps
fig = model.plot_fit(discard=500)

The adjustment is satisfying, but the 95% HPD is very wide because some of the walkers were stuck in local minima far from the solution. A good strategy to reduce the chance that walkers get stuck in local minima would be to tighten the priors are the values we think give a good result. Here we will set new boundaries for \(\eta\) and \(\log \tau\).

# Adjust the boundaries
model.params.update(eta=[0, 25], log_tau=[-15, -5])

model.fit()
fig = model.plot_traces()

100%|██████████| 1000/1000 [00:01<00:00, 533.24it/s]

The stricter priors have allowed all walkers to find a similar stationary state. With these improved parameter chains the fit quality should be improved.

# Plot the fit by discarding the first 500 steps
fig = model.plot_fit(discard=500)

The adjustment is satisfying, and the 95% HPD reasonable if we consider the measurement error bars. We will now visualize the posterior distribution of the Dias model with the plot_corner method.

# Plot the posterior by discarding the first 500 steps
fig = model.plot_corner(discard=500)

The corner plot shows interesting correlations between various parameters. Finally, let’s look at the optimal parameters and their uncertainties.

# Use this utility function to nicely print parameter values
# for the Dias model with LaTeX

# Print the mean and std of the parameters after discarding burn-in samples
values = model.get_param_mean(discard=500)
uncertainties = model.get_param_std(discard=500)

model.print_latex_parameters(model.param_names, values, uncertainties)

$\displaystyle \rho_0: 1.023 \pm 0.003$

$\displaystyle m: 0.688 \pm 0.058$

$\displaystyle \log\tau: -10.315 \pm 0.162$

$\displaystyle \eta: 7.767 \pm 0.838$

$\displaystyle \delta: 0.707 \pm 0.068$

Conclusions

From this experiment, we conclude that the \(\rho_0\) parameter is relatively independent from the others. We also note that \(m\) and \(\tau\) are characterized by a strong correlation coefficient. Most importantly, we find that this correlation makes the range of ‘best’ values for these parameters quite large, indicating that these parameters are not well resolved for this particular data file.

The polynomial decomposition

Introduction

The polynomial decomposition approach was originally proposed by Keery et al. (2012). It consists of performing Debye decomposition of the SIP spectra with the assumption that the relaxation time distribution (RTD) of that spectra can be approximated by a polynomial function. In this approach, the complex resistivity \(\rho^*\) is given by

\begin{equation} \rho^* = \rho_0 \left[ 1 - \sum_{l=1}^{L} m_l\left(1-\frac{1}{1+(i\omega\tau_l)} \right) \right], \end{equation}

where \(\omega\) is the measurement angular frequencies (\(\omega=2\pi f\)) and \(i\) is the imaginary unit. BISIP uses a generalized Cole-Cole polynomial decomposition scheme by including a user-defined Cole-Cole exponent (\(c\)). If \(c = 1\), the equation is equivalent to Debye decomposition. If \(c=0.5\), the equation is known as Warburg decomposition. For any other values we refer to this scheme as the generalized Cole-Cole decomposition. The complex resistivity is thus given by

\begin{equation} \rho^* = \rho_0 \left[ 1 - \sum_{l=1}^{L} m_l\left(1-\frac{1}{1+(i\omega\tau_l)^c} \right) \right], \end{equation}

where \(L\) is an arbitrary real number that defines the discretization level of the RTD. In practice BISIP uses \(L = 2\times n_{freq}\), where \(n_{freq}\) is the number of measurement points in the frequency range. The RTD is defined by the distribution of \(m_l\) values as a function of \(\tau_l\) values. According to Keery et al. (2012), the RTD can be written as

\begin{equation} m_l = \sum_{p=0}^P a_p (\log_{10} \tau_l)^p, \end{equation}

where \(a_p\) are polynomial coefficients and \(P\) is the polynomial degree of the approximation. In practice, \(P\) is often set to 3, 4 or 5. The user should aim to use the lowest degree possible while maintaining a satisfying fit with the data. As a result, \(\rho^*\) depends on \(P+1\) parameters:

• \(\rho_0 \in [0, \infty)\), the direct current resistivity \(\rho_0 = \rho^*(\omega\to 0)\).
• \(a_p\), the \(P\) polynomial coefficients.

From the recovered \(a_p\) and \(\rho_0\) values, the user can then compute the RTD and extract several integrating parameters, such as total chargeability and mean relaxation time. The user is referred to Keery et al. (2012) for the definition of these parameters. In this tutorial, we will compare the total chargeability values obtained with Debye and Warburg decomposition for 5 rock samples.

Debye decomposition

First import the required packages.

import numpy as np
from bisip import PolynomialDecomposition
from bisip import DataFiles

np.random.seed(42)

Let’s start by inspecting the data sets so we can decide which polynomial approximation would be appropriate to fit the observations.

data_files = DataFiles()  # get the dictionary of BISIP data files

for n, p in data_files.items():  # n is the filename, p the filepath
    model = PolynomialDecomposition(filepath=p)
    fig = model.plot_data('phase')
    ax = fig.gca()
    ax.set_title(n)

Samples K389173 and K389176 can be explained with a 3rd degree polynomial approxation, whereas the rest of the samples would likely require a 4th order polynomial approximation. For this example we will use a 4th order approaximation for all data sets.

# Make a dict to hold results
results_debye = {}
# Define MCMC and Debye parameters
params = {'nwalkers': 32,
          'nsteps': 1000,
          'c_exp': 1,
          'poly_deg': 4,
         }

for n, p in data_files.items():  # n is the filename, p the filepath
    # Initialize with 1 file at a time, and pass the params
    model = PolynomialDecomposition(filepath=p, **params)
    model.fit()
    # Store the fitted model in the results dict
    results_debye[n] = model

100%|██████████| 1000/1000 [00:01<00:00, 505.89it/s]
100%|██████████| 1000/1000 [00:01<00:00, 512.17it/s]
100%|██████████| 1000/1000 [00:01<00:00, 508.69it/s]
100%|██████████| 1000/1000 [00:01<00:00, 509.93it/s]
100%|██████████| 1000/1000 [00:01<00:00, 532.47it/s]
100%|██████████| 1000/1000 [00:01<00:00, 522.29it/s]

Now we have a dictionary of fitted models in results_debye. Let’s quickly inspect the fit results while discarding half the steps as burn-in. Check the chains with the plot_traces() method to see why the value is appropriate.

for name, model in results_debye.items():
    fig = model.plot_fit_pa(discard=500)
    ax = fig.gca()
    ax.set_title(n)

We have obtained a good fit and 95% HPD for all data sets. This is interesting because we have a unique solution for each data set. Therefore, our recovered values should be usable for all samples, whereas if we used a Cole-Cole model samples K389173 and K389176 would have a ill-defined solution (see the Pelton tutorial).

We need to do a bit of work to extract the total chargeability. We will use the pandas library to organize the results in a DataFrame.

import pandas as pd

names, models = zip(*results_debye.items())

mean_values = [m.get_param_mean(discard=500) for m in models]

df_debye = pd.DataFrame(mean_values)
df_debye.columns= models[0].param_names
df_debye.index = names
df_debye.head(6)

	r0	a0	a1	a2	a3	a4
SIP-K389170	0.989771	0.014033	-0.002350	-0.004634	-0.000338	0.000219
SIP-K389172	1.018388	0.017680	-0.011578	-0.003509	0.002061	0.000526
SIP-K389173	1.011923	0.003390	-0.000995	-0.000499	0.000396	0.000172
SIP-K389174	1.004851	0.007547	-0.003206	-0.001871	0.000548	0.000242
SIP-K389175	0.997613	0.006870	-0.003937	-0.001338	0.000741	0.000219
SIP-K389176	1.006941	0.002426	-0.001084	-0.000048	0.000564	0.000159

We need to sum the coefficients and log_tau values:

\begin{equation} m_l = \sum_{p=0}^P a_p (\log_{10} \tau_l)^p. \end{equation}

# Take log_tau from any model (considering the
# frequency range is the same for all samples)
log_tau = models[0].log_tau
print(log_tau)

[-6.         -5.79487179 -5.58974359 -5.38461538 -5.17948718 -4.97435897
 -4.76923077 -4.56410256 -4.35897436 -4.15384615 -3.94871795 -3.74358974
 -3.53846154 -3.33333333 -3.12820513 -2.92307692 -2.71794872 -2.51282051
 -2.30769231 -2.1025641  -1.8974359  -1.69230769 -1.48717949 -1.28205128
 -1.07692308 -0.87179487 -0.66666667 -0.46153846 -0.25641026 -0.05128205
  0.15384615  0.35897436  0.56410256  0.76923077  0.97435897  1.17948718
  1.38461538  1.58974359  1.79487179  2.        ]

Define this utility function to get the chargeability value (\(m\)) that corresponds to each \(\tau\) value.

def get_m(df, P):
    m = {}
    for idx in df.index:
        m[idx] = 0
        for p in range(P+1):
            coeff = df_debye.loc[idx, f'a{p}']
            m[idx] += coeff*log_tau**p
        # Add the total chargeability in the df
        df.loc[idx, 'total_m'] = np.sum(m[idx])
    return m

m_debye = get_m(df_debye, params['poly_deg'])
df_debye.head(6)

	r0	a0	a1	a2	a3	a4	total_m
SIP-K389170	0.989771	0.014033	-0.002350	-0.004634	-0.000338	0.000219	1.342871
SIP-K389172	1.018388	0.017680	-0.011578	-0.003509	0.002061	0.000526	1.210898
SIP-K389173	1.011923	0.003390	-0.000995	-0.000499	0.000396	0.000172	0.786556
SIP-K389174	1.004851	0.007547	-0.003206	-0.001871	0.000548	0.000242	0.930805
SIP-K389175	0.997613	0.006870	-0.003937	-0.001338	0.000741	0.000219	0.655028
SIP-K389176	1.006941	0.002426	-0.001084	-0.000048	0.000564	0.000159	0.545386

Let’s inspect the recovered RTDs:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
for n in names:
    ax.plot(10**log_tau, m_debye[n], label=n)
    ax.set_ylabel(r'$m$')
    ax.set_xlabel(r'$\tau$ (s)')
ax.set_title('Debye decomposition')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim([1e-5, 1e1])
ax.set_ylim([1e-3, 1e-1])
ax.legend(ncol=2)

<matplotlib.legend.Legend at 0x116aed160>

Warburg decomposition

Now we will repeat these steps for a Warburg decomposition (\(c = 0.5\)).

# Make a dict to hold results
results_warburg = {}
# Define MCMC and Debye parameters
params = {'nwalkers': 32,
          'nsteps': 1000,
          'c_exp': 0.5,
          'poly_deg': 4,
         }

for n, p in data_files.items():  # n is the filename, p the filepath
    # Initialize with 1 file at a time, and pass the params
    model = PolynomialDecomposition(filepath=p, **params)
    model.fit()
    # Store the fitted model in the results dict
    results_warburg[n] = model

100%|██████████| 1000/1000 [00:03<00:00, 262.07it/s]
100%|██████████| 1000/1000 [00:03<00:00, 265.81it/s]
100%|██████████| 1000/1000 [00:03<00:00, 263.03it/s]
100%|██████████| 1000/1000 [00:03<00:00, 264.32it/s]
100%|██████████| 1000/1000 [00:03<00:00, 269.38it/s]
100%|██████████| 1000/1000 [00:03<00:00, 268.75it/s]

names, models = zip(*results_debye.items())

mean_values = [m.get_param_mean(discard=500) for m in models]

df_warburg = pd.DataFrame(mean_values)
df_warburg.columns= models[0].param_names
df_warburg.index = names

m_warburg = get_m(df_warburg, params['poly_deg'])
df_warburg.head(6)

	r0	a0	a1	a2	a3	a4	total_m
SIP-K389170	0.989771	0.014033	-0.002350	-0.004634	-0.000338	0.000219	1.342871
SIP-K389172	1.018388	0.017680	-0.011578	-0.003509	0.002061	0.000526	1.210898
SIP-K389173	1.011923	0.003390	-0.000995	-0.000499	0.000396	0.000172	0.786556
SIP-K389174	1.004851	0.007547	-0.003206	-0.001871	0.000548	0.000242	0.930805
SIP-K389175	0.997613	0.006870	-0.003937	-0.001338	0.000741	0.000219	0.655028
SIP-K389176	1.006941	0.002426	-0.001084	-0.000048	0.000564	0.000159	0.545386

Add the Debye and Warburg values to a single DataFrame and print it.

df = pd.merge(df_debye, df_warburg,
              left_index=True, right_index=True,
              suffixes=('_Debye', '_Warburg'))
df[['total_m_Debye', 'total_m_Warburg']].head(6)

	total_m_Debye	total_m_Warburg
SIP-K389170	1.342871	1.342871
SIP-K389172	1.210898	1.210898
SIP-K389173	0.786556	0.786556
SIP-K389174	0.930805	0.930805
SIP-K389175	0.655028	0.655028
SIP-K389176	0.545386	0.545386

The values are identical. Finally let’s look at the RTDs to confirm that they are the same.

fig, ax = plt.subplots()
for n in names:
    ax.plot(10**log_tau, m_debye[n], label=n)
    ax.set_ylabel(r'$m$')
    ax.set_xlabel(r'$\tau$ (s)')
ax.set_title('Warburg decomposition')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim([1e-5, 1e1])
ax.set_ylim([1e-3, 1e-1])
ax.legend(ncol=2)

<matplotlib.legend.Legend at 0x116198828>

Conclusions

In this tutorial we have seen how to do batch inversions, organize the results in pandas DataFrames and how to extract integrating parameters from the polynomial decomposition approach. The total chargeability parameter was identical for Debye and Warburg decomposition approaches.

The Shin (2015) model

Introduction

The model proposed by Shin (2015) is an equivalent circuit that aims to reproduce SIP data This model predicts that the complex resistivity spectra \(\rho^*\) of a polarizable rock sample can be described by

\begin{equation} \rho^* = \sum_{i=1}^2 \frac{\rho_i}{(i\omega)^{n_i} \rho_iQ_i + 1} \end{equation}

where \(\omega\) is the measurement angular frequencies (\(\omega=2\pi f\)) and \(i\) is the imaginary unit.

Here, \(\rho^*\) depends on 3 pairs of parameters:

• \(\rho_i \in [0, \infty)\), the resistivity of the resistance element in Shin’s circuit.
• \(Q_i \in [0, \infty)\), the capacitance of the CPE.
• \(n_i \in [0, 1]\), the exponent of the CPE impedance (0 = resistor, 0.5 = warburg, 1.0 = capacitance).

In this tutorial we will perform batch inversion of all SIP data files provided with BISIP with the Shin (2015) and double Cole-Cole models, and we will compare their respective relaxation time (\(\tau\)) parameters.

Exploring the parameter space

First import the required packages.

import numpy as np
from bisip import PeltonColeCole
from bisip import Shin2015
from bisip import DataFiles

np.random.seed(42)

# Load the data file paths
data_files = DataFiles()

results = {'Shin': {},
           'Pelton': {},
          }

nsteps = 1000
for fname, fpath in data_files.items():

        if fname == 'SIP-K389175':
            model = Shin2015(fpath, nsteps=nsteps)
            model.fit()
            results['Shin'][fname] = model

#             model = PeltonColeCole(fpath, nsteps=nsteps, n_modes=2)
#             model.fit()
#             results['Pelton'][fname] = model

100%|██████████| 1000/1000 [00:01<00:00, 724.53it/s]

fig = results['Shin']['SIP-K389175'].plot_traces()

fig = results['Shin']['SIP-K389175'].plot_fit(discard=500)

import matplotlib.pyplot as plt
freq = np.logspace(-2, 100, 10000)
w = 2*np.pi*freq
theta = np.array([8.16E02, 3.12E03, np.log(1.80E-14), np.log(2.25E-06), 0.3098, 0.4584])

Z = model.forward(theta, w)
plt.plot(*Z)

[<matplotlib.lines.Line2D at 0x12336b978>]

Conclusions