Repository: brian-team/brian2modelfitting (language: python)

Issue:
Update tests
Resolves #55.

Few tests for `Inferencer`. 
This needs to be updated but for now I will merge this into `sbi_support` because of the deadline for GSoC.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index fe292b66..8f674e08 100644
--- a/.gitignore
+++ b/.gitignore
@@ -80,3 +80,19 @@ target/
 
 # Installation and auto-generated files
 brian2modelfitting.egg-info
+
+# sbi logs
+examples/sbi-logs/
+brian2modelfitting/tests/sbi-logs/
+
+# other sbi and torch stuff
+examples/*.pt
+examples/*.pth
+brian2modelfitting/tests/*.pt
+brian2modelfitting/tests/*.pth
+
+# numpy files
+examples/*.npy
+examples/*.npz
+brian2modelfitting/tests/*.npy
+brian2modelfitting/tests/*.npz
\ No newline at end of file
diff --git a/README.md b/README.md
index 31e29c70..278c3165 100644
--- a/README.md
+++ b/README.md
@@ -1,28 +1,38 @@
 brian2modelfitting
 ==================
 
+Model fitting toolbox for Brian 2 simulator.
 
-Model Fitting Toolbox for Brian 2 simulator, to allow the user to find the best fit of the parameters for recorded traces and spike trains.
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4601961.svg)](https://doi.org/10.5281/zenodo.4601961)
 
-By default, we support a range of global derivative-free optimization methods, that include popular methods for model fitting, such as: Differential Evolution, Particle Swarm Optimization and Covariance Matrix Adaptation (provided by the Nevergrad) as well as Bayesian Optimization for black box functions (provided by Scikit-Optimize).
 
-Documentation for Brian2modelfitting can be found at http://brian2modelfitting.readthedocs.org
 
-Brian2modelfitting is released under the terms of the CeCILL 2.1 license.
+The package [`brian2modelfitting`](https://brian2modelfitting.readthedocs.io) allows the user to find the best fit of the unknown free parameters for recorded traces and spike trains. It also supports simulation-based inference, where instead of point-estimated parameter values, a full posterior distribution over the parameters is computed.
 
-Please report issues at the github issue tracker (https://github.com/brian-team/brian2modelfitting/issues) or at the brian discussion forum (https://brian.discourse.group)
+By default, the toolbox supports a range of global derivative-free optimization methods, that include popular methods for model fitting: differential evolution, particle swarm optimization and covariance matrix adaptation (provided by the  [`Nevergrad`](https://facebookresearch.github.io/nevergrad/), a gradient-free optimization platform) as well as Bayesian optimization for black box functions (provided by [`scikit-optimize`](https://scikit-optimize.github.io/stable/), a sequential model-based optimization library). On the other hand, simulation-based inference is the process of finding parameters of a simulator from observations by taking a Bayesian approach via sequential neural posterior estimation, likelihood estimation or ration estimation (provided by the [`sbi`](https://www.mackelab.org/sbi/)), where neural densitiy estimator, a deep neural network allowing probabilistic association between the data and underlying parameter space, is trained. After the network is trained, the approximated posterior distribution is available.
 
-[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4601961.svg)](https://doi.org/10.5281/zenodo.4601961)
+Documentation
+-------------
+[![Documentation Status](https://readthedocs.org/projects/ansicolortags/badge/?version=latest)](https://brian2modelfitting.readthedocs.io)
+
+The full documentation is available at http://brian2modelfitting.readthedocs.org.
+
+Testing status
+--------------
+[![Build Status](https://github.com/brian-team/brian2modelfitting/workflows/Tests/badge.svg)](https://github.com/brian-team/brian2modelfitting/actions) 
+[![Coverage Status](https://coveralls.io/repos/github/brian-team/brian2modelfitting/badge.svg?branch=master)](https://coveralls.io/github/brian-team/brian2modelfitting?branch=master)
 
 Installation
 ------------
-Module can be installed with pypi:
+Install `brian2modelfitting` from the Python package index via pip:
 ```
 pip install brian2modelfitting
 ```
 
+License
+-------
+The model fitting toolbox is released under the terms of the CeCILL 2.1 license and is available [here](https://github.com/brian-team/brian2modelfitting/blob/master/LICENSE).
 
-Testing status
---------------
-[![Build Status](https://github.com/brian-team/brian2modelfitting/workflows/Tests/badge.svg)](https://github.com/brian-team/brian2modelfitting/actions) 
-[![Coverage Status](https://coveralls.io/repos/github/brian-team/brian2modelfitting/badge.svg?branch=master)](https://coveralls.io/github/brian-team/brian2modelfitting?branch=master)
+Use
+---
+Please report issues at the GitHub [issue tracker](https://github.com/brian-team/brian2modelfitting/issues) or at the [Brian 2 discussion forum](https://brian.discourse.group).
diff --git a/brian2modelfitting/__init__.py b/brian2modelfitting/__init__.py
index c3c8bbc1..e1994f6a 100644
--- a/brian2modelfitting/__init__.py
+++ b/brian2modelfitting/__init__.py
@@ -6,6 +6,7 @@
 
 from .tests import run as run_test
 from .fitter import *
+from .inferencer import Inferencer
 from .optimizer import *
 from .metric import *
 from .simulator import *
@@ -16,4 +17,5 @@
            'Simulator', 'RuntimeSimulator', 'CPPStandaloneSimulator',
            'MSEMetric', 'Metric', 'GammaFactor', 'FeatureMetric',
            'SpikeMetric', 'TraceMetric', 'get_gamma_factor', 'firing_rate',
-           'Fitter', 'SpikeFitter', 'TraceFitter', 'OnlineTraceFitter']
+           'Fitter', 'SpikeFitter', 'TraceFitter', 'OnlineTraceFitter',
+           'Inferencer']
diff --git a/brian2modelfitting/inferencer.py b/brian2modelfitting/inferencer.py
new file mode 100644
index 00000000..3cc4154b
--- /dev/null
+++ b/brian2modelfitting/inferencer.py
@@ -0,0 +1,1390 @@
+from numbers import Number
+from typing import Mapping
+import warnings
+
+from brian2.core.functions import Function
+from brian2.core.namespace import get_local_namespace
+from brian2.core.network import Network
+from brian2.devices.cpp_standalone.device import CPPStandaloneDevice
+from brian2.devices.device import get_device, device
+from brian2.equations.equations import Equations
+from brian2.groups.neurongroup import NeuronGroup
+from brian2.input.timedarray import TimedArray
+from brian2.monitors.spikemonitor import SpikeMonitor
+from brian2.monitors.statemonitor import StateMonitor
+from brian2.units.fundamentalunits import (DIMENSIONLESS,
+                                           fail_for_dimension_mismatch,
+                                           get_dimensions,
+                                           Quantity)
+from brian2.utils.logger import get_logger
+from brian2modelfitting.fitter import get_spikes
+import numpy as np
+from sbi.utils.get_nn_models import (posterior_nn,
+                                     likelihood_nn,
+                                     classifier_nn)
+from sbi.utils.torchutils import BoxUniform
+import sbi.analysis
+import sbi.inference
+import torch
+
+from .simulator import RuntimeSimulator, CPPStandaloneSimulator
+from .utils import tqdm
+
+
+logger = get_logger(__name__)
+
+
+def configure_simulator():
+    """Return the configured simulator, which can be either
+    `.RuntimeSimulator`, object for the use with `.RuntimeDevice`, or
+    `.CPPStandaloneSimulator`, object for the use with
+    `.CPPStandaloneDevice`.
+
+    Parameters
+    ----------
+    None
+
+    Returns
+    -------
+    brian2modelfitting.simulator.Simulator
+        Either `.RuntimeSimulator` or `.CPPStandaloneSimulator`
+        depending on the currently active ``Device`` object describing
+        the available computational engine.
+    """
+    simulators = {'CPPStandaloneDevice': CPPStandaloneSimulator(),
+                  'RuntimeDevice': RuntimeSimulator()}
+    if isinstance(get_device(), CPPStandaloneDevice):
+        if device.has_been_run is True:
+            build_options = dict(device.build_options)
+            get_device().reinit()
+            get_device().activate(**build_options)
+    return simulators[get_device().__class__.__name__]
+
+
+def get_full_namespace(additional_namespace, level=0):
+    """Return the namespace with added ``additional_namespace``, in
+    which references to external parameters or functions are stored.
+
+    Parameters
+    ----------
+    additional_namespace : dict
+        References to external parameters or functions, where key is
+        the name and value is the value of the external parameter or
+        function.
+    level : int, optional
+        How far to go back to get the locals/globals.
+
+    Returns
+    -------
+    dict
+        Namespace with additional references to the external parameters
+        or functions.
+    """
+    namespace = {key: value
+                 for key, value in get_local_namespace(level=level + 1).items()
+                 if isinstance(value, (Number, Quantity, Function))}
+    namespace.update(additional_namespace)
+    return namespace
+
+
+def get_param_dict(param_values, param_names, n_values):
+    """Return a dictionary compiled of parameter names and values.
+
+    Parameters
+    ----------
+    param_values : numpy.ndarray
+        Parameter values in a 2-dimensional array with the number of
+        rows corresponding to a number of samples and the number of
+        columns corresponding to a number of parameters.
+    param_names : list
+        List containing parameter names.
+    n_values : int
+        Total number of given values for a single parameter.
+
+    Returns
+    -------
+    dict
+        Dictionary containing key-value pairs that correspond to a
+        parameter name and value(s).
+    """
+    param_values = np.array(param_values)
+    param_dict = dict()
+    for name, value in zip(param_names, param_values.T):
+        param_dict[name] = (np.ones((n_values, )) * value)
+    return param_dict
+
+
+def calc_prior(param_names, **params):
+    """Return the prior distribution over given parameters.
+
+    Note that the only currently supported prior distribution is the
+    multi-dimensional uniform distribution defined on a box.
+
+    Parameters
+    ----------
+    param_names : list
+        List containing parameter names.
+    params : dict
+        Dictionary with keys that correspond to parameter names, and
+        the respective values are 2-element lists that hold the upper
+        and the lower bound of a distribution.
+
+    Return
+    ------
+    sbi.utils.torchutils.BoxUniform
+        ``sbi``-compatible object that contains a uniform prior
+        distribution over a given set of parameters.
+    """
+    for param_name in param_names:
+        if param_name not in params:
+            raise TypeError(f'Bounds must be set for parameter {param_name}')
+    prior_min = []
+    prior_max = []
+    for param_name in param_names:
+        prior_min.append(min(params[param_name]).item())
+        prior_max.append(max(params[param_name]).item())
+    prior = BoxUniform(low=torch.as_tensor(prior_min),
+                       high=torch.as_tensor(prior_max))
+    return prior
+
+
+class Inferencer(object):
+    """Class for a simulation-based inference.
+
+    It offers an interface similar to that of the `.Fitter` class but
+    instead of fitting, a neural density estimator is trained using a
+    generative model which ultimately provides the posterior
+    distribution over unknown free parameters.
+
+    To utilize simulation-based inference, this class uses a ``sbi``
+    library, for details see Tejero-Cantero 2020.
+
+    Parameters
+    ----------
+    dt : brian2.units.fundamentalunits.Quantity
+        Integration time step.
+    model : str or brian2.equations.equations.Equations
+        Single cell model equations.
+    input : dict
+        Input traces in the dictionary format where key corresponds to
+        the name of the input variable as defined in ``model``, and
+        value corresponds to an array of input traces.
+    output : dict
+        Dictionary of recorded (or simulated) output data traces, where
+        key corresponds to the name of the output variable as defined
+        in ``model``, and value corresponds to an array of recorded
+        traces.
+    features : dict, optional
+        Dictionary of callables that take a 1-dimensional voltage
+        trace or a spike train and output summary statistics. Keys
+        correspond to output variable names, while values are lists of
+        callables. If ``features`` are set to None, automatic feature
+        extraction process will occur instead either by using the
+        default multi-layer perceptron or by using the custom embedding
+        network.
+    method : str, optional
+        Integration method.
+    threshold : str, optional
+        The condition which produces spikes. It should be a single line
+        boolean expression.
+    reset : str, optional
+        The (possibly multi-line) string with the code to execute on
+        reset.
+    refractory : bool or str, optional
+        Either the length of the refractory period (e.g., ``2*ms``), a
+        string expression that evaluates to the length of the
+        refractory period after each spike, e.g., ``'(1 + rand())*ms'``,
+        or a string expression evaluating to a boolean value, given the
+        condition under which the neuron stays refractory after a spike,
+        e.g., ``'v > -20*mV'``.
+    param_init : dict, optional
+        Dictionary of state variables to be initialized with respective
+        values, i.e., initial conditions for the model.
+
+    References
+    ----------
+    * Tejero-Cantero, A., Boelts, J. et al. "sbi: A toolkit for
+      simulation-based inference" Journal of Open Source Software
+      (JOOS), 5(52):2505. 2020.
+    """
+    def __init__(self, dt, model, input, output, features=None, method=None,
+                 threshold=None, reset=None, refractory=False,
+                 param_init=None):
+        # time scale
+        self.dt = dt
+
+        # model equations
+        if isinstance(model, str):
+            model = Equations(model)
+        else:
+            raise TypeError('Equations must be appropriately formatted.')
+
+        # input data traces
+        if not isinstance(input, Mapping):
+            raise TypeError('`input` argument must be a dictionary mapping'
+                            ' the name of the input variable and `input`.')
+        if len(input) > 1:
+            raise NotImplementedError('Only a single input is supported.')
+        input_var = list(input.keys())[0]
+        input = input[input_var]
+        if input_var not in model.identifiers:
+            raise NameError(f'{input_var} is not an identifier in the model.')
+
+        # output data traces
+        if not isinstance(output, Mapping):
+            raise TypeError('`output` argument must be a dictionary mapping'
+                            ' the name of the output variable and `output`.')
+        output_var = list(output.keys())
+        output = list(output.values())
+        for o_var in output_var:
+            if o_var != 'spikes' and o_var not in model.names:
+                raise NameError(f'{o_var} is not a model variable')
+        self.output_var = output_var
+        self.output = output
+
+        # create variable for parameter names
+        self.param_names = sorted(model.parameter_names)
+
+        # set the simulation time for a given time scale
+        self.n_traces, n_steps = input.shape
+        self.sim_time = dt * n_steps
+
+        # handle multiple output variables
+        self.output_dim = []
+        for o_var, out in zip(self.output_var, self.output):
+            if o_var == 'spikes':
+                self.output_dim.append(DIMENSIONLESS)
+            else:
+                self.output_dim.append(model[o_var].dim)
+                fail_for_dimension_mismatch(out, self.output_dim[-1],
+                                            'The provided target values must'
+                                            ' have the same units as the'
+                                            f' variable {o_var}')
+
+        # add input to equations
+        self.model = model
+        input_dim = get_dimensions(input)
+        input_dim = '1' if input_dim is DIMENSIONLESS else repr(input_dim)
+        input_eqs = f'{input_var} = input_var(t, i % n_traces) : {input_dim}'
+        self.model += input_eqs
+
+        # add output to equations
+        counter = 0
+        for o_var, o_dim in zip(self.output_var, self.output_dim):
+            if o_var != 'spikes':
+                counter += 1
+                output_expr = f'output_var_{counter}(t, i % n_traces)'
+                output_dim = ('1' if o_dim is DIMENSIONLESS else repr(o_dim))
+                output_eqs = f'{o_var}_target = {output_expr} : {output_dim}'
+                self.model += output_eqs
+
+        # create ``TimedArray`` object for input w.r.t. a given time scale
+        self.input_traces = TimedArray(input.transpose(), dt=self.dt)
+
+        # handle initial values for the ODE system
+        if not param_init:
+            param_init = {}
+        for param in param_init.keys():
+            if not (param in self.model.diff_eq_names or
+                    param in self.model.parameter_names):
+                raise ValueError(f'{param} is not a model variable or a'
+                                 ' parameter in the model')
+        self.param_init = param_init
+
+        # handle the rest of optional parameters for the ``NeuronGroup`` class
+        self.method = method
+        self.threshold = threshold
+        self.reset = reset
+        self.refractory = refractory
+
+        # observation the focus is on
+        obs = []
+        if features:
+            for ov, o in zip(self.output_var, self.output):
+                for _o in o:
+                    for feature in features[ov]:
+                        obs.append(feature(_o))
+            x_o = np.array(obs, dtype=np.float32)
+        else:
+            for o in self.output:
+                o = np.array(o)
+                obs.append(o.ravel().astype(np.float32))
+            x_o = np.concatenate(obs)
+        self.x_o = x_o
+        self.features = features
+
+        # additional placeholders
+        self.params = None
+        self.n_samples = None
+        self.samples = None
+        self.posterior = None
+        self.theta = None
+        self.x = None
+        self.sbi_device = 'cpu'
+
+    @property
+    def n_neurons(self):
+        """Return the number of neurons that are used in `.NeuronGroup`
+        class while generating data for training the neural density
+        estimator.
+
+        Unlike the `.Fitter` class, `.Inferencer` does not take the
+        total number of samples directly in the constructor. Thus, this
+        property becomes available only after the simulation is performed.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        int
+            Total number of neurons.
+        """
+        if self.n_samples is None:
+            raise ValueError('Number of samples have not been yet defined.'
+                             ' Call `generate_training_data` method first.')
+        return self.n_traces * self.n_samples
+
+    def setup_simulator(self, network_name, n_neurons, output_var, param_init,
+                        level=1):
+        """Return configured simulator.
+
+        Parameters
+        ----------
+        network_name : str
+            Network name.
+        n_neurons : int
+            Number of neurons which equals to the number of samples
+            times the number of input/output traces.
+        output_var : str
+            Name of the output variable.
+        param_init : dict
+            Dictionary of state variables to be initialized with
+            respective values.
+        level : int, optional
+            How far to go back to get the locals/globals.
+
+        Returns
+        -------
+        brian2modelfitting.simulator.Simulator
+            Configured simulator w.r.t. the available device.
+        """
+        # configure the simulator
+        simulator = configure_simulator()
+
+        # update the local namespace
+        namespace = get_full_namespace({'input_var': self.input_traces,
+                                        'n_traces': self.n_traces},
+                                       level=level+1)
+        counter = 0
+        for o_var, out in zip(self.output_var, self.output):
+            if o_var != 'spikes':
+                counter += 1
+                namespace[f'output_var_{counter}'] = TimedArray(out.T,
+                                                                dt=self.dt)
+
+        # setup neuron group
+        kwds = {}
+        if self.method is not None:
+            kwds['method'] = self.method
+        model = (self.model
+                 + Equations('iteration : integer (constant, shared)'))
+        neurons = NeuronGroup(N=n_neurons,
+                              model=model,
+                              threshold=self.threshold,
+                              reset=self.reset,
+                              refractory=self.refractory,
+                              dt=self.dt,
+                              namespace=namespace,
+                              name='neurons',
+                              **kwds)
+
+        # create a network of neurons
+        network = Network(neurons)
+        if isinstance(output_var, str):
+            output_var = [output_var]
+        if 'spikes' in output_var:
+            network.add(SpikeMonitor(neurons, name='spikemonitor'))
+        record_vars = [v for v in output_var if v != 'spikes']
+        if len(record_vars):
+            network.add(StateMonitor(source=neurons, variables=record_vars,
+                                     record=True, dt=self.dt,
+                                     name='statemonitor'))
+
+        # initialize the simulator
+        simulator.initialize(network, param_init, name=network_name)
+        return simulator
+
+    def init_prior(self, **params):
+        """Return the prior uniform distribution over the parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary with keys that correspond to parameter names,
+            and the respective values are 2-element lists that hold
+            the upper and the lower bound of a distribution.
+
+        Returns
+        -------
+        sbi.utils.torchutils.BoxUniform
+            ``sbi``-compatible object that contains a uniform prior
+            distribution over a given set of parameters.
+        """
+        for param in params:
+            if param not in self.param_names:
+                raise ValueError(f'Parameter {param} must be defined as a'
+                                 ' model\'s parameter')
+        prior = calc_prior(self.param_names, **params)
+        self.params = params
+        return prior
+
+    def generate_training_data(self, n_samples, prior):
+        """Return sampled prior given the total number of samples.
+
+        Parameter
+        ---------
+        n_samples : int
+            The number of samples.
+        prior : sbi.utils.BoxUniform
+            Uniformly distributed prior over given parameters.
+
+        Returns
+        -------
+        numpy.ndarray
+            Sampled prior with the number of rows that corresponds to
+            the ``n_samples``, while the number of columns depends on
+            the number of free parameters.
+        """
+        # set n_samples to class variable to be able to call self.n_neurons
+        self.n_samples = n_samples
+
+        # sample from the prior distribution
+        theta = prior.sample((n_samples, ))
+        theta = np.atleast_2d(theta.numpy())
+        return theta
+
+    def extract_summary_statistics(self, theta, level=0):
+        """Return the summary statistics for the process of training
+        of the neural density estimator.
+
+        Parameters
+        ----------
+        theta : numpy.ndarray
+            Sampled prior with ``n_samples`` rows, and the number of
+            columns corresponds to the number of free parameters.
+        level : int, optional
+            How far to go back to get the locals/globals.
+
+        Returns
+        -------
+        numpy.ndarray
+            Summary statistics.
+        """
+        # repeat each row for how many input/output different trace are there
+        _theta = np.repeat(theta, repeats=self.n_traces, axis=0)
+
+        # create a dictionary with repeated sampled prior
+        d_param = get_param_dict(_theta, self.param_names, self.n_neurons)
+
+        # set up and run the simulator
+        network_name = 'infere'
+        simulator = self.setup_simulator(network_name=network_name,
+                                         n_neurons=self.n_neurons,
+                                         output_var=self.output_var,
+                                         param_init=self.param_init,
+                                         level=level+1)
+        simulator.run(self.sim_time, d_param, self.param_names, iteration=0,
+                      name=network_name)
+
+        # extract features for each output variable and each trace
+        try:
+            obs = simulator.statemonitor.recorded_variables
+        except KeyError:
+            logger.warn('The state monitor object is not defined.',
+                        name_suffix='statemonitor_definition')
+        if 'spikes' in self.output_var:
+            spike_trains = list(simulator.spikemonitor.spike_trains().values())
+        x = []
+        if self.features:
+            for ov in self.output_var:
+                print(f'Extracting features for \'{ov}\'...')
+                summary_statistics = []
+                if ov != 'spikes':
+                    o = obs[ov].get_value().T
+                # TODO: should be vectorized
+                for i in range(self.n_neurons):
+                    if ov != 'spikes':
+                        for feature in self.features[ov]:
+                            summary_statistics.append(feature(o[i, :]))
+                    else:
+                        for feature in self.features[ov]:
+                            summary_statistics.append(feature(spike_trains[i]))
+                _x = np.array(summary_statistics, dtype=np.float32)
+                _x = _x.reshape(self.n_samples, -1)
+                x.append(_x)
+            x = np.hstack(x)
+        else:
+            print('Aranging traces for automatic feature extraction...')
+            for ov in self.output_var:
+                o = obs[ov].get_value().T
+                x.append(o.reshape(self.n_samples, -1).astype(np.float32))
+            x = np.hstack(x)
+        return x
+
+    def save_summary_statistics(self, f, theta=None, x=None):
+        """Save sampled prior and the extracted summary statistics into
+        a single compressed ``.npz`` file.
+
+        Parameters
+        ----------
+        f : str or os.PathLike
+            Path to a file either as string or ``os.PathLike`` object
+            that contains file name.
+        theta : numpy.ndarray, optional
+            Sampled prior.
+        x : numpy.ndarray, optional
+            Summary statistics.
+
+        Returns
+        -------
+        None
+        """
+        if theta is not None:
+            t = theta
+        elif self.theta is not None:
+            t = self.theta
+        else:
+            raise AttributeError('Provide sampled prior or call'
+                                 ' `infer_step` method first.')
+        if x is not None:
+            pass
+        elif self.x is not None:
+            x = self.x
+        else:
+            raise AttributeError('Provide summary feautures or call '
+                                 ' `infer_step` method first.')
+        np.savez_compressed(f, theta=t, x=x)
+
+    def load_summary_statistics(self, f):
+        """Load samples from a prior and the extracted summary
+        statistics from a compressed ``.npz`` file.
+
+        Parameters
+        ----------
+        f : str or os.PathLike
+            Path to file either as string or ``os.PathLike`` object
+            that contains file name.
+
+        Returns
+        -------
+        tuple
+            Sampled prior and the summary statistics arrays.
+        """
+        loaded = np.load(f, allow_pickle=True)
+        if set(loaded.files) == {'theta', 'x'}:
+            theta = loaded['theta']
+            x = loaded['x']
+        self.theta = theta
+        self.x = x
+        return theta, x
+
+    def init_inference(self, inference_method, density_estimator_model, prior,
+                       sbi_device='cpu', **inference_kwargs):
+        """Return instantiated inference object.
+
+        Parameters
+        ----------
+        inference_method : str
+            Inference method. Either SNPE, SNLE or SNRE.
+        density_estimator_model : str
+            The type of density estimator to be created. Either
+            ``mdn``, ``made``, ``maf``, ``nsf`` for SNPE and SNLE, or
+            ``linear``, ``mlp``, ``resnet`` for SNRE.
+        prior : sbi.utils.BoxUniform
+            Uniformly distributed prior over free parameters.
+        sbi_device : str, optional
+            Device on which the ``sbi`` will operate. By default this
+            is set to ``cpu`` and it is advisable to remain so for most
+            cases. In cases where the user provides custom embedding
+            network through ``inference_kwargs`` argument, which will
+            be trained more efficiently by using GPU, device should be
+            set accordingly to either ``gpu`` or ``cuda``.
+        inference_kwargs : dict, optional
+            Additional keyword arguments for different density
+            estimator builder functions:
+            ``sbi.utils.get_nn_models.posterior_nn`` for SNPE,
+            ``sbi.utils.get_nn_models.classifier_nn`` for SNRE, and
+            ``sbi.utils.get_nn_models.likelihood_nn`` for SNLE. For
+            details check the official ``sbi`` documentation. A single
+            highlighted keyword argument is a custom embedding network
+            that serves a purpose to learn features from potentially
+            high-dimensional simulation outputs. By default multi-layer
+            perceptron is used if no custom embedding network is
+            provided. For SNPE and SNLE, the user may pass an embedding
+            network for simulation outputs through ``embedding_net``
+            argument, while for SNRE, the user may pass two embedding
+            networks, one for parameters through
+            ``embedding_net_theta`` argument, and the other for
+            simulation outputs through ``embedding_net_x`` argument.
+
+        Returns
+        -------
+        sbi.inference.NeuralInference
+            Instantiated inference object.
+        """
+        try:
+            inference_method = str.upper(inference_method)
+            inference_method_fun = getattr(sbi.inference, inference_method)
+        except AttributeError:
+            raise NameError(f'Inference method {inference_method} is not'
+                            ' supported. Choose between SNPE, SNLE or SNRE.')
+        finally:
+            if inference_method == 'SNPE':
+                density_estimator_builder = posterior_nn(
+                    model=density_estimator_model, **inference_kwargs)
+            elif inference_method == 'SNLE':
+                density_estimator_builder = likelihood_nn(
+                    model=density_estimator_model, **inference_kwargs)
+            else:
+                density_estimator_builder = classifier_nn(
+                    model=density_estimator_model, **inference_kwargs)
+
+        sbi_device = str.lower(sbi_device)
+        if sbi_device in ['cuda', 'gpu']:
+            if torch.cuda.is_available():
+                sbi_device = 'gpu'
+                self.sbi_device = 'cuda'
+            else:
+                logger.warn(f'Device {sbi_device} is not available.'
+                            ' Falling back to CPU.')
+                sbi_device = 'cpu'
+                self.sbi_device = sbi_device
+        else:
+            sbi_device = 'cpu'
+            self.sbi_device = sbi_device
+
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore')
+            inference = inference_method_fun(prior, density_estimator_builder,
+                                             device=sbi_device,
+                                             show_progress_bars=True)
+        self.inference = inference
+        return inference
+
+    def train(self, inference, theta, x, *args, **train_kwargs):
+        """Return the trained neural inference object.
+
+        Parameters
+        ----------
+        inference : sbi.inference.NeuralInference
+            Instantiated inference object with stored paramaters and
+            simulation outputs prepared for the training process.
+        theta : torch.tensor
+            Sampled prior.
+        x : torch.tensor
+            Summary statistics.
+        args : tuple, optional
+            Contains a uniformly distributed proposal. Used only for
+            SNPE, for SNLE and SNRE, proposal should not be passed to
+            inference object, thus ``args`` should not be passed. The
+            additional arguments should be passed only if the
+            parameters were not sampled from the prior, e.g., during
+            the multi-round inference. For SNLE and SNRE, this can be
+            the number of round from which the data is stemmed from,
+            e.g., 0 means from the prior. This is used only if the
+            ``discard_prior_samples`` is set to True inside the
+            ``train_kwargs``.
+        train_kwargs : dict, optional
+            Additional keyword arguments for ``train`` method in the
+            ``sbi.inference.NeuralInference`` class. The user is able
+            to gain the full control over the training process by
+            tuning hyperparameters, e.g., the batch size (by specifiying
+            ``training_batch_size`` argument), the learning rate
+            (``learning_rate``), the validation fraction
+            (``validation_fraction``), the number of training epochs
+            (``max_num_epochs``), etc. For details, check the official
+            ``sbi`` documentation.
+
+        Returns
+        -------
+        sbi.inference.NeuralInference
+            Trained inference object.
+        """
+        inference = inference.append_simulations(theta, x, *args)
+        _ = inference.train(**train_kwargs)
+        return inference
+
+    def build_posterior(self, inference, **posterior_kwargs):
+        """Return the updated inference and the neural posterior
+        objects.
+
+        Parameters
+        ----------
+        inference : sbi.inference.NeuralInference
+            Instantiated inference object with stored paramaters and
+            simulation outputs.
+        posterior_kwargs : dict, optional
+            Additional keyword arguments for ``build_posterior`` method
+            in all ``sbi.inference.NeuralInference``-type classes. For
+            details, check the official ``sbi`` documentation.
+
+        Returns
+        -------
+        tuple
+            ``sbi.inference.NeuralInference`` object with stored
+            paramaters and simulation outputs prepared for training and
+            the neural posterior object.
+        """
+        posterior = inference.build_posterior(**posterior_kwargs)
+        return inference, posterior
+
+    def infer_step(self, proposal, inference, n_samples=None, theta=None,
+                   x=None, train_kwargs={}, posterior_kwargs={}, *args):
+        """Return the trained neural density estimator.
+
+        Parameters
+        ----------
+        proposal : ``sbi.utils.torchutils.BoxUniform``
+            Prior over parameters for the current round of inference.
+        inference : ``sbi.inference.NeuralInference``
+            Inference object obtained via `.init_inference` method.
+        n_samples : int, optional
+            The number of samples.
+        theta : numpy.ndarray, optional
+            Sampled prior.
+        x : numpy.ndarray, optional
+            Summary statistics.
+        train_kwargs : dict, optional
+            Additional keyword arguments for `.train`.
+        posterior_kwargs : dict, optional
+            Additional keyword arguments for `.build_posterior`.
+        args : list, optional
+            Additional arguments for `.train`.
+
+        Returns
+        -------
+        sbi.inference.posteriors.base_posterior.NeuralPosterior
+            Trained posterior.
+        """
+        # extract the training data and make adjustments for the ``sbi``
+        if theta is None:
+            if n_samples is None:
+                raise ValueError('Either provide `theta` or `n_samples`.')
+            else:
+                theta = self.generate_training_data(n_samples, proposal)
+        self.theta = theta
+        theta = torch.tensor(theta, dtype=torch.float32)
+
+        # extract the summary statistics and make adjustments for the ``sbi``
+        if x is None:
+            if n_samples is None:
+                raise ValueError('Either provide `x` or `n_samples`.')
+            else:
+                x = self.extract_summary_statistics(theta, level=2)
+        self.x = x
+        x = torch.tensor(x)
+
+        # pass the simulated data to the inference object and train it
+        inference = self.train(inference, theta, x, *args, **train_kwargs)
+
+        # use the density estimator to build the posterior
+        inference, posterior = self.build_posterior(inference,
+                                                    **posterior_kwargs)
+        self.inference = inference
+        self.posterior = posterior
+        return posterior
+
+    def infer(self, n_samples=None, theta=None, x=None, n_rounds=1,
+              inference_method='SNPE', density_estimator_model='maf',
+              inference_kwargs={}, train_kwargs={}, posterior_kwargs={},
+              restart=False, sbi_device='cpu', **params):
+        """Return the trained posterior.
+
+        Note that if ``theta`` and ``x`` are not provided,
+        ``n_samples`` has to be defined. Otherwise, if ``n_samples`` is
+        provided, neither ``theta`` nor ``x`` are needed and will be
+        ignored.
+
+        Parameters
+        ----------
+        n_samples : int, optional
+            The number of samples.
+        theta : numpy.ndarray, optional
+            Sampled prior.
+        x : numpy.ndarray, optional
+            Summary statistics.
+        n_rounds : int, optional
+            If ``n_rounds`` is set to 1, amortized inference will be
+            performed. Otherwise, if ``n_rounds`` is integer larger
+            than 1, multi-round inference will be performed. This is
+            only valid if posterior has not been defined manually.
+            Otherwise, if this method is called after posterior has
+            already been built, multi-round inference is performed.
+        inference_method : str, optional
+            Inference method. Either SNPE, SNLE or SNRE.
+        density_estimator_model : str, optional
+            The type of density estimator to be created. Either
+            ``mdn``, ``made``, ``maf``, ``nsf`` for SNPE and SNLE, or
+            ``linear``, ``mlp``, ``resnet`` for SNRE.
+        inference_kwargs : dict, optional
+            Additional keyword arguments for the `.init_inference`.
+        train_kwargs : dict, optional
+            Additional keyword arguments for `.train`.
+        posterior_kwargs : dict, optional
+            Additional keyword arguments for `.build_posterior`.
+        restart : bool, optional
+            When the method is called for a second time, set to True if
+            amortized inference should be performed. If False,
+            multi-round inference with the existing posterior will be
+            performed.
+        sbi_device : str, optional
+            Device on which the ``sbi`` will operate. By default this
+            is set to ``cpu`` and it is advisable to remain so for most
+            cases. In cases where the user provide custom embedding
+            network through ``inference_kwargs`` argument, which will
+            be trained more efficiently by using GPU, device should be
+            set accordingly to either ``gpu`` or ``cuda``.
+        params : dict
+            Bounds for each parameter. Keys should correspond to names
+            of parameters as defined in the model equations, while
+            values are lists with the lower and the upper bound with
+            corresponding quantities of the parameter.
+
+        Returns
+        -------
+        sbi.inference.posteriors.base_posterior.NeuralPosterior
+            Approximated posterior distribution over parameters.
+        """
+        if restart:
+            self.posterior = None
+        if self.posterior is None:
+            # handle the number of rounds
+            if not isinstance(n_rounds, int):
+                raise ValueError('`n_rounds` has to be a positive integer.')
+
+            # handle inference methods
+            try:
+                inference_method = str.upper(inference_method)
+            except ValueError as e:
+                print(e, '\nInference method should be defined as string.')
+            if inference_method not in ['SNPE', 'SNLE', 'SNRE']:
+                raise ValueError(f'Inference method {inference_method} is not'
+                                 ' supported.')
+
+            # initialize prior
+            prior = self.init_prior(**params)
+
+            # extract the training data
+            if theta is None:
+                if n_samples is None:
+                    raise ValueError('Either provide `theta` or `n_samples`.')
+                else:
+                    theta = self.generate_training_data(n_samples, prior)
+            self.theta = theta
+
+            # extract the summary statistics
+            if x is None:
+                if n_samples is None:
+                    raise ValueError('Either provide `x` or `n_samples`.')
+                else:
+                    x = self.extract_summary_statistics(theta, level=1)
+            self.x = x
+
+            # initialize inference object
+            _ = self.init_inference(inference_method, density_estimator_model,
+                                    prior, sbi_device, **inference_kwargs)
+
+            # args for SNPE in `.train`
+            # args for SNRE and SNLE are not supported here, if needed the user
+            # could provide them by using more flexible inferface via
+            # `.infer_step` method
+            if inference_method == 'SNPE':
+                args = [prior, ]
+            else:
+                args = [None, ]
+        else:  # `.infer_step` has been called manually
+            x_o = torch.tensor(self.x_o, dtype=torch.float32)
+            prior = self.posterior.set_default_x(x_o)
+            if self.posterior._method_family == 'snpe':
+                args = [prior, ]
+            else:
+                args = [None, ]
+            # generate data if the posterior already exist given proposal
+            self.theta = self.generate_training_data(self.n_samples, prior)
+            self.x = self.extract_summary_statistics(self.theta, level=1)
+
+        # allocate empty list of posteriors
+        posteriors = []
+
+        # set a proposal
+        proposal = prior
+
+        # main inference loop
+        if self.posterior or n_rounds > 1:
+            tqdm_desc = f'{n_rounds}-round focused inference'
+        else:
+            tqdm_desc = 'Amortized inference'
+        for round in tqdm(range(n_rounds), desc=tqdm_desc):
+            # inference step
+            posterior = self.infer_step(proposal, self.inference,
+                                        n_samples, self.theta, self.x,
+                                        train_kwargs, posterior_kwargs, *args)
+
+            # append the current posterior to the list of posteriors
+            posteriors.append(posterior)
+
+            # update the proposal given the observation
+            if n_rounds > 1 and round < n_rounds - 1:
+                x_o = torch.tensor(self.x_o, dtype=torch.float32)
+                proposal = posterior.set_default_x(x_o)
+                if posterior._method_family == 'snpe':
+                    args = [proposal, ]
+                else:
+                    args = [None, ]
+                self.theta = self.generate_training_data(self.n_samples,
+                                                         proposal)
+                self.x = self.extract_summary_statistics(self.theta, level=1)
+
+        self.posterior = posterior
+        return posterior
+
+    def save_posterior(self, f):
+        """Save the density estimator state dictionary to a disk file.
+
+        Parameters
+        ----------
+        posterior: neural posterior object, optional
+            Posterior distribution over parameters.
+        f : str or os.PathLike
+            Path to file either as string or ``os.PathLike`` object
+            that contains file name.
+
+        Returns
+        -------
+        None
+        """
+        torch.save(self.posterior, f)
+
+    def load_posterior(self, f):
+        """Loads the density estimator state dictionary from a disk
+        file.
+
+        Parameters
+        ----------
+        f : str or os.PathLike
+            Path to file either as string or ``os.PathLike`` object
+            that contains file name.
+
+        Returns
+        -------
+        sbi.inference.posteriors.base_posterior.NeuralPosterior
+            Loaded neural posterior with defined method family, density
+            estimator state dictionary, the prior over parameters and
+            the output shape of the simulator.
+        """
+        p = torch.load(f)
+        self.posterior = p
+        return p
+
+    def sample(self, shape, posterior=None, **kwargs):
+        """Return samples from posterior distribution.
+
+        Parameters
+        ----------
+        shape : tuple
+            Desired shape of samples that are drawn from posterior.
+        posterior: neural posterior object, optional
+            Posterior distribution.
+        kwargs : dict, optional
+            Additional keyword arguments for ``sample`` method of
+            the neural posterior object.
+        Returns
+        -------
+        numpy.ndarray
+            Samples taken from the posterior of the shape as given in
+            ``shape``.
+        """
+        if posterior:
+            p = posterior
+        elif posterior is None and self.posterior:
+            p = self.posterior
+        else:
+            raise ValueError('Need to provide posterior argument if no'
+                             ' posterior has computed built yet.')
+        x_o = torch.tensor(self.x_o, dtype=torch.float32)
+        samples = p.sample(shape, x=x_o, **kwargs)
+        self.samples = samples
+        return samples.numpy()
+
+    def pairplot(self, samples=None, points=None, limits=None, subset=None,
+                 labels=None, ticks=None, **kwargs):
+        """Plot samples in a 2-dimensional grid with marginals and
+        pairwise marginals.
+
+        Check ``sbi.analysis.plot.pairplot`` for more details.
+
+        Parameters
+        ----------
+        samples : list or numpy.ndarray, optional
+            Samples used to build the pairplot.
+        points : dict, optional
+            Additional points to scatter, e.g., true parameter values,
+            if known.
+        limits : dict, optional
+            Limits for each parameter. Keys correspond to parameter
+            names as defined in the model, while values are lists with
+            limits defined as the Brian 2 quantity objects. If None,
+            min and max of the given samples will be used.
+        subset : list, optional
+            The names as strings of parameters to plot.
+        labels : dict, optional
+            Names for each parameter. Keys correspond to parameter
+            names as defined in the model, while values are lists of
+            strings.
+        ticks : dict, optional
+            Position of the ticks. Keys correspond to parameter names
+            as defined in the model, while values are lists with ticks
+            defined as the Brian 2 quantity objects. If None, default
+            ticks positions will be used.
+        kwargs : dict, optional
+            Additional keyword arguments for the
+            ``sbi.analysis.pairplot`` function.
+
+        Returns
+        -------
+        tuple
+            Figure and axis of the posterior distribution plot.
+        """
+        if samples is not None:
+            s = samples
+        else:
+            if self.samples is not None:
+                s = self.samples
+            else:
+                raise ValueError('Samples are not available.')
+        if points:
+            for param_name in points.keys():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            points = np.array([[points[param_name].item()
+                               for param_name in self.param_names]])
+        if limits:
+            for param_name, lim_vals in limits.items():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+                if len(lim_vals) != 2:
+                    raise ValueError('Invalid limits for parameter: '
+                                     f'{param_name}')
+            limits = [[limits[param_name][0].item(),
+                       limits[param_name][1].item()]
+                      for param_name in self.param_names]
+        if subset:
+            for param_name in subset:
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            subset = [self.param_names.index(param_name)
+                      for param_name in subset]
+        if ticks:
+            for param_name, lim_vals in ticks.items():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+                if len(lim_vals) != 2:
+                    raise ValueError('Invalid limits for parameter: '
+                                     f'{param_name}')
+            ticks = [[ticks[param_name][0].item(),
+                      ticks[param_name][1].item()]
+                     for param_name in self.param_names]
+        else:
+            ticks = []
+        if labels:
+            for param_name in labels.keys():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            labels = [labels[param_name] for param_name in self.param_names]
+        fig, axes = sbi.analysis.pairplot(samples=s,
+                                          points=points,
+                                          limits=limits,
+                                          subset=subset,
+                                          labels=labels,
+                                          ticks=ticks,
+                                          **kwargs)
+        return fig, axes
+
+    def conditional_pairplot(self, condition, density=None, points=None,
+                             limits=None, subset=None, labels=None,
+                             ticks=None, **kwargs):
+        """Plot conditional distribution given all other parameters.
+
+        The conditionals can be interpreted as slices through the
+        density at a location given by condition.
+
+        Check ``sbi.analysis.conditional_pairplot`` for more details.
+
+        Parameters
+        ----------
+        condition : numpy.ndarray
+            Condition that all but the one/two regarded parameters are
+            fixed to.
+        density : neural posterior object, optional
+            Posterior probability density.
+        points : dict, optional
+            Additional points to scatter, e.g., true parameter values,
+            if known.
+        limits : dict, optional
+            Limits for each parameter. Keys correspond to parameter
+            names as defined in the model, while values are lists with
+            limits defined as the Brian 2 quantity objects. If None,
+            min and max of the given samples will be used.
+        subset : list, optional
+            The names as strings of parameters to plot.
+        labels : dict, optional
+            Names for each parameter. Keys correspond to parameter
+            names as defined in the model, while values are lists of
+            strings.
+        ticks : dict, optional
+            Position of the ticks. Keys correspond to parameter names
+            as defined in the model, while values are lists with ticks
+            defined as the Brian 2 quantity objects. If None, default
+            ticks positions will be used.
+        kwargs : dict, optional
+            Additional keyword arguments for the
+            ``sbi.analysis.conditional_pairplot`` function.
+
+        Returns
+        -------
+        tuple
+            Figure and axis of conditional pairplot.
+        """
+        condition = torch.from_numpy(condition)
+        if density is not None:
+            d = density
+        else:
+            if self.posterior is not None:
+                if self.posterior.default_x is None:
+                    x_o = torch.tensor(self.x_o, dtype=torch.float32)
+                    d = self.posterior.set_default_x(x_o)
+                else:
+                    d = self.posterior
+            else:
+                raise ValueError('Density is not available.')
+        if points:
+            for param_name in points.keys():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            points = np.array([[points[param_name].item()
+                               for param_name in self.param_names]])
+        if limits:
+            for param_name, lim_vals in limits.items():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+                if len(lim_vals) != 2:
+                    raise ValueError('Invalid limits for parameter: '
+                                     f'{param_name}')
+            limits = [[limits[param_name][0].item(),
+                       limits[param_name][1].item()]
+                      for param_name in self.param_names]
+        else:
+            limits = [[self.params[param_name][0].item(),
+                       self.params[param_name][1].item()]
+                      for param_name in self.param_names]
+        if subset:
+            for param_name in subset:
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            subset = [self.param_names.index(param_name)
+                      for param_name in subset]
+        if ticks:
+            for param_name, lim_vals in ticks.items():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+                if len(lim_vals) != 2:
+                    raise ValueError('Invalid limits for parameter: '
+                                     f'{param_name}')
+            ticks = [[ticks[param_name][0].item(),
+                      ticks[param_name][1].item()]
+                     for param_name in self.param_names]
+        else:
+            ticks = []
+        if labels:
+            for param_name in labels.keys():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            labels = [labels[param_name] for param_name in self.param_names]
+        fig, axes = sbi.analysis.conditional_pairplot(density=d,
+                                                      condition=condition,
+                                                      limits=limits,
+                                                      points=points,
+                                                      subset=subset,
+                                                      labels=labels,
+                                                      ticks=ticks,
+                                                      **kwargs)
+        return fig, axes
+
+    def conditional_corrcoeff(self, condition, density=None, limits=None,
+                              subset=None, **kwargs):
+        """Return the conditional correlation matrix of a distribution.
+
+        All but two parameters are conditioned with the condition as
+        defined in the ``condition`` argument and the Pearson
+        correlation coefficient is computed between the remaining two
+        parameters under the distribution. This is performed for all
+        pairs of parameters given whose names are defined in the
+        ``subset`` argument. The conditional correlation matrix is
+        stored in the 2-dimenstional array.
+
+        Check ``sbi.analysis.conditional_density.conditional_corrcoeff``
+        for more details.
+
+        Parameters
+        ----------
+        condition : numpy.ndarray
+            Condition that all but the one/two regarded parameters are
+            fixed to.
+        density : neural posterior object, optional
+            Posterior probability density.
+        limits : dict, optional
+            Limits for each parameter. Keys correspond to parameter
+            names as defined in the model, while values are lists with
+            limits defined as the Brian 2 quantity objects. If None,
+            min and max of the given samples will be used.
+        subset : list, optional
+            Parameters that are taken for conditional distribution, if
+            None all parameters are considered.
+        kwargs : dict, optional
+            Additional keyword arguments for the
+            ``sbi.analysis.conditional_corrcoeff`` function.
+
+        Returns
+        -------
+        numpy.ndarray
+            Average conditional correlation matrix.
+        """
+        condition = torch.from_numpy(condition)
+        if density is not None:
+            d = density
+        else:
+            if self.posterior is not None:
+                if self.posterior.default_x is None:
+                    x_o = torch.tensor(self.x_o, dtype=torch.float32)
+                    d = self.posterior.set_default_x(x_o)
+                else:
+                    d = self.posterior
+            else:
+                raise ValueError('Density is not available.')
+        if limits:
+            for param_name, lim_vals in limits.items():
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+                if len(lim_vals) != 2:
+                    raise ValueError('Invalid limits for parameter: '
+                                     f'{param_name}')
+            limits = [[limits[param_name][0].item(),
+                       limits[param_name][1].item()]
+                      for param_name in self.param_names]
+        else:
+            limits = [[self.params[param_name][0].item(),
+                       self.params[param_name][1].item()]
+                      for param_name in self.param_names]
+        limits = torch.tensor(limits)
+        if subset:
+            for param_name in subset:
+                if param_name not in self.param_names:
+                    raise AttributeError(f'Invalid parameter: {param_name}')
+            subset = [self.param_names.index(param_name)
+                      for param_name in subset]
+        cond_coeff = sbi.analysis.conditional_corrcoeff(density=d,
+                                                        limits=limits,
+                                                        condition=condition,
+                                                        subset=subset,
+                                                        **kwargs)
+        return cond_coeff.numpy()
+
+    def generate_traces(self, n_samples=1, posterior=None, output_var=None,
+                        param_init=None, level=0):
+        """Generates traces for a single drawn sample from the trained
+        posterior and all inputs.
+
+        Parameters
+        ----------
+        n_samples : int, optional
+            The number of parameters samples. If n_samples is larger
+            than 1, the mean value of sampled set will be used.
+        posterior : neural posterior object, optional
+            Posterior distribution.
+        output_var : str or list
+            Name of the output variable to be monitored, it can also be
+            a list of names to record multiple variables.
+        param_init : dict
+            Dictionary of initial values for the model.
+        level : int, optional
+            How far to go back to get the locals/globals.
+
+        Returns
+        -------
+        brian2.units.fundamentalunits.Quantity or dict
+            If a single output variable is observed, 2-dimensional
+            array of traces generated by using a set of parameters
+            sampled from the trained posterior distribution of shape
+            (``self.n_traces``, number of time steps). Otherwise, a
+            dictionary with keys set to names of output variables, and
+            values to generated traces of respective output variables.
+        """
+        if not isinstance(n_samples, int):
+            raise ValueError('`n_samples` argument needs to be a positive'
+                             ' integer.')
+
+        # sample a single set of parameters from posterior distribution
+        if posterior:
+            p = posterior
+        elif posterior is None and self.posterior:
+            p = self.posterior
+        else:
+            raise ValueError('Need to provide posterior argument if no'
+                             ' posterior has been calculated by the `infere`'
+                             ' method.')
+        x_o = torch.tensor(self.x_o, dtype=torch.float32)
+        params = p.sample((n_samples, ), x=x_o).mean(axis=0)
+
+        # set output variable that is monitored
+        if output_var is None:
+            output_var = self.output_var
+
+        # set up initial values
+        if param_init is None:
+            param_init = self.param_init
+        else:
+            param_init = dict(self.param_init)
+            self.param_init.update(param_init)
+
+        # create a dictionary with repeated sampled prior
+        d_param = get_param_dict(params, self.param_names, self.n_traces)
+
+        # set up and run the simulator
+        network_name = 'generate_traces'
+        simulator = self.setup_simulator('generate_traces',
+                                         self.n_traces,
+                                         output_var=output_var,
+                                         param_init=param_init,
+                                         level=level+1)
+        simulator.run(self.sim_time, d_param, self.param_names, iteration=0,
+                      name=network_name)
+
+        # create dictionary of traces for multiple observed output variables
+        if len(output_var) > 1:
+            for ov in output_var:
+                if ov == 'spikes':
+                    trace = get_spikes(simulator.spikemonitor, 1,
+                                       self.n_traces)[0]
+                    traces = {ov: trace}
+                else:
+                    try:
+                        trace = getattr(simulator.statemonitor, ov)[:]
+                        traces = {ov: trace}
+                    except KeyError:
+                        logger.warn('No state monitor object found.'
+                                    ' Call again with specified `output_var`.')
+        else:
+            traces = getattr(simulator.statemonitor, output_var[0])[:]
+        return traces
diff --git a/brian2modelfitting/utils.py b/brian2modelfitting/utils.py
index 984e0a6b..534c6960 100644
--- a/brian2modelfitting/utils.py
+++ b/brian2modelfitting/utils.py
@@ -2,7 +2,7 @@
 
 from brian2 import have_same_dimensions
 from brian2.units.fundamentalunits import Quantity
-from tqdm.autonotebook import tqdm
+from tqdm.auto import tqdm
 
 
 def _format_quantity(v, precision=3):
@@ -12,7 +12,8 @@ def _format_quantity(v, precision=3):
         return f'{v:.{precision}g}'
 
 
-def callback_text(params, errors, best_params, best_error, index, additional_info):
+def callback_text(params, errors, best_params, best_error, index,
+                  additional_info):
     """Default callback print-out for Fitters"""
     params = []
     for p, v in sorted(best_params.items()):
@@ -25,14 +26,11 @@ def callback_text(params, errors, best_params, best_error, index, additional_inf
         for error, normed_error, varname in zip(additional_info['objective_errors'],
                                                 additional_info['objective_errors_normalized'],
                                                 additional_info['output_var']):
-
             if not have_same_dimensions(error, normed_error) or error != normed_error:
                 raw_error_str = f', unnormalized error: {_format_quantity(error)}'
             else:
                 raw_error_str = ''
-
             errors.append(f'{_format_quantity(normed_error)} ({varname}{raw_error_str})')
-
         error_sum = ' + '.join(errors)
         print(f"{round}Best parameters {param_str}\n"
               f"{' '*len(round)}Best error: {best_error_str} = {error_sum}")
@@ -53,8 +51,8 @@ def callback_text(params, errors, best_params, best_error, index, additional_inf
             print(f"{' ' * len(round)}Best error: {best_error_str} ({additional_info['output_var'][0]})")
 
 
-
-def callback_none(params, errors, best_params, best_error, index, additional_info):
+def callback_none(params, errors, best_params, best_error, index,
+                  additional_info):
     """Non-verbose callback"""
     pass
 
@@ -83,14 +81,12 @@ def callback_setup(set_type, n_rounds):
     elif type(set_type) is FunctionType:
         callback = set_type
     else:
-        raise TypeError("callback has to be a str ('text' or 'progressbar'), "
-                        "callable or None")
-
+        raise TypeError('callback has to be a str (`text` or `progressbar`), '
+                        'allable or None')
     return callback
 
 
 def make_dic(names, values):
     """Create dictionary based on list of strings and 2D array"""
     result_dict = {name: value for name, value in zip(names, values)}
-
     return result_dict
diff --git a/docs_sphinx/_static/hh_sbi_4params_4features_pairplot.png b/docs_sphinx/_static/hh_sbi_4params_4features_pairplot.png
new file mode 100644
index 00000000..af5ab24c
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_4features_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_9features_pairplot.png b/docs_sphinx/_static/hh_sbi_4params_9features_pairplot.png
new file mode 100644
index 00000000..e1425b43
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_9features_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_default_automatic_ext_pairplot.png b/docs_sphinx/_static/hh_sbi_4params_default_automatic_ext_pairplot.png
new file mode 100644
index 00000000..debe0121
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_default_automatic_ext_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_inf_traces.png b/docs_sphinx/_static/hh_sbi_4params_inf_traces.png
new file mode 100644
index 00000000..5859865b
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_inf_traces.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_multiround_4features_pairplot.png b/docs_sphinx/_static/hh_sbi_4params_multiround_4features_pairplot.png
new file mode 100644
index 00000000..3d8393bd
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_multiround_4features_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_spike_train.png b/docs_sphinx/_static/hh_sbi_4params_spike_train.png
new file mode 100644
index 00000000..06cf9897
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_spike_train.png differ
diff --git a/docs_sphinx/_static/hh_sbi_4params_yulenet_pairplot.png b/docs_sphinx/_static/hh_sbi_4params_yulenet_pairplot.png
new file mode 100644
index 00000000..dc4f6904
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_4params_yulenet_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_recording_traces.png b/docs_sphinx/_static/hh_sbi_recording_traces.png
new file mode 100644
index 00000000..bd0f15b7
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_recording_traces.png differ
diff --git a/docs_sphinx/_static/hh_sbi_toy_inf_traces.png b/docs_sphinx/_static/hh_sbi_toy_inf_traces.png
new file mode 100644
index 00000000..21e46760
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_toy_inf_traces.png differ
diff --git a/docs_sphinx/_static/hh_sbi_toy_pairplot.png b/docs_sphinx/_static/hh_sbi_toy_pairplot.png
new file mode 100644
index 00000000..d989c668
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_toy_pairplot.png differ
diff --git a/docs_sphinx/_static/hh_sbi_toy_pairplot_zoom_in.png b/docs_sphinx/_static/hh_sbi_toy_pairplot_zoom_in.png
new file mode 100644
index 00000000..97c22068
Binary files /dev/null and b/docs_sphinx/_static/hh_sbi_toy_pairplot_zoom_in.png differ
diff --git a/docs_sphinx/api/brian2modelfitting.inferencer.rst b/docs_sphinx/api/brian2modelfitting.inferencer.rst
new file mode 100644
index 00000000..30cf9c6a
--- /dev/null
+++ b/docs_sphinx/api/brian2modelfitting.inferencer.rst
@@ -0,0 +1,7 @@
+brian2modelfitting.inferencer module
+======================================
+
+.. automodule:: brian2modelfitting.inferencer
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs_sphinx/api/brian2modelfitting.rst b/docs_sphinx/api/brian2modelfitting.rst
index d09c81a9..4a10106c 100644
--- a/docs_sphinx/api/brian2modelfitting.rst
+++ b/docs_sphinx/api/brian2modelfitting.rst
@@ -8,9 +8,10 @@ Subpackages and -modules
 
 .. toctree::
 
-   brian2modelfitting.tests
-   brian2modelfitting.metric
    brian2modelfitting.fitter
+   brian2modelfitting.inferencer
+   brian2modelfitting.metric
    brian2modelfitting.optimizer
    brian2modelfitting.simulator
+   brian2modelfitting.tests
    brian2modelfitting.utils
diff --git a/docs_sphinx/conf.py b/docs_sphinx/conf.py
index 6cb48469..fc49303a 100644
--- a/docs_sphinx/conf.py
+++ b/docs_sphinx/conf.py
@@ -10,6 +10,7 @@
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
+
 import os
 import sys
 sys.path.insert(0, os.path.abspath('..'))
@@ -17,11 +18,12 @@
 needs_sphinx = '1.8'
 
 
-brian2modelfitting_dir = os.path.abspath(os.path.join(os.path.dirname(__file__),
-                                               '..', 'brian2modelfitting'))
+brian2modelfitting_dir = os.path.abspath(os.path.join(
+    os.path.dirname(__file__), '..', 'brian2modelfitting')
+    )
 # -- Project information -----------------------------------------------------
 
-project = 'Brian2modelfitting'
+project = 'brian2modelfitting'
 copyright = '2019, brian-team'
 author = 'brian-team'
 
@@ -66,7 +68,8 @@
 autodoc_default_options = {'inherited-members': True}
 
 # -- Options for HTML output -------------------------------------------------
-# on_rtd is whether we are on readthedocs.org, this line of code grabbed from docs.readthedocs.org
+# on_rtd is whether we are on readthedocs.org, this line of code grabbed from
+# docs.readthedocs.org
 on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
 
 if not on_rtd:
@@ -92,10 +95,12 @@
 
 
 # -- Extension configuration -------------------------------------------------
-intersphinx_mapping = {'python': ('https://docs.python.org/', None),
-                       'brian2': ('https://brian2.readthedocs.org/en/stable/', None),
-                       'matplotlib': ('http://matplotlib.org/', None),
-                       'numpy': ('http://docs.scipy.org/doc/numpy/', None),
-                       'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None),
-                       'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
-                       'lmfit': ('https://lmfit.github.io/lmfit-py/', None)}
+intersphinx_mapping = {
+    'python': ('https://docs.python.org/', None),
+    'brian2': ('https://brian2.readthedocs.org/en/stable/', None),
+    'matplotlib': ('http://matplotlib.org/', None),
+    'numpy': ('http://docs.scipy.org/doc/numpy/', None),
+    'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None),
+    'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
+    'lmfit': ('https://lmfit.github.io/lmfit-py/', None),
+    }
diff --git a/docs_sphinx/examples/base.rst b/docs_sphinx/examples/base.rst
index a6f4c7df..9c034184 100644
--- a/docs_sphinx/examples/base.rst
+++ b/docs_sphinx/examples/base.rst
@@ -1,7 +1,8 @@
 Simple Examples
 ===============
 
-Following pieces of code show an example of Fitter class calls with possible inputs.
+Following pieces of code show an example of two ``Fitter`` class calls and an
+``Inferencer`` class call with possible inputs.
 
 
 TraceFitter
@@ -27,8 +28,6 @@ TraceFitter
                               g_na=[1*mS*cm**-2 * area, 2000*mS*cm**-2 * area],
                               g_kd=[1*mS*cm**-2 * area, 1000*mS*cm**-2 * area])
 
-
-
 SpikeFitter
 -----------
 
@@ -51,3 +50,31 @@ SpikeFitter
                    metric=metric,
                    gL=[20*nS, 40*nS],
                    C = [0.5*nF, 1.5*nF])
+
+Inferencer
+----------
+
+.. code:: python
+
+  v_features = [
+      lambda x: max(x[(t > t_start) & (t < t_end)]),  # AP max
+      lambda x: mean(x[(t > t_start) & (t < t_end)]),  # AP mean
+      lambda x: std(x[(t > t_start) & (t < t_end)]),  # AP std
+      lambda x: mean(x[(t > .25 * t_start) & (t < .75 * t_start)]),  # resting
+  ]
+  s_features = [lambda x: x.size]  # number of spikes in a spike train
+
+  inferencer = Inferencer(model=eqs, dt=0.1*ms,
+                          input={'I': inp_traces},
+                          output={'v': out_traces, 'spike': spike_times}, 
+                          features={'v': v_features, 'spikes': s_features},
+                          method='exponential_euler',
+                          threshold='m > 0.5',
+                          refractory='m > 0.5',
+                          param_init={'v': 'VT'})
+
+  posterior = inferencer.infer(n_samples=1_000,
+                               n_rounds=2,
+                               inference_method='SNPE',
+                               gL=[20*nS, 40*nS],
+                               C = [0.5*nF, 1.5*nF])
diff --git a/docs_sphinx/examples/hh_multirun.rst b/docs_sphinx/examples/hh_multirun.rst
index e155a3b4..d1ed8955 100644
--- a/docs_sphinx/examples/hh_multirun.rst
+++ b/docs_sphinx/examples/hh_multirun.rst
@@ -1,5 +1,5 @@
-Multirun of Hodgkin-Huxley
-==========================
+Multirun fitting of Hodgkin-Huxley
+==================================
 
 Here you can download the data:
 :download:`input_traces <../../examples/input_traces_hh.csv>`
diff --git a/docs_sphinx/examples/hh_sbi_flex.rst b/docs_sphinx/examples/hh_sbi_flex.rst
new file mode 100644
index 00000000..37b58eaa
--- /dev/null
+++ b/docs_sphinx/examples/hh_sbi_flex.rst
@@ -0,0 +1,188 @@
+Inference on Hodgin-Huxley model: flexible interface
+====================================================
+
+You can also download and run this example by clicking here:
+:download:`hh_sbi_simple_interface.py <../../examples/hh_sbi_flexible_interface.py>`
+
+Here you can download the data:
+:download:`input traces <../../examples/input_traces_hh.csv>`
+:download:`output traces <../../examples/output_traces_hh.csv>`
+
+.. code:: python
+
+  from brian2 import *
+  from brian2modelfitting import *
+  import pandas as pd
+
+
+To load the data, use the following:
+
+.. code:: python
+  
+  df_inp_traces = pd.read_csv('input_traces_hh.csv')
+  df_out_traces = pd.read_csv('output_traces_hh.csv')
+  inp_traces = df_inp_traces.to_numpy()
+  inp_traces = inp_traces[[0, 1, 3], 1:]
+  out_traces = df_out_traces.to_numpy()
+  out_traces = out_traces[[0, 1, 3], 1:]
+
+The model used for this example is the Hodgkin-Huxley neuron model.
+The parameters of the model are defined as follows:
+
+.. code:: python
+
+  area = 20_000*um**2
+  El = -65*mV
+  EK = -90*mV
+  ENa = 50*mV
+  VT = -63*mV
+  dt = 0.01*ms
+  eqs = '''
+      dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
+      dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
+              (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
+              (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
+      dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
+              (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
+      dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
+
+      # unknown parameters
+      g_na : siemens (constant)
+      g_kd : siemens (constant)
+      gl   : siemens (constant)
+      Cm   : farad (constant)
+      '''
+
+Now, let's define the time domain and start with the inferencer procedure
+manually:
+
+.. code:: python
+    
+  t = arange(0, out_traces.shape[1]*dt/ms, dt/ms)
+  t_start, t_end = t[where(inp_traces[0, :] != 0)[0][[0, -1]]]
+
+  # Start with the regular instatiation of the class
+  inferencer = Inferencer(dt=dt, model=eqs,
+                          input={'I': inp_traces*amp},
+                          output={'v': out_traces*mV},
+                          features={'v': [lambda x: max(x),
+                                          lambda x: mean(x[(t > t_start) & (t < t_end)]),
+                                          lambda x: std(x[(t > t_start) & (t < t_end)])]},
+                          method='exponential_euler',
+                          threshold='m > 0.5',
+                          refractory='m > 0.5',
+                          param_init={'v': 'VT'})
+
+The prior should be initialized by defining the upper and lower bounds for
+each unknown parameter:
+
+.. code:: python
+  
+  prior = inferencer.init_prior(gl=[1e-09*siemens, 1e-07*siemens],
+                                g_na=[2e-06*siemens, 2e-04*siemens],
+                                g_kd=[6e-07*siemens, 6e-05*siemens],
+                                Cm=[0.1*uF*cm**-2*area, 2*uF*cm**-2*area])
+
+If the input and output data for the training of the neural density estimator
+already exists, we can load it as follows:
+
+.. code:: python
+
+  path_to_data = ...
+  theta, x = inferencer.load_summary_statistics(path_to_data)
+
+Otherwise, we have to generate training data and summary statistics from a
+given list of features:
+
+.. code:: python
+
+  theta = inferencer.generate_training_data(n_samples=10_000,
+                                            prior=prior)
+  x = inferencer.extract_summary_statistics(theta)
+
+And the data can be saved for the later use:
+
+.. code:: python
+ 
+  inferencer.save_summary_statistics(path_to_data, theta, x)
+
+Finally, let's get our hands dirty and let's perform a single step of
+inference:
+
+.. code:: python
+
+  # amortized inference
+  inference = inferencer.init_inference(inference_method='SNPE',
+                                        density_estimator_model='mdn',
+                                        prior=prior)
+  # first round of inference where no observation data is set to posterior
+  posterior_amortized = inferencer.infer_step(proposal=prior,
+                                              inference=inference,
+                                              theta=theta, x=x)
+
+After the posterior has been built, it can be stored as follows:                                              
+
+.. code:: python
+  
+  # storing the trained posterior without a default observation
+  path_to_posterior = ...
+  inferencer.save_posterior(path_to_posterior)
+
+Now, as in the simple interface example, sampling can be performed via
+``sample`` method where it is enough to define a number of parameters to
+be drawn from the posterior:
+
+.. code:: python
+
+  inferencer.sample((10_000, ))
+
+
+Creating the pairwise relationship visualizations using the approximated
+posterior distribution
+
+.. code:: python
+
+  # define the label for each parameter
+  labels = {'gl': r'$\overline{g}_\mathrm{l}$',
+            'g_na': r'$\overline{g}_\mathrm{Na}$',
+            'g_kd': r'$\overline{g}_\mathrm{K}$',
+            'Cm': r'$\overline{C}_{m}$'}
+  inferencer.pairplot(labels=labels)
+
+
+It is possible to continue with the focused inference (to draw parameters
+from the posterior and to perform the training of a neural network to
+estimate the posterior distribution by focusing on a particular observation)
+by using a standard approach through ``infer`` method:
+
+.. code:: python
+
+  posterior_focused = inferencer.infer()
+
+For every future call of ``inferencer``, the last trained posterior will be
+used by default, e.g., when generating traces by using a single sample of
+parameters from a now non-amortized approximated posterior distribution:
+
+.. code:: python
+
+  inf_traces = inferencer.generate_traces()
+  nrows = 2
+  ncols = out_traces.shape[0]
+  fig, axs = subplots(nrows, ncols, sharex=True,
+                      gridspec_kw={'height_ratios': [3, 1]},
+                      figsize=(ncols * 3, 3))
+  for idx in range(ncols):
+      axs[0, idx].plot(t, out_traces[idx, :].T, 'C3-', lw=3, label='recordings')
+      axs[0, idx].plot(t, inf_traces[idx, :].T/mV, 'k--', lw=2,
+                       label='sampled traces')
+      axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimuli')
+     axs[1, idx].set_xlabel('$t$, ms')
+      if idx == 0:
+          axs[0, idx].set_ylabel('$V$, mV')
+          axs[1, idx].set_ylabel('$I$, nA')
+  handles, labels = [(h + l) for h, l
+                     in zip(axs[0, idx].get_legend_handles_labels(),
+                     axs[1, idx].get_legend_handles_labels())]
+  fig.legend(handles, labels)
+  tight_layout()
+  show()
diff --git a/docs_sphinx/examples/hh_sbi_simple.rst b/docs_sphinx/examples/hh_sbi_simple.rst
new file mode 100644
index 00000000..ef8d650d
--- /dev/null
+++ b/docs_sphinx/examples/hh_sbi_simple.rst
@@ -0,0 +1,161 @@
+Inference on Hodgin-Huxley model: simple interface
+==================================================
+
+You can also download and run a similar example available here:
+:download:`hh_sbi_simple_interface.py <../../examples/hh_sbi_simple_interface.py>`
+
+Here you can download the data:
+:download:`input traces <../../examples/input_traces_hh.csv>`
+:download:`output traces <../../examples/output_traces_hh.csv>`
+
+.. code:: python
+
+  from brian2 import *
+  from brian2modelfitting import *
+  import pandas as pd
+
+To load the data, use the following:
+
+.. code:: python
+
+  df_inp_traces = pd.read_csv('input_traces_hh.csv')
+  df_out_traces = pd.read_csv('output_traces_hh.csv')
+  inp_traces = df_inp_traces.to_numpy()
+  inp_traces = inp_traces[[0, 1], 1:]
+  out_traces = df_out_traces.to_numpy()
+  out_traces = out_traces[[0, 1], 1:]
+
+Then we have to define the model and its parameters:
+
+.. code:: python
+
+  area = 20_000*um**2
+  El = -65*mV
+  EK = -90*mV
+  ENa = 50*mV
+  VT = -63*mV
+  dt = 0.01*ms
+  eqs = '''
+      dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
+      dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
+              (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
+              (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
+      dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
+              (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
+      dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
+
+      # free parameters
+      g_na : siemens (constant)
+      g_kd : siemens (constant)
+      gl   : siemens (constant)
+      Cm   : farad (constant)
+  '''
+
+Let's also specify time domain for more convenient plotting afterwards:
+
+.. code:: python
+
+  t = arange(0, out_traces.shape[1]*dt/ms, dt/ms)
+  stim_start, stim_end = t[where(inp_traces[0, :] != 0)[0][[0, -1]]]
+
+Now, we have to define features in order to create a summary statistics representation of the output data traces:
+
+.. code:: python
+  
+    list_of_features = [
+        lambda x: max(x[(t > stim_start) & (t < stim_end)]),  # max active potential
+        lambda x: mean(x[(t > stim_start) & (t < stim_end)]),  # mean active potential
+        lambda x: std(x[(t > stim_start) & (t < stim_end)]),  # std active potential
+        lambda x: mean(x[(t > .25 * stim_start) & (t < .75 * stim_start)]),  # resting
+    ]
+
+We have to instantiate the object by using the class ``Inferencer`` in which the data and the list of features should be passed:
+
+.. code:: python
+
+  inferencer = Inferencer(dt=dt, model=eqs,
+                          input={'I': inp_traces*amp},
+                          output={'v': out_traces*mV},
+                          features={'v': list_of_features},
+                          method='exponential_euler',
+                          threshold='m > 0.5',
+                          refractory='m > 0.5',
+                          param_init={'v': 'VT'})
+
+
+
+Be sure that the names of parameters passed to the ``infer`` method correspond to the names of unknown parameters defined as constatns in the model equations.
+
+.. code:: python
+
+  posterior = inferencer.infer(n_samples=5_000,
+                               n_rounds=3,
+                               inference_method='SNPE',
+                               density_estimator_model='mdn',
+                               gl=[1e-09*siemens, 1e-07*siemens],
+                               g_na=[2e-06*siemens, 2e-04*siemens],
+                               g_kd=[6e-07*siemens, 6e-05*siemens],
+                               Cm=[0.1*uF*cm**-2*area, 2*uF*cm**-2*area])
+
+After the training of the neural density estimator stored accessible through ``posterior`` is done, we can draw samples from the approximated posterior distribution as follows:
+
+.. code:: python
+
+  samples = inferencer.sample((5_000, ))
+
+In order to analyze the sampled data further, we can use the embedded ``pairplot`` method which visualizes the pairwise relationship between each two parameters:
+
+.. code:: python
+
+  limits = {'gl': [1e-9*siemens, 1e-07*siemens],
+            'g_na': [2e-06*siemens, 2e-04*siemens],
+            'g_kd': [6e-07*siemens, 6e-05*siemens],
+            'Cm': [0.1*uF*cm**-2*area, 2*uF*cm**-2*area]}
+  labels = {'gl': r'$\overline{g}_{l}$',
+            'g_na': r'$\overline{g}_{Na}$',
+            'g_kd': r'$\overline{g}_{K}$',
+            'Cm': r'$C_{m}$'}
+  inferencer.pairplot(limits=limits,
+                      labels=labels,
+                      ticks=limits,
+                      figsize=(6, 6))
+  condition = inferencer.sample((1, ))
+  inferencer.conditional_pairplot(condition=condition,
+                                  limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  figsize=(6, 6))
+
+To obtain a simulated trace from a single sample of parameters drawn from posterior distribution, use the following code:
+
+.. code:: python
+
+  inf_traces = inferencer.generate_traces(output_var='v')
+
+Let us now visualize the recordings and simulated traces:
+
+.. code:: python
+
+  inf_traces = inferencer.generate_traces(output_var='v')
+
+  nrows = 2
+  ncols = out_traces.shape[0]
+  fig, axs = subplots(nrows, ncols, sharex=True,
+                      gridspec_kw={'height_ratios': [3, 1]}, figsize=(9, 3))
+  for idx in range(ncols):
+      spike_idx = in1d(t, spike_times[idx]).nonzero()[0]
+      spike_v = (out_traces[idx, :].min(), out_traces[idx, :].max())
+      axs[0, idx].plot(t, out_traces[idx, :].T, 'C3-', lw=3, label='recordings')
+      axs[0, idx].plot(t, inf_traces[idx, :].T/mV, 'k--', lw=2,
+                       label='sampled traces')
+      axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimuli')
+      axs[1, idx].set_xlabel('$t$, ms')
+      if idx == 0:
+          axs[0, idx].set_ylabel('$V$, mV')
+          axs[1, idx].set_ylabel('$I$, nA')
+  handles, labels = [(h + l) for h, l
+                     in zip(axs[0, idx].get_legend_handles_labels(),
+                     axs[1, idx].get_legend_handles_labels())]
+  fig.legend(handles, labels)
+  tight_layout()
+  show()
diff --git a/docs_sphinx/examples/index.rst b/docs_sphinx/examples/index.rst
index 25f08309..957a08ad 100644
--- a/docs_sphinx/examples/index.rst
+++ b/docs_sphinx/examples/index.rst
@@ -6,3 +6,5 @@ Examples
 
    base
    hh_multirun
+   hh_sbi_simple
+   hh_sbi_flex
diff --git a/docs_sphinx/features/index.rst b/docs_sphinx/features/index.rst
index 9fc32457..806ec6ba 100644
--- a/docs_sphinx/features/index.rst
+++ b/docs_sphinx/features/index.rst
@@ -1,8 +1,10 @@
 Advanced Features
 =================
 
-This part of documentation list other features provided alongside or inside `~brian2modelfitting.fitter.Fitter`
-objects, to help the user with easier and more flexible applications.
+This part of documentation lists other features provided alongside or inside
+`~brian2modelfitting.fitter.Fitter` and `~brian2modelfitting.inferencer.Inferencer`
+objects to allow users easier and a more flexible development when working on
+their own problems.
 
 .. contents::
     :local:
@@ -11,108 +13,199 @@ objects, to help the user with easier and more flexible applications.
 Parameters initialization
 -------------------------
 
-Whilst running `~brian2modelfitting.fitter.Fitter` user can specify values with which model evaluation
-of differential equations start.
+Whilst running `~brian2modelfitting.fitter.Fitter` or `~brian2modelfitting.inferencer.Inferencer`,
+the user is able to pass the values of the parameters and variables that will
+be used as initial conditions when solving the differential equations defined
+in the neuron model.
 
-The fitting functions accept additional dictionary input to address that. To do so,
-dictionary argument has to be added to `~brian2modelfitting.fitter.Fitter` initialization:
+Initial conditions should be passed by using an additional dictionary to the
+constructor:
 
 .. code:: python
 
-   param_init = {'v': -30*mV}
+  init_conds = {'v': -30*mV}
 
 
 .. code:: python
 
-   fitter = TraceFitter(..., param_init = {'v': -30*mV})
+  fitter = TraceFitter(..., param_init = init_conds)
 
+or
+
+.. code:: python
+
+  inferencer = Inferencer(..., param_init=init_conds)
 
 Restart
 -------
-By default any `~brian2modelfitting.fitter.Fitter` works in continuous optimization mode between run, where all of the
-parameters drawn are being evaluated.
+By default any `~brian2modelfitting.fitter.Fitter` object works in continuous
+optimization mode between run, where all of the parameters drawn are being
+evaluated.
 
-Through changing the input flag in `~brian2modelfitting.fitter.Fitter.fit()`: ``restart`` to ``True``, user can reset the optimizer and
-start the optimization from scratch.
+By setting the ``restart`` argument in `~brian2modelfitting.fitter.Fitter.fit()`
+to ``True``, the user can restart the optimizer and the optimization will
+start from scratch.
 
-Used by Fitter optimizer and metric can only be changed when the flat is `True`.
+Used by Fitter optimizer and metric can only be changed when the flat is
+``True``.
 
+The previously outlined ``restart`` argument is used in the similar fashion
+in `~brian2modelfitting.inferencer.Inferencer.infer()` method. It is set to
+``False`` by default, and each following re-call of the method will result
+in the multi-round inference. If the user wants amortized inference without
+using any knowledge from the previous round of optimization instead, the
+``restart`` argument should be set to ``True``.
 
-Multiobjective normalization
+Multi-objective optimization
 ----------------------------
-It is possible to fit more than one output variable at the same time, by combining the errors for each variable. To do
-so, users can specify several output variables during the initialization::
+In the case of `~brian2modelfitting.fitter.Fitter` classes, it is possible to
+fit more than one output variable at the same time by combining the errors for
+each variable. To do so, the user can specify several output variables during
+the initialization as follows:
+
+.. code:: python
 
-    fitter = TraceFitter(..., output={'x': target_x,
-                                      'y': target_y},
-                         ...)
+  fitter = TraceFitter(...,
+                       output={'x': target_x,
+                               'y': target_y})
 
-If the fit uses a single metric, it applies to both variables. Note that this requires that the resulting error has the
-same units for both variables – for example, it would not be possible to use the same `.MSEMetric` on variables with
-different units, since the errors cannot be simply added up. As a more general solution, users can specify a metric for
-each variable and use their normalization arguments to make the units compatible (most commonly by turning both errors
-into dimensionless quantities). This normalization also defines the relative weights of both errors. For example, if the
-variable ``x`` has dimensions of volt and the variable ``y`` is dimensionless, the following metrics can be used to make
-an error of 10mV in ``x`` to be weighed as much as an error of 0.1 in ``y``::
+If the fitter function uses a single metric, it is applied to both variables.
+
+.. note::
+ 
+  This approach requires that the resulting error has the same units for all
+  variables, i.e., it would not be possible to use the same `.MSEMetric` on
+  variables with different units, since the errors cannot be simply added up.
+
+As a more general solution, the user can specify a metric for each variable
+and utilize their normalization arguments to make the units compatible (most
+commonly by turning both errors into dimensionless quantities). The
+normalization also defines the relative weights of all errors. For example, if
+the variable ``x`` has dimensions of mV and the variable ``y`` is
+dimensionless, the following metrics can be used to make an error of 10 mV in
+``x`` to be weighed as much as an error of 0.1 in ``y``
+
+.. code:: python
 
-    metrics = {'x': MSEMetric(normalization=10*mV)
-               'y': MSEMetric(normalization=0.1)}
+  metrics = {'x': MSEMetric(normalization=10*mV),
+             'y': MSEMetric(normalization=0.1)}
+
+This has to be passed as the ``metric`` argument of the `~brian2modelfitting.fitter.Fitter.fit`
+function.
+
+In the case of the `~brian2modelfitting.inferencer.Inferencer` class,
+switching from a single- to multi-objective optimization is seamless. The user
+has to provide multiple output variables during the initialization process the
+same way as for `~brian2modelfitting.fitter.Fitter` classes:
+
+.. code:: python
+
+  inferencer = Inferencer(...,
+                          output={'x': target_x,
+                                  'y': target_y})
+
+Later, during the inference process, the user has to define feautres for each
+output variable as follows:
+
+.. code:: python
 
-This dictionary then has to be provided as the ``metric`` argument of the `~.Fitter.fit` function.
+  posterior = inferencer.infer(...,
+                               features={'x': list_of_features_for_x,
+                                         'y': list_of_features_for_y})
 
+If the user prefers automatic feature extraction, the ``features`` argument
+should not be defined (it should stay set to None).
+
+.. warning::
+ 
+  If the user chooses to define a list of features for extracting the summary
+  features, it is important to keep in mind that the total number of features
+  will be increased as many times as there are output variables set for
+  multi-objective optimization.
 
 Callback function
 -----------------
 
-To visualize the progress of the optimization we provided few possibilities of feedback
-inside `~brian2modelfitting.fitter.Fitter`.
-
+To visualize the progress of the optimization we provided few possibilities of
+the feedback inside the `~brian2modelfitting.fitter.Fitter`.
 
 The 'callback' input provides few default options, updated in each round:
- - ``'text'`` (default) that prints out the parameters of the best fit and corresponding error
- - ``'progressbar'`` that uses ``tqdm.autonotebook`` to provide a progress bar
- - ``None`` for non-verbose option
+ - ``'text'`` (default) - prints out the parameters of the best fit and
+   corresponding error;
+ - ``'progressbar'`` - uses ``tqdm.autonotebook`` to provide a progress bar;
+ - ``None`` - non-verbose;
 
 as well as **customized feedback option**. User can provide
-a *callable* (i.e. function), that will provide an output or printout. If callback returns
-``True`` the fitting execution is interrupted.
+a *callable* (i.e., a function), that ensures either returning an output or
+printout. If callback returns ``True``, the fitting execution will be
+interrupted.
 
 User gets four arguments to customize over:
- - ``params`` - set of parameters from current round
- - ``errors`` - set of errors from current round
- - ``best_params`` - best parameters globally, from all rounds
- - ``best_error`` - best parameters globally, from all rounds
- - ``index`` - index of current round
+ - ``params`` - set of parameters from current round;
+ - ``errors`` - set of errors from current round;
+ - ``best_params`` - best parameters globally, from all rounds;
+ - ``best_error`` - best parameters globally, from all rounds;
+ - ``index`` - index of current round.
 
-An example function:
+An example callback function:
 
 .. code:: python
 
- def callback(params, errors, best_params, best_error, index):
-    print('index {} errors minimum: {}'.format(index, min(errors)))
-
-.. code:: python
+  def callback_fun(params, errors, best_params, best_error, index):
+      print('index {} errors minimum: {}'.format(index, min(errors)))
 
-   fitter = TraceFitter(...)
-   result, error  = fitter.fit(..., callback=...)
+  ...
 
+  fitter = TraceFitter(...)
+  result, error  = fitter.fit(..., callback=callback_fun)
 
+OnlineTraceFitter
+-----------------
+  
+`~brian2modelfitting.fitter.OnlineTraceFitter` was created to work with long
+traces or large-scale optimization problems. This `~brian2modelfitting.fitter.Fitter`
+class uses online mean square error as a metric.
+When the `~brian2modelfitting.fitter.OnlineTraceFitter.fit()` method is called
+there is no need of specifying a metric, which is by default set to None.
+The errors are instead calculated with `~brian2.groups.group.Group.run_regularly`
+for each simulation.
+  
+.. code:: python
+  
+  fitter = OnlineTraceFitter(model=model,
+                             input={'I': inp_traces},
+                             output={'v': out_traces},
+                             dt=0.1*ms,
+                             n_samples=5)
+  
+  result, error = fitter.fit(optimizer=optimizer,
+                             n_rounds=1,
+                             gl=[1e-8*siemens*cm**-2 * area, 1e-3*siemens*cm**-2 * area])
+  
+Reference the target values in the equations
+--------------------------------------------
+  
+A model can refer to the target output values within the equations. For
+example, if the membrane potential trace *v* (i.e. `output_var='v'`) is used
+for the optimization, equations can refer to the target trace as `v_target`.
+This allows adding a coupling term such as: `coupling*(v_target - v)` to
+the equation that corresponds to state variable `v`, pulling the trajectory
+towards the correct solution.
 
 Generate Traces
 ---------------
 
-With the same `~brian2modelfitting.fitter.Fitter` class user can also generate the traces with newly
-optimized parameters.
+`~brian2modelfitting.fitter.Fitter` and `~brian2modelfitting.inferencer.Inferencer`
+classes allow the user can to generate the traces with optimized parameters.
 
-To simulate and visualize the traces or spikes for the parameters of choice.
-For a quick access to best fitted set of parameters Fitter classes provided
-ready to use functions:
+For a quick access to best fitted set of parameters `~brian2modelfitting.fitter.Fitter`
+classes provide ready to use functions:
 
- - `~brian2modelfitting.fitter.TraceFitter.generate_traces` inside `~brian2modelfitting.fitter.TraceFitter`
- - `~brian2modelfitting.fitter.SpikeFitter.generate_spikes` inside `~brian2modelfitting.fitter.SpikeFitter`
+ - `~brian2modelfitting.fitter.TraceFitter.generate_traces` inside `~brian2modelfitting.fitter.TraceFitter`;
+ - `~brian2modelfitting.fitter.SpikeFitter.generate_spikes` inside `~brian2modelfitting.fitter.SpikeFitter`.
 
-Functions can be called after fitting in the following manner, without
-any input arguments:
+These functions can be called after the fitting procedure is finalized in the
+following manner, without any input arguments:
 
 .. code:: python
 
@@ -126,27 +219,46 @@ any input arguments:
     results, error = fitter.fit(...)
     spikes = fitter.generate_spikes()
 
+On the other hand, since the `~brian2modelfitting.inferencer.Inferencer` class
+is able to perform the inference of the unknown parameter distribution by
+utilizing output traces and spike trains simultaneously, ``generate_traces``
+is used for both.
 
-Custom generate
-~~~~~~~~~~~~~~~
+Once the approximated posterior distribution is built, the user is allowed to
+call ``generate_traces`` on `~brian2modelfitting.inferencer.Inferencer`
+object. If only one output variable is used for the optimization of the
+parameters, the user does not have to specifiy output variable in the 
+``generate_traces`` method through ``output_var`` argument. If, for example,
+the multi-objective optimization is performed by using both output traces and
+spike trains and the user is interested in only times of spike events,
+``output_var`` should be set to ``'spike'``. Otherwise, if the user specifies
+a list of names or the ``output_var`` is not specified, a dictionary with keys
+set to output variable names and with their respective values, will be
+returned instead.
+
+
+Customize the ``generate`` method for `~brian2modelfitting.fitter.Fitter`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-To create traces for other parameters, or generate traces after spike
-train fitting, user can call the - `~brian2modelfitting.fitter.Fitter.generate` call, that takes in following
-arguments:
+To create traces for other parameters, or generate traces after the spike
+train fitting, user can call the `~brian2modelfitting.fitter.Fitter.generate`
+method, which takes in the following arguments:
 
 .. code:: python
 
-  fitter.generate(params=None, output_var=None, param_init=None, level=0)
+  fitter.generate(params=..., output_var=..., param_init=..., level=0)
 
-Where ``params`` is a dictionary of parameters for which the traces we generate.
-``output_var`` provides an option to pick one or more variable for visualization. With
-``param_init``, user can define the initial values for differential equations.
-``level`` allows for specification of namespace level from which we get
-the constant parameters of the model.
+where ``params`` should be a dictionary of parameters for which we generate
+the traces; ``output_var`` provides an option to pick one or more variables
+for visualization; with ``param_init``, the user is able to define the initial
+values for differential equations in the model; and ``level`` allows for
+specification of the namespace level from which we are able to get the
+constant parameters of the model.
 
-If ``output_var`` is the name of a single variable name (or the special name ``'spikes'``), a single `~.Quantity`
-(for normal variables) or a list of spikes time arrays (for ``'spikes'``) will be returned. If a list of names is
-provided, then the result is a dictionary with all the results.
+If ``output_var`` is the name of a single variable name (or the special name
+``'spikes'``), a single `~.Quantity` (for variables) or a list of spikes time
+arrays (for ``'spikes'``) will be returned. If a list of names is provided,
+then the result is a dictionary with all the results.
 
 .. code:: python
 
@@ -155,43 +267,44 @@ provided, then the result is a dictionary with all the results.
     traces = fitter.generate(output_var=['v', 'h', 'n', 'm'])
     v_trace = traces['v']
     h_trace = traces['h']
-    ...
-
 
 Results
 -------
 
-Fitter class stores all of the parameters examined by the optimizer as well
-as the corresponding error. To retrieve them you can call the - `~brian2modelfitting.fitter.Fitter.results`.
+`~brian2modelfitting.fitter.Fitter` classes store all of the parameters
+used by the optimizer as well as the corresponding errors. To retrieve them
+you can call the `~brian2modelfitting.fitter.Fitter.results`.
 
 
 .. code:: python
 
-    fitter = TraceFitter(...)
-    ...
-    traces = fitter.generate_traces()
+  fitter = TraceFitter(...)
+  ...
+  traces = fitter.generate_traces()
 
 .. code:: python
 
-    fitter = SpikeFitter(...)
-    ...
-    results = fitter.results(format='dataframe')
+  fitter = SpikeFitter(...)
+  ...
+  results = fitter.results(format='dataframe')
 
 
 Results can be returned in one of the following formats:
 
- - ``'list'`` (default) returns a list of dictionaries with corresponding parameters (including units) and errors
- - ``'dict'`` returns a dictionary of arrays with corresponding parameters (including units) and errors
- - ``'dataframe'`` returns a `~pandas.DataFrame` (without units)
+ - ``'list'`` (default) - returns a list of dictionaries with corresponding
+   parameters (including units) and errors;
+ - ``'dict'`` - returns a dictionary of arrays with corresponding parameters
+   (including units) and errors;
+ - ``'dataframe'`` - returns a `~pandas.DataFrame` (without units).
 
-The use of units (only relevant for formats ``'list'`` and ``'dict'``) can be switched
-on or off with the ``use_units`` argument. If it is not specified, it will default to
-the value used during the initialization of the `Fitter` (which itself defaults to
-``True``).
+The use of units (only relevant for formats ``'list'`` and ``'dict'``) can be
+switched on or off with the ``use_units`` argument. If it is not specified, it
+will default to the value used during the initialization of the `Fitter`
+(which itself defaults to ``True``).
 
 Example output:
 ~~~~~~~~~~~~~~~
-``'list'``:
+- ``format='list'``:
 
 .. code:: python
 
@@ -200,7 +313,7 @@ Example output:
    ...]
 
 
-``'dict'``:
+- ``format='dict'``:
 
 .. code:: python
 
@@ -210,68 +323,122 @@ Example output:
    'errors': array([1.16788502, 0.5253008 , ...])}
 
 
-``'dataframe'``:
+- ``format='dataframe'``:
 
 .. code:: python
 
-   g_na            gl      g_kd    errors
-   0  0.000280  8.870238e-08  0.000047  0.521425
-   1  0.000192  1.121861e-07  0.000118  0.387140
-   ...
+     g_na      gl            g_kd      errors
+  0  0.000280  8.870238e-08  0.000047  0.521425
+  1  0.000192  1.121861e-07  0.000118  0.387140
+  ...
 
 
+Posterior distribution analysis
+-------------------------------
+
+Unlike `~brian2modelfitting.fitter.Fitter` classes, the `~brian2modelfitting.inferencer.Inferencer`
+class does not keep track of all parameter values. Rather, it stores all
+training data for neural density estimator which will later be used for
+building the posterior distribution of each unknown parameter. Thus, the `~brian2modelfitting.inferencer.Inferencer`
+does not returns best-fit values and corresponding errors, but the entire
+posterior distribution that can be used to draw samples from, compute
+descriptive statistics of parameters, analyize pairwise relationship between
+each to parameters, etc.
+
+There are three methods that enable the comprehensive analysis of the
+posterior:
+
+- `~brian2modelfitting.inferencer.Inferencer.pairplot` - returns axes of drawn
+  samples from the posterior in a 2-dimenstional grid with marginals and
+  pairwise marginals. Using this method, the user is able to inspect the
+  relationship for all combinations of distributions for each parameter;
+- `~brian2modelfitting.inferencer.Inferencer.conditional_pairplot` -
+  visualizes the conditional pairplot;
+- `~brian2modelfitting.inferencer.Inferencer.conditional_corrcoeff` - returns
+  the correlation matrix of a distribution conditioned with the user-specified
+  condition.
+
+To see this in action, go to our tutorial page and learn how to use each of
+these methods.
 
 Standalone mode
 ---------------
 
-Just like with regular Brian script, modelfitting computations can be performed in
-``Runtime`` mode (default) or ``Standalone`` mode.
-(https://brian2.readthedocs.io/en/stable/user/computation.html)
+Just like with regular Brian 2 scripts, all computations in the toolbox can be
+performed in ``Runtime`` mode (default) or ``Standalone`` mode. For details,
+please check the official Brian 2 documentation: https://brian2.readthedocs.io/en/stable/user/computation.html
 
-To enable this mode, add the following line after your Brian import, but before your simulation code:
+To enable the ``Standalone`` mode, and to allowthe source code generation to
+C++ code, add the following code right after Brian 2 is imported, but before
+the simulation code:
 
 .. code:: python
 
   set_device('cpp_standalone')
 
-
 Important notes:
 ~~~~~~~~~~~~~~~~
 
 .. warning::
-    In standlone mode one script can not be used to contain multiple - `~brian2modelfitting.fitter.Fitter`, use separate scripts!
+
+  In the ``Standalone`` mode, a single script should not contain multiple
+  `~brian2modelfitting.fitter.Fitter` or `~brian2modelfitting.inferencer.Inferencer`
+  classes. Please, use separate scripts.
 
 Note that the generation of traces or spikes via `~brian2modelfitting.fitter.Fitter.generate`
 will always use runtime mode, even when the fitting procedure uses standalone mode.
 
+Embedding network for automatic feature extraction
+--------------------------------------------------
 
-OnlineTraceFitter
------------------
-
-`~brian2modelfitting.fitter.OnlineTraceFitter` was created to work with long traces or big optimization.
-This `~brian2modelfitting.fitter.Fitter` uses online Mean Square Error as a metric.
-When `~brian2modelfitting.fitter.Fitter.fit()` is called there is no need of specifying a metric, that is by
-default set to None. Instead the errors are calculated with use of brian's `~brian2.groups.group.Group.run_regularly`,
-with each simulation.
+If the ``features`` argument of the `~brian2modelfitting.inferencer.Inferencer`
+class is not defined, automatic feature extraction from the given output
+traces will occur. By default, this is done by using the multi-layer
+perceptron that is trained in parallel with the neural density estimator of
+choice during the inference process. If the user wants to specify their own
+custom embedding network, it is possible to do so by creating a neural
+network by using ``PyTorch`` library and passing the instance of that neural
+network as an additional keyword argument as follows:
 
 .. code:: python
+  
+  import torch
+  from torch import nn
+  
+  ...
 
-  fitter = OnlineTraceFitter(model=model,
-                             input={'I': inp_traces},
-                             output={'v': out_traces},
-                             dt=0.1*ms,
-                             n_samples=5)
+  class CustomEmbeddingNet(nn.Module):
 
-  result, error = fitter.fit(optimizer=optimizer,
-                             n_rounds=1,
-                             gl=[1e-8*siemens*cm**-2 * area, 1e-3*siemens*cm**-2 * area],)
+      def __init__(self, in_features, out_features, ...):
+          ...
 
+      def forward(self, x):
+          ...
 
-Reference the target values in the equations
---------------------------------------------
+  
+  in_features = out_traces.shape[1]
+  out_features = ...
+  embedding_net = CustomEmbeddingNet(in_features, out_features, ...)
+
+  ...
+
+  inferencer = Inferencer(...)
+  inferencer.infer(...,
+                   inference_kwargs={'embedding_net': embedding_net})
+
+GPU usage for inference
+-----------------------
+
+It is possible to use the GPU for training the sdensity estimator. It is enough
+to specify the ``sbi_device`` to ``'gpu'`` or ``'cuda'``.  Otherwise, if not
+specified, or if set to ``'cpu'``, training will be done by using the CPU.
+
+.. note::
+
+  For default density estimators that are used either for SNPE, SNLE and SNRE,
+  there are no significant speed-ups expected if the training is translocated
+  to the GPU.
 
-A model can refer to the target output values within the equations. For example, if you
-are fitting a membrane potential trace *v* (i.e. `output_var='v'`), then the equations
-can refer to the target trace as `v_target`. This allows you for example to add a coupling
-term like `coupling*(v_target - v)` to the equation for `v`, pulling the trajectory towards the
-correct solution.
\ No newline at end of file
+It is, however, possible to achieve a significant speed-up if the custom
+embedding network relies on convolutions to extract feautres. Such operations
+are known to achieve improvement in compuation time multifold.
diff --git a/docs_sphinx/index.rst b/docs_sphinx/index.rst
index a4876c93..644b7d83 100644
--- a/docs_sphinx/index.rst
+++ b/docs_sphinx/index.rst
@@ -1,25 +1,41 @@
 brian2modelfitting
 ==================
 
-The package `.brian2modelfitting` is a tool for parameter fitting of neuron
-models in the `Brian 2 simulator <https://brian2.readthedocs.org>`_.
+The package `.brian2modelfitting` is a tool for parameter identification of
+neuron models defined in the `Brian 2 simulator <https://brian2.readthedocs.org>`_.
 
-Please report bugs at the `github issue tracker <https://github.com/brian-team/brian2modelfitting/issues>`_
-or at the `Brian discussion forum <https://groups.google.com/forum/#!forum/briansupport>`_. The latter is
-also a could place to discuss feature requests or potential contributions.
+Please report bugs at the `GitHub issue tracker <https://github.com/brian-team/brian2modelfitting/issues>`_
+or at the `Brian 2 discussion forum <https://groups.google.com/forum/#!forum/briansupport>`_.
+The latter is also a place to discuss feature requests or potential
+contributions.
 
 
 Model fitting
 -------------
 
-The ``brian2modelfitting`` toolbox offers allows the user to perform data driven
-optimization for custom neuronal models specified with
-`Brian 2 <https://brian2.readthedocs.org>`_.
+This toolbox allows the user to find the best fit of the unknown free
+parameters for recorded traces and spike trains.
+It also supports simulation-based inference, where instead of point-estimated
+parameter values, a full posterior distribution over the parameters is
+computed.
 
-The toolbox allows the user to find the best fit of the parameters for recorded traces and
-spike trains. Just like Brian itself, the Model Fitting Toolbox is designed to
-be easy to use and save time through automatic parallelization of the
-simulations using code generation.
+By default, the toolbox supports a range of global derivative-free
+optimization methods, that include popular methods for model fitting:
+differential evolution, particle swarm optimization and covariance matrix
+adaptation (provided by the ``nevergrad``, a gradient-free optimization
+platform) as well as Bayesian optimization for black box functions (provided
+by ``scikit-optimize``, a sequential model-based optimization library). On the
+other hand, simulation-based inference is the process of finding parameters of
+a simulator from observations by taking a Bayesian approach, in our case, via
+sequential neural posterior estimation, likelihood estimation or ratio
+estimation (provided by the ``sbi``), where neural densitiy estimator, a deep
+neural network allowing probabilistic association between the data and
+underlying parameter space, is trained. After the network is trained, the
+approximated posterior distribution is available.
+
+Just like Brian 2 simulator itself, the `.brian2modelfitting` toolbox is
+designed to be easy to use and to save time through automatic parallelization
+of the simulations using code generation.
 
 
 Contents
@@ -32,6 +48,7 @@ Contents
    introduction/index
    optimizer/index
    metric/index
+   inferencer/index
    features/index
    examples/index
 
@@ -48,4 +65,3 @@ Indices and tables
 
 * :ref:`genindex`
 * :ref:`modindex`
-* :ref:`search`
diff --git a/docs_sphinx/inferencer/index.rst b/docs_sphinx/inferencer/index.rst
new file mode 100644
index 00000000..f3fbc8d2
--- /dev/null
+++ b/docs_sphinx/inferencer/index.rst
@@ -0,0 +1,83 @@
+Inferencer
+==========
+
+Unlike more traditional inverse identification procedures that rely either on
+gradient or gradient-free methods, the `~brian2modelfitting.inferencer.Inferencer`
+class supports simulation-based inference that has been established as a
+powerful alternative approach.
+
+The simulation-based inference is data-driven procedure supported by the
+`sbi <https://www.mackelab.org/sbi/>`_, ``PyTorch``-based toolbox by Macke
+lab, [Tejero-Cantero2020]_.
+
+In general, this method yields twofold improvement over point-estimate fitting
+procedures:
+
+#. Simulation-based inference acts as if the actual statistical inference is
+   performed, even in cases of extremly complex models with untractable
+   likelihood function. Thus, instead of returning a single set of optimal
+   parameters, it results in the approximated posterior distribution over
+   unknown parameters. This is achieved by training a neural density estimator,
+   details of which will be explained in depth later in the documentation.
+#. Simulation-based inference uses prior system knowledge sparsely, using 
+   only the most important features to identify mechanistic models that are 
+   consistent with the recordings. This is achieved either by providing the 
+   predifend set of features, or by automatically extraciting summary features 
+   by using deep neural networks which is trained in parallel with neural 
+   density estimator. 
+
+The `~brian2modelfitting.inferencer.Inferencer` class, in its core, is a fancy
+wrapper around the ``sbi`` package,  where the focus is put on inferring the
+unknown parameters of the single-cell neuron models defined in Brian 2
+simulator.
+
+Neural density estimator
+------------------------
+
+There are three main estimation techniques supported in ``sbi`` that the user
+can take the full control over seamlesly by using the `~brian2modelfitting.inferencer.Inferencer`:
+
+#. sequential neural posterior estimation (SNPE)
+#. sequential neural likelihood estimation (SNLE)
+#. sequential neural ratio estimator (SNRE)
+
+Simulation-based inference workflow
+----------------------------------- 
+
+The inferencer procedure is defined via three main steps:
+
+#. step.
+   Prior over unknown parameters needs to be defined, where the simplest
+   choice would be uniform distribution given lower and upper bound
+   (currently, this is only prior distribution supported through
+   ``brian2modelfitting`` toolbox).
+   After that, simulated data are generated given a mechanistic model with
+   unknown parameters set as constants.
+   Instead of taking the full output of the model, the neural network takes
+   in summary data statistics of the output, e.g. instead of voltage trace as
+   the output from a neuron model, we would feed a neural network with
+   relevant electrophysiology features that outline the gist of the output
+   sufficiently well.
+#. step.
+   A neural network learns association between the summary data statistics
+   and unknown parameters (given the prior distribution over parameters).
+   The learning method is heavily dependent on the choice of the inference
+   technique.
+#. step.
+   The trained neural network is applied to the empirical data to infer
+   posterior distribution over unknown parameters.  Optionally, this process
+   can be repeated by using the trained posterior distribution over parameters
+   as the prior distribution proposal for a refined optimization.
+
+Implementation
+--------------
+
+Go to :doc:`the tutorial section <../introduction/tutorial_sbi>` for the
+in-depth implementation analysis.
+
+References
+----------
+
+.. [Tejero-Cantero2020] Tejero-Cantero, A., Boelts, J. et al. "sbi: A toolkit
+                        for simulation-based inference" Journal of Open Source
+                        Software 5(52):2505. 2020. doi: `10.21105/joss.02505 <https://doi.org/10.21105/joss.02505>`_
\ No newline at end of file
diff --git a/docs_sphinx/introduction/howitworks.rst b/docs_sphinx/introduction/howitworks.rst
index 39466c0f..c4c8bfee 100644
--- a/docs_sphinx/introduction/howitworks.rst
+++ b/docs_sphinx/introduction/howitworks.rst
@@ -1,10 +1,14 @@
+============
 How it works
 ============
 
+Fitting
+=======
+
 Model fitting script requires three components:
- - A **Fitter** of choice: object that will perform the optimization
- - A **metric**: to compare results and decide which one is the best
- - An **optimization** algorithm: to decide which parameter combinations to try
+ - a **fitter**: object that will perform the optimization
+ - a **metric**: objective function
+ - an **optimizer**: optimization algorithm
 
 All of which need to be initialized for fitting application.
 Each optimization works with a following scheme:
@@ -13,61 +17,63 @@ Each optimization works with a following scheme:
 
   opt = Optimizer()
   metric = Metric()
-
   fitter = Fitter(...)
   result, error = fitter.fit(metric=metric, optimizer=opt, ...)
 
 
-The proposed solution is developed using a modular approach, where both the optimization
-method and metric to be optimized can be easily swapped out by a custom implementation.
+The proposed solution is developed using a modular approach, where both the
+optimization method and the objective function can be easily swapped out by a
+user-defined custom implementation.
 
-`~.Fitter` objects require 'model' defined as an
+`~.Fitter` objects require a model defined as an
 `~brian2.equations.equations.Equations` object or as a string, that has
 parameters that will be optimized specified as constants in the following way:
 
 .. code:: python
 
   model = '''
-  ...
-  g_na : siemens (constant)
-  g_kd : siemens (constant)
-  gl   : siemens (constant)
-  '''
+      ...
+      g_na : siemens (constant)
+      g_kd : siemens (constant)
+      gl   : siemens (constant)
+      '''
 
 Initialization of Fitter requires:
-  - ``dt`` - time step
-  - ``input`` - a dictionary with the name of the input variable and a set of input traces (list or array)
-  - ``output`` - a dictionary with the name of the output variable(s) and a set of goal output (traces/spike trains) (list or array)
-  - ``n_samples`` - number of samples to draw in each round (limited by method)
-  - ``reset``, and ``threshold`` in case of spiking neurons (can take refractory as well)
+  - ``dt`` - the time step
+  - ``input`` - a dictionary with the name of the input variable and a set of
+  input traces (list or array)
+  - ``output`` - a dictionary with the name of the output variable(s) and a
+  set of goal output (traces/spike trains) (list or array)
+  - ``n_samples`` - a number of samples to draw in each round (limited by
+  method)
+  - ``reset`` and ``threshold`` in case of spiking neurons (can take
+  refractory as well)
 
 Additionally, upon call of `~brian2modelfitting.fitter.Fitter.fit()`,
 object requires:
 
- - ``n_rounds`` - number of rounds to optimize over
+ - ``n_rounds`` - a number of rounds to optimize over
  - parameters with ranges to be optimized over
 
-...as well as an ``optimizer`` and a ``metric``
-
 Each free parameter of the model that shall be fitted is defined by two values:
 
 .. code:: python
 
-  param_name = [min, max]
+  param_name = [lower_bound, upper_bound]
 
 Ready to use elements
 ---------------------
 
-Alongside three optimization classes:
+Optimization classes:
  - `~brian2modelfitting.fitter.TraceFitter`
  - `~brian2modelfitting.fitter.SpikeFitter`
  - `~brian2modelfitting.fitter.OnlineTraceFitter`
 
-We also provide ready optimizers:
+Optimization algorithms:
  - `~brian2modelfitting.optimizer.NevergradOptimizer`
  - `~brian2modelfitting.optimizer.SkoptOptimizer`
 
-and metrics:
+Metrics:
  - `~brian2modelfitting.metric.MSEMetric` (for `~brian2modelfitting.fitter.TraceFitter`)
  - `~brian2modelfitting.metric.GammaFactor` (for `~brian2modelfitting.fitter.SpikeFitter`)
 
@@ -85,11 +91,11 @@ Example of `~brian2modelfitting.fitter.TraceFitter` with all of the necessary ar
   result, error = fitter.fit(optimizer=optimizer,
                              metric=metric,
                              n_rounds=1,
-                             gl=[1e-8*siemens*cm**-2 * area, 1e-3*siemens*cm**-2 * area],)
+                             gl=[1e-8*siemens*cm**-2 * area, 1e-3*siemens*cm**-2 * area])
 
 Remarks
 -------
-- After performing first fitting, user can continue the optimization
+- After performing first fitting round, user can continue the optimization
   with another `~brian2modelfitting.fitter.Fitter.fit()` run.
 
 - Number of samples can not be changed between rounds or `~brian2modelfitting.fitter.Fitter.fit()`
@@ -101,3 +107,163 @@ Remarks
 
 - When using the `~brian2modelfitting.fitter.TraceFitter`, users can use a standard
   curve fitting algorithm for refinement by calling `~brian2modelfitting.fitter.TraceFitter.refine`.
+
+Simulation-based inference
+==========================
+
+The `~brian2modelfitting.inferencer.Inferencer` class has to be initialized
+within the script that will perform a simulation-based inference procedure.
+
+Initialization of `~brian2modelfitting.inferencer.Inferencer` requires:
+
+- ``dt`` - the time step in Brian 2 units.
+- ``model`` - single cell model equations, defined as either string or as 
+  ``brian2.Equation`` object.
+- ``input`` - a dictionary where key corresponds to the name of the input
+  variable as defined in ``model`` and value corresponds to an array of
+  input traces.
+- ``output`` - a dictionary where key corresponds to the name of the output
+  variable as defined in ``model`` and value corresponds to an array of
+  recorded traces and/or spike trains.
+
+.. code:: python
+
+  inferencer = Inferencer(dt=0.1*ms, model=eqs,
+                          input={'I': inp_traces*amp},
+                          output={'v': out_traces*mV})
+
+Optionally, arguments to be passed to the constructor are:
+
+- ``features`` - a dictionary of callables that take the voltage
+  trace and/or spike trains and output summary statistics. Keys correspond to
+  output variable names, while values are lists of callables. If features are
+  not provided, automatic feature extraction will be performed either by using
+  the default multi-layer perceptron or by using the user-provided embedding
+  network.
+- ``method`` - a string that defines an integration method.
+- ``threshold`` - optional string that defines the condition which produces
+  spikes. It should be a single line boolean expression.
+- ``reset`` - an optional (multi-line) string that that holds the code to
+  execute on reset.
+- ``refractory`` - can be either Boolean expression or string. Defines either
+  the length of the refractory period (e.g., ``2*ms``), a string expression
+  that evaluates to the length of the refractory period after each spike,
+  e.g., ``'(1 + rand())*ms'``, or a string expression evaluating to a boolean
+  value, given the condition under which the neuron stays refractory after a
+  spike, e.g., ``'v > -20*mV'``.
+- ``param_init`` - a dictionary of state variables to be initialized with
+  respective values, i.e., initial conditions.
+
+.. code:: python
+
+  inferencer = Inferencer(dt=dt, model=eqs_inf,
+                          input={'I': inp_trace*amp},
+                          output={'v': out_traces*mV},
+                          features={'v': voltage_feature_list},
+                          method='exponential_euler',
+                          threshold='v > -50*mV',
+                          reset='v = -70*mV',
+                          param_init={'v': -70*mV})
+
+Inference
+---------
+
+After the `~brian2modelfitting.inferencer.Inferencer` class is
+instantiated, the simplest and the most convenient way to start with the
+inferencer procedure is by calling `~brian2modelfitting.inferencer.Inferencer.infer`
+method on `~brian2modelfitting.inferencer.Inferencer` object.
+
+In the nutshell, infer method returns the trained neural posterior object,
+which may or may not be used by the user, but it has to exist. There are two
+possible approaches:
+
+- amortized inference
+- multi-round inference
+
+If the number of inference rounds is 1, then amortized inference will be
+performed. Otherwise if the number of inference rounds is 2 or above, the
+focused multi-round inference will be performed. Multi-round inference,
+unlike the amortized one, is focused on a particular observation, where in
+each new round of inference, samples are drawn from the posterior distribution
+conditioned exactly by this observation. This process can be repeated
+aribtrarily many times to get increasingly better approximations of the the
+posterior distribution.
+
+The infer method requires:
+
+- ``n_samples`` - the number of samples from which the neural posterior will
+  be learnt.
+  
+or:
+
+- ``theta`` -  sampled prior.
+- and ``x`` - summary statistics.
+
+along with the:
+
+- ``params`` - a dictionary of bounds for each free parameter defined in the
+  ``model``. Keys should correspond to names of parameters as defined in the
+  model equations, values are lists with lower and upper bounds with
+  quantities of respective parameter. 
+
+The simplest way to start the inference process is by calling:
+
+.. code:: python
+
+  posterior = inferencer.infer(n_samples=1000,
+                               gl=[10*nS, 100*nS],
+                               C=[0.1*nF, 10*nF])
+
+Optionally, user can defined the following arguments:
+
+- ``n_rounds`` - if it is set to 1, amortized inference will be performed.
+  Otherwise, if ``n_rounds`` is integer larger than 1, multi-round inference
+  will be performed. This is only valid if the posterior has not yet been
+  defined. Otherwise, if this method is called after the posterior has already
+  been built, multi-round inference is performed, e.g. repeated calling of
+  ``~brian2modelfitting.inferencer.Inferencer.infer`` method or manually
+  building the posterior by approaching the inference with flexible inference.
+- ``inference_method`` - either SNPE, SNLE or SNRE.
+- ``density_estimator_model`` - string that defines the type of density
+  estimator to be created. Either ``mdn``, ``made``, ``maf``, ``nsf`` for SNPE
+  and SNLE, or ``linear``, ``mlp``, ``resnet`` for SNRE.
+- ``inference_kwargs`` - a dictionary that holds additional keyword arguments
+  for the `~brian2modelfitting.inferencer.Inferencer.init_inference`.
+- ``train_kwargs`` - a dictionary that holds additional keyword arguments for
+  `~brian2modelfitting.inferencer.Inferencer.train`.
+- ``posterior_kwargs`` - a dictionary that holds additional keyword arguments
+  for `~brian2modelfitting.inferencer.Inferencer.build_posterior`.
+- ``restart`` - when the method is called for a second time, set to True if
+  amortized inference should be performed. If False, multi-round inference
+  with the existing posterior will be performed.
+- ``sbi_device`` a string that defines the device on which the ``sbi`` and
+  subseqently the ``torch`` will operate. By default this is set to ``cpu``
+  and it is advisable to remain so for most cases. In cases where the user
+  provides custom embedding network through ``inference_kwargs`` argument,
+  which will be trained more efficiently by using GPU, device should be set
+  accordingly to ``gpu``.
+
+A bit more comprehensive specification of the infer call is showcased below:
+
+.. code:: python
+
+  posterior = inferencer.infer(n_samples=5_000,
+                               n_rounds=3,
+                               inference_method='SNPE',
+                               density_estimator_model='mdn',
+                               restart=True,
+                               sbi_device='cpu',
+                               gl=[10*nS, 100*nS],
+                               C=[0.1*nF, 10*nF])
+
+Remarks
+-------
+
+For a better understanding, please go through examples that go step-by-step
+through the entire process. Currently, there are two tutorials: the one that
+is covering `simple interface <https://brian2modelfitting.readthedocs.io/en/stable/examples/hh_sbi_simple.html>`_,
+appropriate for the regular user, and the one that goes a bit more in-depth by
+using `flexible interface <https://brian2modelfitting.readthedocs.io/en/stable/examples/hh_sbi_flex.html>`_,
+and shows how to manually go through the process of inference, storing/loading
+the training data and the trained neural density estimator, parameter space
+visualization, conditioning, etc.
\ No newline at end of file
diff --git a/docs_sphinx/introduction/index.rst b/docs_sphinx/introduction/index.rst
index fedf7e19..139929e2 100644
--- a/docs_sphinx/introduction/index.rst
+++ b/docs_sphinx/introduction/index.rst
@@ -7,3 +7,4 @@ Model Fitting
    introduction
    howitworks
    tutorial_hh
+   tutorial_sbi
diff --git a/docs_sphinx/introduction/introduction.rst b/docs_sphinx/introduction/introduction.rst
index ae1feb0f..a849dc8e 100644
--- a/docs_sphinx/introduction/introduction.rst
+++ b/docs_sphinx/introduction/introduction.rst
@@ -6,11 +6,19 @@ The `.brian2modelfitting` toolbox provides three optimization classes:
  - `~brian2modelfitting.fitter.SpikeFitter`
  - `~brian2modelfitting.fitter.OnlineTraceFitter`
 
-The classes expect a model and data as an input and returns the best fit of
-parameters and the corresponding error. The toolbox can optimize over multiple
-traces (e.g. input currents) at the same time.
+and a simulation-based inference class:
+ - `~brian2modelfitting.inferencer.Inferencer`
 
-In following documentation we assume that ``brian2modelfitting`` has been imported like this:
+All classes expect a model and the data as an input and return either the best
+fit of each parameter with the corresponding error, or a posterior
+distribution over unknown parameters.
+The toolbox can optimize over multiple traces (e.g. input currents) at the
+same time.
+It also allows the possiblity of simultaneous fitting/inferencing by taking
+into account multiple output variables including spike trains.
+
+In following documentation we assume that `.brian2modelfitting` has been
+installed and imported as follows:
 
 .. code:: python
 
@@ -20,10 +28,9 @@ In following documentation we assume that ``brian2modelfitting`` has been import
 Installation
 ------------
 
-To install Model Fitting alongside Brian2 you can use pip, by using
-a pip utility:
+To install the toolbox alongside Brian 2 simulator, use ``pip`` as follows:
 
-.. code:: python
+.. code::
 
   pip install brian2modelfitting
 
@@ -33,9 +40,9 @@ Testing Model Fitting
 
 Version on master branch gets automatically tested with Travis services.
 To test the code yourself, you will need to have ``pytest`` installed and run
-the command:
+the following command inside the `.brian2modelfitting` root directory:
 
 
-.. code:: python
+.. code::
 
     pytest
diff --git a/docs_sphinx/introduction/tutorial_sbi.rst b/docs_sphinx/introduction/tutorial_sbi.rst
new file mode 100644
index 00000000..a848499a
--- /dev/null
+++ b/docs_sphinx/introduction/tutorial_sbi.rst
@@ -0,0 +1,725 @@
+Tutorial: Inferencer
+====================
+
+In this tutorial, we use simulation-based inference on the Hodgkin-Huxley
+neuron model, [Hodgkin1952]_, where different scenarios are considered. The
+data are synthetically generated from the model which enables the comparisson
+of the optimized parameter values with the ground truth parameter values as
+used in the generative model.
+
+We start with importing basic packages and modules. Note that
+``tutorial_sbi_helpers`` contains three different functions that are used
+extensively throughout this tutorial. Namely, ``plot_traces`` is used for
+visualization of the input current and output voltage traces, and optionally
+for visualization of spike trains or sampled traces obtained by using an
+estimated posterior. In order to detect and capture spike events in a given
+voltage trace, we use ``spike_times`` function. Finally,
+``plot_cond_coeff_mat`` is used to visualize conditional coerrelation matrix.
+For a detailed overview of the mechanism of each of the functions, download
+the script: :download:`tutorial_sbi_helpers.py <../../examples/tutorial_sbi_helpers.py>`
+
+.. code:: python
+
+    from brian2 import *
+    from brian2modelfitting import Inferencer
+    from scipy.stats import kurtosis as kurt
+    from tutorial_sbi_helpers import (spike_times,
+                                      plot_traces,
+                                      plot_cond_coeff_mat)
+
+Now, let's load the input current and output voltage traces by using ``NumPy``.
+Note that all functions available in ``NumPy``, as well as in ``Matplotlib``,
+are implicitly available after ``brian2``  was imported.
+
+.. code:: python
+
+    inp_traces = load('input_traces_synthetic.npy').reshape(1, -1)
+    out_traces = load('output_traces_synthetic.npy').reshape(1, -1)
+
+The data is generated by running ``generate_traces_synthetic.py`` script:
+:download:`generate_traces_synthetic.py <../../examples/generate_traces_synthetic.py>`
+
+By setting the time step, we are able to set up the time domain. From the
+array of time steps, we also define the exact time when the stimulus starts
+and when it ends by computing ``stim_start`` and ``stim_end``, respectively.
+This will come handy later during the feature extraction process.
+
+.. code:: python
+
+    dt = 0.05*ms
+    t = arange(0, inp_traces.size*dt/ms, dt/ms)
+    stim_start, stim_end = t[where(inp_traces[0, :] != 0)[0][[0, -1]]]
+
+By calling ``plot_traces`` and passing the array of time steps, the input
+current and output voltage traces, we obtain the visualization of a synthetic
+neural activity recordings:
+
+.. code:: python
+
+    fig, ax = plot_traces(t, inp_traces, out_traces)
+
+.. image:: ../_static/hh_sbi_recording_traces.png
+    :width: 70 %
+
+
+Toy-example: infer two free parameters
+--------------------------------------
+
+The first scenario we cover is a simple inference procedure of two unknown
+parameters in the Hodgkin-Huxley neuron model. The parameters to infer are
+the maximal value of sodium, :math:`\overline{g}_\mathrm{Na}`, and potassium
+electrical conductances, :math:`\overline{g}_\mathrm{K}`.
+
+By following the standard practice from Brian 2 simulator, we have to define
+parameters of the model, initial conditions for differential equations that
+describe the model, and the model itself:
+
+.. code:: python
+
+    # set parameters of the model
+    E_Na = 53*mV
+    E_K = -107*mV
+    E_l = -70*mV
+    VT = -60.0*mV
+    g_l = 10*nS
+    Cm = 200*pF
+
+    # set ground truth parameters, which are unknown from the model's perspective
+    ground_truth_params = {'g_Na': 32*uS,
+                           'g_K': 1*uS}
+
+    # define initial conditions
+    init_conds = {'v': 'E_l',
+                  'm': '1 / (1 + beta_m / alpha_m)',
+                  'h': '1 / (1 + beta_h / alpha_h)',
+                  'n': '1 / (1 + beta_n / alpha_n)'}
+
+    # define the Hodgkin-Huxley neuron model
+    eqs = '''
+        # non-linear set of ordinary differential equations
+        dv/dt = - (g_Na * m ** 3 * h * (v - E_Na)
+                   + g_K * n ** 4 * (v - E_K)
+                   + g_l * (v - E_l) - I) / Cm : volt
+        dm/dt = alpha_m * (1 - m) - beta_m * m : 1
+        dn/dt = alpha_n * (1 - n) - beta_n * n : 1
+        dh/dt = alpha_h * (1 - h) - beta_h * h : 1
+        
+        # time independent rate constants for a channel activation/inactivation
+        alpha_m = ((-0.32 / mV) * (v - VT - 13.*mV))
+                   / (exp((-(v - VT - 13.*mV)) / (4.*mV)) - 1) / ms : Hz
+        beta_m = ((0.28/mV) * (v - VT - 40.*mV))
+                  / (exp((v - VT - 40.*mV) / (5.*mV)) - 1) / ms : Hz
+        alpha_h = 0.128 * exp(-(v - VT - 17.*mV) / (18.*mV)) / ms : Hz
+        beta_h = 4 / (1 + exp((-(v - VT - 40.*mV)) / (5.*mV))) / ms : Hz
+        alpha_n = ((-0.032/mV) * (v - VT - 15.*mV))
+                   / (exp((-(v - VT - 15.*mV)) / (5.*mV)) - 1) / ms : Hz
+        beta_n = 0.5 * exp(-(v - VT - 10.*mV) / (40.*mV)) / ms : Hz
+            
+        # free parameters
+        g_Na : siemens (constant)
+        g_K : siemens (constant)
+        '''
+
+Since the output of the model is extremely high-dimensional, and since we are
+interested only in a few hand-picked features that will capture the gist of
+the neuronal activity, we start the inference process by defining a list of
+summary functions.
+
+Each summary feature is obtained by calling a single-valued function on each
+trace generated by using a sampled prior distribution over unknown parameters.
+In this case, we consider the maximal value, mean and standard deviatian of
+action potential, and the membrane resting potential.
+
+.. code:: python
+
+    v_features = [
+        # max action potential
+        lambda x: np.max(x[(t > stim_start) & (t < stim_end)]),
+        # mean action potential
+        lambda x: np.mean(x[(t > stim_start) & (t < stim_end)]),
+        # std of action potential
+        lambda x: np.std(x[(t > stim_start) & (t < stim_end)]),
+        # membrane resting potential
+        lambda x: np.mean(x[(t > 0.1 * stim_start) & (t < 0.9 * stim_start)])
+    ]
+
+
+Inferencer
+~~~~~~~~~~
+
+The minimum set of arguments for the `~brian2modelfitting.inferencer.Inferencer`
+class constructor are the time step, ``dt``, input data traces, ``input``,
+output data traces, ``output``, and the model that will be used for the
+inference process, ``model``. Input and output traces should have the number
+of rows that corresponds to the number of observed traces, and the number of
+columns should be equal to the number of time steps in each trace. 
+
+Here, we define additional arguments such as: ``method`` to define an
+integration technique used for solving the set of differential equations,
+``threshold`` to define a condition that produce a single spike, ``refractory``
+to define a condition under which a neuron remains refractory, and
+``param_init`` to define a set of initial conditions. We also define the set of
+summary features that is used to represent the data instead of using the entire
+trace. Summary features are passed to the inference algorithm via ``features``
+argument.
+
+.. code:: python
+
+    inferencer = Inferencer(dt=dt, model=eqs,
+                            input={'I': inp_traces*amp},
+                            output={'v': out_traces*mV},
+                            features={'v': v_features},
+                            method='exponential_euler',
+                            threshold='m > 0.5',
+                            refractory='m > 0.5',
+                            param_init=init_conds)
+
+After the ``inferencer`` is instantiated, we may begin the inference process
+by calling `~brian2modelfitting.inferencer.Inferencer.infer` and defining the
+total number of samples that are used for the training of a neural density
+estimator. We use the sequential neural posterior estimation algorithm (SNPE),
+proposed in [Greenberg2019]_. 
+
+
+Posterior
+~~~~~~~~~
+
+Neural density estimator learns the probabilistic mapping of the input data,
+i.e., sampled parameter values given a prior distribution, and the output data,
+i.e., summary features extracted from the traces, obtained by solving the model
+with the corresponding set of sampled parameters from the input data.
+
+We can choose the inference method and the estimator model , but only
+arguments that `~brian2modelfitting.inferencer.Inferencer.infer` requires are
+the number of samples (in case of running the inference process for the first
+time), ``n_samples``, and upper and lower bounds for each unknown parameter.
+
+.. code:: python
+
+    posterior = inferencer.infer(n_samples=15_000,
+                                 n_rounds=1,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS])
+
+After inference is completed, the estimated posterior distribution can be
+analyzed by observing the pairwise relationship between each pair of
+parameters. But before, we have to draw samples from the estimated posterior
+as follows:
+
+.. code:: python
+
+    samples = inferencer.sample((10_000, ))
+
+The samples are stored inside the `~brian2modelfitting.inferencer.Inferencer`
+object and are available through ``samples`` variable.
+We create a visual representation of the pairwise relationship of the posterior
+as follows:
+
+.. code:: python
+
+    limits = {'g_Na': [1*uS, 100*uS],
+              'g_K': [0.1*uS, 10*uS]}
+    labels = {'g_Na': r'$\overline{g}_{Na}$',
+              'g_K': r'$\overline{g}_{K}$'}
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+
+.. image:: ../_static/hh_sbi_toy_pairplot.png
+    :width: 70 %
+
+Let's zoom in a bit:
+
+.. image:: ../_static/hh_sbi_toy_pairplot_zoom_in.png
+    :width: 70 %
+
+The inferred posterior is plotted against the ground truth parameters, and as
+can be seen, the ground truth parameters are located in high-probability
+regions of the estimated distribution.
+
+To further substantiate this, let's now see the traces simulated from a single
+set of parameters sampled from the posterior:
+
+.. code:: python
+
+    inf_traces = inferencer.generate_traces()
+
+We again use the ``plot_traces`` function as follows:
+
+.. code:: python
+
+    fig, ax = plot_traces(t, inp_traces, out_traces, inf_traces=array(inf_traces/mV))
+
+.. image:: ../_static/hh_sbi_toy_inf_traces.png
+    :width: 70 %
+
+
+Additional free parameters
+--------------------------
+
+The simple scenarios where only 2 parameters are considered works quite well
+using synthetic data traces. What if we have a larger number of unkown
+parameters? Let's now consider additional unkown parameters for the same model
+as before. In addition to the unknown maximal values of the
+electrical conductance of the sodium and potassium channels, the membrane
+capacity, :math:`C_m`, and the maximal value of the electrical conductance of
+the leakage ion channel, :math:`\overline{g}_{l}`, are also unknown.
+
+We can try to do the same as before with a bit more training data:
+
+.. code:: python
+    
+    del Cm, g_l
+
+    # set parameters, initial condition and the model
+    E_Na = 53*mV
+    E_K = -107*mV
+    E_l = -70*mV
+    VT = -60.0*mV
+
+    ground_truth_params = {'g_Na': 32*uS,
+                           'g_K': 1*uS,
+                           'g_l': 10*nS,
+                           'Cm': 200*pF}
+
+    init_conds = {'v': 'E_l',
+                  'm': '1 / (1 + beta_m / alpha_m)',
+                  'h': '1 / (1 + beta_h / alpha_h)',
+                  'n': '1 / (1 + beta_n / alpha_n)'}
+
+    eqs = '''
+        # non-linear set of ordinary differential equations
+        dv/dt = - (g_Na * m ** 3 * h * (v - E_Na)
+                   + g_K * n ** 4 * (v - E_K)
+                   + g_l * (v - E_l) - I) / Cm : volt
+        dm/dt = alpha_m * (1 - m) - beta_m * m : 1
+        dn/dt = alpha_n * (1 - n) - beta_n * n : 1
+        dh/dt = alpha_h * (1 - h) - beta_h * h : 1
+        
+        # time independent rate constants for activation and inactivation
+        alpha_m = ((-0.32 / mV) * (v - VT - 13.*mV))
+                   / (exp((-(v - VT - 13.*mV)) / (4.*mV)) - 1) / ms : Hz
+        beta_m = ((0.28/mV) * (v - VT - 40.*mV))
+                  / (exp((v - VT - 40.*mV) / (5.*mV)) - 1) / ms : Hz
+        alpha_h = 0.128 * exp(-(v - VT - 17.*mV) / (18.*mV)) / ms : Hz
+        beta_h = 4 / (1 + exp((-(v - VT - 40.*mV)) / (5.*mV))) / ms : Hz
+        alpha_n = ((-0.032/mV) * (v - VT - 15.*mV))
+                   / (exp((-(v - VT - 15.*mV)) / (5.*mV)) - 1) / ms : Hz
+        beta_n = 0.5 * exp(-(v - VT - 10.*mV) / (40.*mV)) / ms : Hz
+            
+        # free parameters
+        g_Na : siemens (constant)
+        g_K : siemens (constant)
+        g_l : siemens (constant)
+        Cm : farad (constant)
+        '''
+
+    # infer the posterior using the same configuration as before
+    inferencer = Inferencer(dt=dt, model=eqs,
+                            input={'I': inp_traces*amp},
+                            output={'v': out_traces*mV},
+                            features={'v': v_features},
+                            method='exponential_euler',
+                            threshold='m > 0.5',
+                            refractory='m > 0.5',
+                            param_init=init_conds)
+
+    posterior = inferencer.infer(n_samples=20_000,
+                                 n_rounds=1,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS],
+                                 g_l=[1*nS, 100*nS],
+                                 Cm=[20*pF, 2*nF])
+
+    # finally, sample and visualize the posterior distribution
+    samples = inferencer.sample((10_000, ))
+
+    limits = {'g_Na': [1*uS, 100*uS],
+              'g_K': [0.1*uS, 10*uS],
+              'g_l': [1*nS, 100*nS],
+              'Cm': [20*pF, 2*nF]}
+    labels = {'g_Na': r'$\overline{g}_{Na}$',
+              'g_K': r'$\overline{g}_{K}$',
+              'g_l': r'$\overline{g}_{l}$',
+              'Cm': r'$C_{m}$'}
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+
+.. image:: ../_static/hh_sbi_4params_4features_pairplot.png
+    :width: 70 %
+
+This could have been expected. The posterior distribution is estimated poorly
+using a simple approach as in the toy example.
+
+Yes, we can play around with hyper-parameters and tuning the neural density
+estimator, but with this apporach we will not get far.
+
+We can, however, try with the non-amortized (or focused) approach, meaning we
+perform multi-round inference, where each following round will use the
+posterior from the previous one to sample new input data for the training,
+rather than using the same prior distribution as defined in the beginning.
+This approach yields additional advantage - the number of samples may be
+considerably lower, but it will lead to the posterior that is no longer being
+amortized - it is accurate only for a specific observation.
+
+.. code:: python
+
+    # note that the only difference is the number of rounds of inference
+    posterior = inferencer.infer(n_samples=10_000,
+                                 n_rounds=2,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 restart=True,
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS],
+                                 g_l=[1*nS, 100*nS],
+                                 Cm=[20*pF, 2*nF])
+
+    samples = inferencer.sample((10_000, ))
+
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+                        
+.. image:: ../_static/hh_sbi_4params_multiround_4features_pairplot.png
+    :width: 70 %
+
+This seems as a promising approach for parameters that already have the
+high-probability regions of the posterior distribution around ground-truth
+values. For other parameters, this leads to further deterioration of posterior
+estimates.
+
+So, we may wonder, how else can we improve the neural density estimator
+accuracy?
+
+Currently, we use only four features to describe extremely complex output of a
+neural model and should probably create a more comprehensive and more
+descriptive set of summary features. If we want to include data related to
+spikes in summary statistics, it is necessary to perform multi-objective
+optimization since we will observe spike trains as an output in addition to
+voltage traces.
+
+Multi-objective optimization
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In order to use spikes, we have to have some observation to pass to the
+`~brian2modelfitting.inferencer.Inferencer`. We can utilize ``spike_times`` as
+follows:
+
+.. code:: python
+
+    spike_times_list = [spike_times(t, out_trace) for out_trace in out_traces]
+
+To visually prove that spike times are indeed correct, we use ``plot_traces``:
+
+.. code:: python
+
+    fig, ax = plot_traces(t, inp_traces, out_traces, spike_times_list[0])
+
+.. image:: ../_static/hh_sbi_4params_spike_train.png
+    :width: 70 %
+
+Now, let's create additional features that will be applied to voltage traces,
+and a few features that will be applied to spike trains:
+
+.. code:: python
+
+    def voltage_deflection(x):
+        voltage_base = np.mean(x[t < stim_start])
+        stim_end_idx = np.where(t >= stim_end)[0][0]
+        steady_state_voltage_stimend = np.mean(x[stim_end_idx-10:stim_end_idx-5])
+        return steady_state_voltage_stimend - voltage_base
+
+
+    v_features = [
+        # max action potential
+        lambda x: np.max(x[(t > stim_start) & (t < stim_end)]),
+        # mean action potential
+        lambda x: np.mean(x[(t > stim_start) & (t < stim_end)]),
+        # std of action potential
+        lambda x: np.std(x[(t > stim_start) & (t < stim_end)]),
+        # kurtosis of action potential
+        lambda x: kurt(x[(t > stim_start) & (t < stim_end)], fisher=False),
+        # membrane resting potential
+        lambda x: np.mean(x[(t > 0.1 * stim_start) & (t < 0.9 * stim_start)]),
+        # the voltage deflection between base and steady-state voltage
+        voltage_deflection,
+    ]
+
+    s_features = [
+        # number of spikes in a train
+        lambda x: x.size,
+        # mean inter-spike interval
+        lambda x: 0. if np.diff(x).size == 0 else np.mean(diff(x)),
+        # time to first spike
+        lambda x: 0. if x.size == 0 else x[0]
+    ]
+
+The rest of the inference process stays pretty much the same as before:
+
+.. code:: python
+
+    inferencer = Inferencer(dt=dt, model=eqs,
+                            input={'I': inp_traces*amp},
+                            output={'v': out_traces*mV, 'spikes': spike_times_list},
+                            features={'v': v_features, 'spikes': s_features},
+                            method='exponential_euler',
+                            threshold='m > 0.5',
+                            refractory='m > 0.5',
+                            param_init=init_conds)
+
+    posterior = inferencer.infer(n_samples=20_000,
+                                 n_rounds=1,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS],
+                                 g_l=[1*nS, 100*nS],
+                                 Cm=[20*pF, 2*nF])
+
+    samples = inferencer.sample((10_000, ))
+
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+
+.. image:: ../_static/hh_sbi_4params_9features_pairplot.png
+    :width: 70 %
+
+Let's also visualize the sampled trace, this time using the mean of ten
+thousands drawn samples:
+
+.. code:: python
+
+    inf_traces = inferencer.generate_traces(n_samples=10_000, output_var='v')
+
+    fig, ax = plot_traces(t, inp_traces, out_traces, inf_traces=array(inf_traces/mV))
+
+.. image:: ../_static/hh_sbi_4params_inf_traces.png
+    :width: 70 %
+
+Okay, now we are clearly getting somewhere and this should be a strong
+indicatior of the importance of crafting quality summary statistics.
+
+Still, the summary statistics can be a huge bottleneck and can set back the
+training of a neural density estimator. For this reason automatic feature
+extraction can be considered instead.
+
+
+Automatic feature extraction
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To enable automatic feature extraction, ``features`` argument simly should
+not be defined when instantiating an inferencer object. And that's it.
+Everything else happens behind the scenes without any need for additional user
+intervention. If the user wants to gain additional control over the extraction
+process, in addition to changing the hyperparameters, they can also define
+their own embedding neural network.
+
+Default settings
+################
+
+.. code:: python
+
+    inferencer = Inferencer(dt=dt, model=eqs,
+                            input={'I': inp_traces*amp},
+                            output={'v': out_traces*mV},
+                            method='exponential_euler',
+                            threshold='m > 0.5',
+                            refractory='m > 0.5',
+                            param_init=init_conds)
+
+    posterior = inferencer.infer(n_samples=20_000,
+                                 n_rounds=1,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS],
+                                 g_l=[1*nS, 100*nS],
+                                 Cm=[20*pF, 2*nF])
+
+    samples = inferencer.sample((10_000, ))
+
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+
+.. image:: ../_static/hh_sbi_4params_default_automatic_ext_pairplot.png
+    :width: 70 %
+
+
+Custom embedding network
+########################
+
+Here, we demonstrate how to build a custom summary feature extractor and how
+to exploit the GPU processing power to speed up the inference process.
+
+Note that the use of the GPU will result in the speed-up of computation time
+only if a custom automatic feature extractor uses techniques that are actually
+faster to compute on the GPU.
+
+For this case, we use the YuleNet, a convolutional neural network, proposed in
+[Rodrigues2020]_. The authors outline impresive results where the automatic
+feature extraction by using the YuleNet is capable of outperforming carefully
+hand-crafted features.
+
+.. code:: python
+
+    import torch
+    from torch import nn
+
+
+    class YuleNet(nn.Module):
+        """The summary feature extractor proposed in Rodrigues 2020.
+        
+        Parameters
+        ----------
+        in_features : int
+            Number of input features should correspond to the size of a
+            single output voltage trace.
+        out_features : int
+            Number of the features that are used for the inference process.
+            
+        Returns
+        -------
+        None
+        
+        References
+        ----------
+        * Rodrigues, P. L. C. and Gramfort, A. "Learning summary features
+          of time series for likelihood free inference" 3rd Workshop on
+          Machine Learning and the Physical Sciences (NeurIPS 2020). 2020.
+        """
+        def __init__(self, in_features, out_features):
+            super().__init__()
+            self.conv1 = nn.Conv1d(in_channels=1, out_channels=8, kernel_size=64,
+                                   stride=1, padding=32, bias=True)
+            self.relu1 = nn.ReLU()
+            pooling1 = 16
+            self.pool1 = nn.AvgPool1d(kernel_size=pooling1)
+
+            self.conv2 = nn.Conv1d(in_channels=8, out_channels=8, kernel_size=64,
+                                   stride=1, padding=32, bias=True)
+            self.relu2 = nn.ReLU()
+            pooling2 = int((in_features // pooling1) // 16)
+            self.pool2 = nn.AvgPool1d(kernel_size=pooling2)
+
+            self.dropout = nn.Dropout(p=0.50)
+            linear_in = 8 * in_features // (pooling1 * pooling2) - 1
+            self.linear = nn.Linear(in_features=linear_in,
+                                    out_features=out_features)
+            self.relu3 = nn.ReLU()
+
+        def forward(self, x):
+            if x.ndim == 1:
+                x = x.view(1, 1, -1)
+            else:
+                x = x.view(len(x), 1, -1)
+            x_conv1 = self.conv1(x)
+            x_relu1 = self.relu1(x_conv1)
+            x_pool1 = self.pool1(x_relu1)
+
+            x_conv2 = self.conv2(x_pool1)
+            x_relu2 = self.relu2(x_conv2)
+            x_pool2 = self.pool2(x_relu2)
+            
+            x_flatten = x_pool2.view(len(x), 1, -1)
+            x_dropout = self.dropout(x_flatten)
+
+            x = self.relu3(self.linear(x_dropout))
+            return x.view(len(x), -1)
+
+In the following code, we also demonstrate how to control the hyperparameters
+of the density estimator by using additional keyword arguments in `~brian2modelfitting.inferencer.Inferencer.infer`
+method:
+
+.. code:: python
+
+    in_features = out_traces.shape[1]
+    out_features = 10
+
+    inferencer = Inferencer(dt=dt, model=eqs,
+                            input={'I': inp_traces*amp},
+                            output={'v': out_traces*mV},
+                            method='exponential_euler',
+                            threshold='m > 0.5',
+                            refractory='m > 0.5',
+                            param_init=init_conds)
+
+    posterior = inferencer.infer(n_samples=20_000,
+                                 n_rounds=1,
+                                 inference_method='SNPE',
+                                 density_estimator_model='maf',
+                                 inference_kwargs={'embedding_net': YuleNet(in_features, out_features)},
+                                 train_kwargs={'num_atoms': 10,
+                                              'training_batch_size': 100,
+                                              'use_combined_loss': True,
+                                              'discard_prior_samples': True},
+                                 sbi_device='gpu',
+                                 g_Na=[1*uS, 100*uS],
+                                 g_K=[0.1*uS, 10*uS],
+                                 g_l=[1*nS, 100*nS],
+                                 Cm=[20*pF, 2*nF])
+
+    samples = inferencer.sample((10_000, ))
+
+    fig, ax = inferencer.pairplot(limits=limits,
+                                  labels=labels,
+                                  ticks=limits,
+                                  points=ground_truth_params,
+                                  points_offdiag={'markersize': 5},
+                                  points_colors=['C3'],
+                                  figsize=(6, 6))
+
+.. image:: ../_static/hh_sbi_4params_yulenet_pairplot.png
+    :width: 70 %
+
+
+Next steps
+----------
+
+To learn more read the reference API and check out more examples available
+:doc:`here <../examples/index>`.
+
+References
+----------
+
+.. [Hodgkin1952] Hodgkin, A. L., and Huxley, A. F. "A quantitative description
+                 of membrane current and its application to conduction and
+                 excitation in nerve" Journal of Physiology 117(4):500-44. 1952.
+                 doi: `10.1113/jphysiol.1952.sp004764 <https://doi.org/10.1113/jphysiol.1952.sp004764>`_
+
+.. [Greenberg2019] Greenberg, D. S., Nonnenmacher, M. et al. "Automatic
+                   posterior transformation for likelihood-free inference"
+                   36th International Conference on Machine Learning
+                   (ICML 2019). 2019. Available online: `arXiv:1905.07488 <https://arxiv.org/abs/1905.07488>`_
+
+.. [Rodrigues2020] Rodrigues, P. L. C. and Gramfort, A. "Learning summary
+                   features of time series for likelihood free inference" 3rd
+                   Workshop on Machine Learning and the Physical Sciences
+                   (NeurIPS 2020). 2020. Available online: `arXiv:2012.02807 <https://arxiv.org/abs/2012.02807>`_
\ No newline at end of file
diff --git a/examples/IF_sbi.py b/examples/IF_sbi.py
new file mode 100644
index 00000000..3dc182f8
--- /dev/null
+++ b/examples/IF_sbi.py
@@ -0,0 +1,110 @@
+from brian2 import *
+from brian2modelfitting import *
+
+
+# Set parameters for integrate-and-fire model
+dt = 0.1*ms
+sim_time = 60*ms
+El = -70*mV
+VT = -50*mV
+DeltaT = 2*mV
+ground_truth_params = {'gl': 30*nS, 'C': 1*nF}  # parameters to be inferred
+
+# Set input trace
+inp_trace = hstack([zeros(int(round(5*ms/dt))),
+                    ones(int(round(25*ms/dt))),
+                    zeros(int(round(30*ms/dt)))]) * 10*nA
+I = TimedArray(inp_trace, dt=dt)
+
+# Set model equations
+eqs = '''
+    dv/dt = (gl*(El - v) + gl*DeltaT*exp((v - VT) / DeltaT) + I(t)) / C : volt
+    gl    : siemens (constant)
+    C     : farad (constant)
+    '''
+
+# Run a model and create synthetic voltage traces and obtain spike trains
+neurons = NeuronGroup(1, eqs, dt=dt,
+                      threshold='v > -50 * mV',
+                      reset='v = -70 * mV',
+                      method='exponential_euler')
+neurons.v = -70*mV
+neurons.set_states(ground_truth_params)
+state_monitor = StateMonitor(neurons, 'v', record=True)
+spike_monitor = SpikeMonitor(neurons, record=True)
+run(sim_time)
+out_trace = state_monitor.v  # voltage traces in mV
+spike_times = array(spike_monitor.t_)  # spike events in s
+
+# Inference
+start_scope()
+El = -70*mV
+VT = -50*mV
+DeltaT = 2*mV
+eqs_inf = '''
+    dv/dt = (gl*(El - v) + gl*DeltaT*exp((v - VT) / DeltaT) + I) / C : volt
+    gl    : siemens (constant)
+    C     : farad (constant)
+    '''
+
+# Instantiate the inferencer object
+spike_features_list = [
+    lambda x: x.size,  # number of spikes
+    lambda x: 0. if diff(x).size == 0 else mean(diff(x)),  # mean ISI
+    ]
+inferencer = Inferencer(dt=dt, model=eqs_inf,
+                        input={'I': inp_trace.reshape(1, -1)},
+                        output={'spikes': [spike_times]},
+                        features={'spikes': spike_features_list},
+                        method='exponential_euler',
+                        threshold='v > -50*mV',
+                        reset='v = -70*mV',
+                        param_init={'v': -70*mV})
+
+# Infer parameter posteriors for given parameters
+posterior = inferencer.infer(n_samples=3_000,
+                             n_rounds=3,
+                             gl=[10*nS, 100*nS],
+                             C=[0.1*nF, 10*nF])
+
+# Sample from the posterior
+samples = inferencer.sample((3_000, ))
+
+# Create a pairplot
+labels = {'gl': r'$\overline{g}_{l}$', 'C': r'$C$'}
+ticks = {'gl': [10*nS, 100*nS], 'C': [0.1*nF, 10*nF]}
+inferencer.pairplot(labels=labels,
+                    limits=ticks,
+                    ticks=ticks,
+                    points=ground_truth_params,
+                    points_offdiag={'markersize': 9},
+                    points_colors=['r'],
+                    figsize=(6, 6))
+
+# Visualize the trace by sampling trained posterior
+t = arange(0, inp_trace.size*dt/ms, dt/ms)
+t_start, t_end = t[where(inp_trace != 0)[0][[0, -1]]]
+
+# Identify spike events
+spike_v = (out_trace.min()/mV, out_trace.max()/mV)
+spike_i = []
+for spike_time in spike_times:
+    spike_i.append(where(isclose(spike_time * 1000, t))[0].item())
+
+# Generate traces from a single sample of parameters
+inf_trace = inferencer.generate_traces(output_var='v')
+nrows = 2
+ncols = 1
+fig, axs = subplots(nrows, ncols, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]}, figsize=(7, 3))
+axs[0].plot(t, out_trace.T/mV, 'C3-', lw=3, label='recordings')
+axs[0].plot(t, inf_trace.T/mV, 'k--', lw=2, label='sampled traces')
+axs[0].vlines(t[spike_i], *spike_v, lw=3, color='C0', label='spike', zorder=3)
+axs[0].set_ylabel('$V$, mV')
+axs[0].legend()
+axs[1].plot(t, inp_trace.T/nA, lw=3, c='k', label='stimuli')
+axs[1].set_xlabel('$t$, ms')
+axs[1].set_ylabel('$I$, nA')
+axs[1].legend()
+tight_layout()
+show()
diff --git a/examples/generate_traces_synthetic.py b/examples/generate_traces_synthetic.py
new file mode 100644
index 00000000..97e6fa6d
--- /dev/null
+++ b/examples/generate_traces_synthetic.py
@@ -0,0 +1,52 @@
+from brian2 import *
+
+
+defaultclock.dt = 0.1*ms
+
+
+I1 = 0.1*nA
+I2 = 0.5*nA
+
+t_on = 100*ms
+t_total = 1000*ms
+t_off = t_total - t_on
+
+C = 200*pF
+VT = -60.0*mV
+E_Na = 53*mV
+g_Na = 32*uS
+E_K = -107*mV
+g_K = 1*uS
+E_l = -70*mV
+g_l = 10*nS
+E_l = -70*mV
+
+eqs = '''
+    dVm/dt = - (g_Na * m ** 3 * h * (Vm - E_Na) 
+                + g_K * n ** 4 * (Vm - E_K)
+                + g_l * (Vm - E_l) - I_inj) / C : volt
+    I_inj = ((I2 - I1) / (t_off - t_on) * (t - t_on) + I1) * int((t >= t_on) and (t < t_off)) : amp (shared)
+    dm/dt = alpha_m * (1-m) - beta_m * m : 1
+    dn/dt = alpha_n * (1-n) - beta_n * n : 1
+    dh/dt = alpha_h * (1-h) - beta_h * h : 1
+    alpha_m = ((-0.32/mV) * (Vm - VT - 13.*mV)
+               / (exp((-(Vm - VT - 13.*mV)) / (4.*mV)) - 1))/ms : Hz
+    beta_m = ((0.28/mV) * (Vm - VT - 40.*mV)
+              / (exp((Vm - VT - 40.*mV) / (5.*mV)) - 1))/ms : Hz
+    alpha_h = 0.128 * exp(-(Vm - VT - 17.*mV) / (18.*mV))/ms : Hz
+    beta_h = 4 / (1 + exp((-(Vm - VT - 40.*mV)) / (5.*mV)))/ms : Hz
+    alpha_n = ((-0.032/mV) * (Vm - VT - 15.*mV)
+               / (exp((-(Vm - VT - 15.*mV)) / (5.*mV)) - 1))/ms : Hz
+    beta_n = 0.5 * exp(-(Vm - VT - 10.*mV) / (40.*mV))/ms : Hz
+    '''
+neurons = NeuronGroup(1, eqs, threshold='m>0.5', refractory='m>0.5',
+                      method='exponential_euler', name='neurons')
+neurons.Vm = 'E_l'
+neurons.m = '1 / (1 + beta_m / alpha_m)'
+neurons.h = '1 / (1 + beta_h / alpha_h)'
+neurons.n = '1 / (1 + beta_n / alpha_n)'
+Vm_mon = StateMonitor(neurons, ['Vm', 'I_inj'], record=True, name='Vm_mon')
+run(t_total)
+
+save('input_traces_synthetic.npy', Vm_mon.I_inj.ravel()/amp)
+save('output_traces_synthetic.npy', Vm_mon.Vm.ravel()/mV)
diff --git a/examples/hh_sbi_flexible_interface.py b/examples/hh_sbi_flexible_interface.py
new file mode 100644
index 00000000..bdc95b15
--- /dev/null
+++ b/examples/hh_sbi_flexible_interface.py
@@ -0,0 +1,121 @@
+import os
+
+from brian2 import *
+from brian2modelfitting import *
+import pandas as pd
+
+
+# Load input and output data traces
+df_inp_traces = pd.read_csv('input_traces_hh.csv')
+df_out_traces = pd.read_csv('output_traces_hh.csv')
+inp_traces = df_inp_traces.to_numpy()
+inp_traces = inp_traces[[0, 1, 3], 1:]
+out_traces = df_out_traces.to_numpy()
+out_traces = out_traces[[0, 1, 3], 1:]
+
+# Model and its parameters
+area = 20_000*um**2
+El = -65*mV
+EK = -90*mV
+ENa = 50*mV
+VT = -63*mV
+dt = 0.01*ms
+eqs = '''
+    dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
+    dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
+            (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
+            (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
+    dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
+            (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
+    dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
+    g_na : siemens (constant)
+    g_kd : siemens (constant)
+    gl   : siemens (constant)
+    Cm   : farad (constant)
+    '''
+
+# Time domain
+t = arange(0, out_traces.shape[1]*dt/ms, dt/ms)
+t_start, t_end = t[where(inp_traces[0, :] != 0)[0][[0, -1]]]
+
+# Step-by-step inference
+# Start with the regular instatiation of the class
+inferencer = Inferencer(dt=dt, model=eqs,
+                        input={'I': inp_traces*amp},
+                        output={'v': out_traces*mV},
+                        features={'v': [lambda x: max(x),
+                                        lambda x: mean(x[(t > t_start) & (t < t_end)]),
+                                        lambda x: std(x[(t > t_start) & (t < t_end)])]},
+                        method='exponential_euler',
+                        threshold='m > 0.5',
+                        refractory='m > 0.5',
+                        param_init={'v': 'VT'})
+
+# Training data generation
+# Initializing the prior
+prior = inferencer.init_prior(gl=[1e-09*siemens, 1e-07*siemens],
+                              g_na=[2e-06*siemens, 2e-04*siemens],
+                              g_kd=[6e-07*siemens, 6e-05*siemens],
+                              Cm=[0.1*uF*cm**-2*area, 2*uF*cm**-2*area])
+# Prepare training data
+path_to_data = __file__[:-3] + '_data.npz'
+if os.path.exists(path_to_data):
+    theta, x = inferencer.load_summary_statistics(path_to_data)
+else:
+    # Generate training data
+    theta = inferencer.generate_training_data(n_samples=10_000,
+                                              prior=prior)
+    # Extract summary stats
+    x = inferencer.extract_summary_statistics(theta)
+    # Save the data for later use
+    inferencer.save_summary_statistics(path_to_data, theta, x)
+
+# Amortized inference
+# Training the neural density estimator
+inference = inferencer.init_inference(inference_method='SNPE',
+                                      density_estimator_model='mdn',
+                                      prior=prior)
+# First round of inference where no observation data is set to posterior
+posterior_amortized = inferencer.infer_step(proposal=prior,
+                                            inference=inference,
+                                            theta=theta, x=x)
+# Storing the trained posterior without a default observation
+path_to_posterior = __file__[:-3] + '_posterior_amortized.pth'
+inferencer.save_posterior(path_to_posterior)
+
+# Sampling from the posterior given observations, ...
+inferencer.sample((10_000, ))
+# ...visualize the samples, ...
+labels = {'gl': r'$\overline{g}_\mathrm{l}$',
+          'g_na': r'$\overline{g}_\mathrm{Na}$',
+          'g_kd': r'$\overline{g}_\mathrm{K}$',
+          'Cm': r'$\overline{C}_{m}$'}
+inferencer.pairplot(labels=labels)
+
+# ...and optionally, continue the multiround inference using ``infer`` method
+posterior_multiround = inferencer.infer()
+inferencer.sample((10_000, ))
+inferencer.pairplot(labels=labels)
+
+# Generate traces from a single sample of parameters
+inf_traces = inferencer.generate_traces()
+nrows = 2
+ncols = out_traces.shape[0]
+fig, axs = subplots(nrows, ncols, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]},
+                    figsize=(ncols * 3, 3))
+for idx in range(ncols):
+    axs[0, idx].plot(t, out_traces[idx, :].T, 'C3-', lw=3, label='recordings')
+    axs[0, idx].plot(t, inf_traces[idx, :].T/mV, 'k--', lw=2,
+                     label='sampled traces')
+    axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimuli')
+    axs[1, idx].set_xlabel('$t$, ms')
+    if idx == 0:
+        axs[0, idx].set_ylabel('$V$, mV')
+        axs[1, idx].set_ylabel('$I$, nA')
+handles, labels = [(h + l) for h, l
+                   in zip(axs[0, idx].get_legend_handles_labels(),
+                   axs[1, idx].get_legend_handles_labels())]
+fig.legend(handles, labels)
+tight_layout()
+show()
diff --git a/examples/hh_sbi_simple_interface.py b/examples/hh_sbi_simple_interface.py
new file mode 100644
index 00000000..5fb40396
--- /dev/null
+++ b/examples/hh_sbi_simple_interface.py
@@ -0,0 +1,191 @@
+from brian2 import *
+from brian2modelfitting import *
+import pandas as pd
+from scipy.stats import kurtosis as kurt
+
+
+# Load input and output data traces
+df_inp_traces = pd.read_csv('input_traces_hh.csv')
+df_out_traces = pd.read_csv('output_traces_hh.csv')
+inp_traces = df_inp_traces.to_numpy()
+inp_traces = inp_traces[[0, 1, 3, 4], 1:]
+out_traces = df_out_traces.to_numpy()
+out_traces = out_traces[[0, 1, 3, 4], 1:]
+
+# Define model and its parameters
+area = 20_000*um**2
+El = -65*mV
+EK = -90*mV
+ENa = 50*mV
+VT = -63*mV
+dt = 0.01*ms
+eqs = '''
+    dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
+    dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
+            (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
+            (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
+    dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
+            (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
+    dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
+
+    g_na : siemens (constant)
+    g_kd : siemens (constant)
+    gl   : siemens (constant)
+    Cm   : farad (constant)
+    '''
+
+# Time domain
+t = arange(0, out_traces.shape[1]*dt/ms, dt/ms)
+t_start, t_end = t[where(inp_traces[0, :] != 0)[0][[0, -1]]]
+
+# Visualize the recordings
+nrows = 2
+ncols = out_traces.shape[0]
+fig, axs = subplots(nrows, ncols, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]}, figsize=(12, 3))
+for idx in range(ncols):
+    axs[0, idx].plot(t, out_traces[idx, :].T, lw=3, c='C3', label='recordings')
+    axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimulus')
+    axs[1, idx].set_xlabel('$t$, ms')
+    if idx == 0:
+        axs[0, idx].set_ylabel('$V$, mV')
+        axs[1, idx].set_ylabel('$I$, nA')
+handles, labels = [(h + l) for h, l
+                   in zip(axs[0, idx].get_legend_handles_labels(),
+                   axs[1, idx].get_legend_handles_labels())]
+fig.legend(handles, labels)
+tight_layout()
+
+
+# Obtain spike times manually from recordings
+def get_spike_times(x):
+    x = x.copy()
+    # put everything to -40 mV that is below -40 mV or has negative slope
+    ind = where(x < -0.04)
+    x[ind] = -0.04
+    ind = where(diff(x) < 0)
+    x[ind] = -0.04
+
+    # remaining negative slopes are at spike peaks
+    ind = where(diff(x) < 0)
+    spike_times = array(t)[ind]
+
+    # number of spikes
+    if spike_times.shape[0] > 0:
+        spike_times = spike_times[
+            append(1, diff(spike_times)) > 0.5
+        ]
+    return spike_times / 1000  # in seconds
+
+
+# store spike times for each trace into the list
+spike_times = [get_spike_times(out_trace) for out_trace in out_traces]
+
+# Visualize the recordings and spikes while stimulus is on
+nrows = 2
+ncols = out_traces.shape[0]
+fig, axs = subplots(nrows, ncols, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]}, figsize=(12, 3))
+for idx in range(ncols):
+    spike_idx = in1d(t, spike_times[idx] * 1000).nonzero()[0]
+    spike_v = (out_traces[idx, :].min(), out_traces[idx, :].max())
+    axs[0, idx].plot(t, out_traces[idx, :].T, lw=3, c='C3',
+                     label='recordings', zorder=1)
+    axs[0, idx].vlines(t[spike_idx], *spike_v, lw=3, color='C0',
+                       label='spike event', zorder=2)
+    axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimulus')
+    axs[1, idx].set_xlabel('$t$, ms')
+    if idx == 0:
+        axs[0, idx].set_ylabel('$V$, mV')
+        axs[1, idx].set_ylabel('$I$, nA')
+handles, labels = [(h + l) for h, l
+                   in zip(axs[0, idx].get_legend_handles_labels(),
+                   axs[1, idx].get_legend_handles_labels())]
+fig.legend(handles, labels)
+tight_layout()
+
+# Define features to create a summary statistics representation of traces
+v_features = [
+    lambda x: max(x[(t > t_start) & (t < t_end)]),  # AP max
+    lambda x: mean(x[(t > t_start) & (t < t_end)]),  # AP mean
+    lambda x: std(x[(t > t_start) & (t < t_end)]),  # AP std
+    lambda x: kurt(x[(t > t_start) & (t < t_end)], fisher=False),  # AP kurt
+    lambda x: mean(x[(t > .25 * t_start) & (t < .75 * t_start)]),  # resting
+    lambda x: mean(x[(t > t_end) & (t <= max(t))]),  # steady-state
+]
+s_features = [lambda x: x.size]  # the number of spikes in a spike train
+
+# Simulation-based inference object instantiation
+inferencer = Inferencer(dt=dt, model=eqs,
+                        input={'I': inp_traces*amp},
+                        output={'v': out_traces*mV, 'spikes': spike_times},
+                        features={'v': v_features, 'spikes': s_features},
+                        method='exponential_euler',
+                        threshold='m > 0.5',
+                        refractory='m > 0.5',
+                        param_init={'v': 'VT'})
+
+# Multi-round inference
+posterior = inferencer.infer(n_samples=10_000,
+                             n_rounds=2,
+                             inference_method='SNPE',
+                             density_estimator_model='mdn',
+                             gl=[1e-09*siemens, 1e-07*siemens],
+                             g_na=[2e-06*siemens, 2e-04*siemens],
+                             g_kd=[6e-07*siemens, 6e-05*siemens],
+                             Cm=[0.1*uF*cm**-2*area, 2*uF*cm**-2*area])
+
+# Draw samples from posterior
+samples = inferencer.sample((10_000, ))
+
+# Visualize conditional pairwise relationships between parameter distributions
+limits = {'gl': [1e-9*siemens, 1e-07*siemens],
+          'g_na': [2e-06*siemens, 2e-04*siemens],
+          'g_kd': [6e-07*siemens, 6e-05*siemens],
+          'Cm': [0.1*uF*cm**-2*area, 2*uF*cm**-2*area]}
+labels = {'gl': r'$\overline{g}_{l}$',
+          'g_na': r'$\overline{g}_{Na}$',
+          'g_kd': r'$\overline{g}_{K}$',
+          'Cm': r'$C_{m}$'}
+inferencer.pairplot(limits=limits,
+                    labels=labels,
+                    ticks=limits,
+                    figsize=(6, 6))
+condition = inferencer.sample((1, ))
+inferencer.conditional_pairplot(condition=condition,
+                                limits=limits,
+                                labels=labels,
+                                ticks=limits,
+                                figsize=(6, 6))
+
+# Construct and visualize conditional coefficient matrix
+cond_coeff_mat = inferencer.conditional_corrcoeff(condition=condition,
+                                                  limits=limits)
+fig, ax = subplots(1, 1, figsize=(4, 4))
+im = imshow(cond_coeff_mat, clim=[-1, 1])
+_ = fig.colorbar(im)
+
+# Generate traces and visualize from a single sample of parameters
+inf_traces = inferencer.generate_traces(n_samples=10_000, output_var='v')
+
+nrows = 2
+ncols = out_traces.shape[0]
+fig, axs = subplots(nrows, ncols, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]}, figsize=(12, 3))
+for idx in range(ncols):
+    spike_idx = in1d(t, spike_times[idx]).nonzero()[0]
+    spike_v = (out_traces[idx, :].min(), out_traces[idx, :].max())
+    axs[0, idx].plot(t, out_traces[idx, :].T, 'C3-', lw=3, label='recordings')
+    axs[0, idx].plot(t, inf_traces[idx, :].T/mV, 'k--', lw=2,
+                     label='sampled traces')
+    axs[1, idx].plot(t, inp_traces[idx, :].T/nA, lw=3, c='k', label='stimuli')
+    axs[1, idx].set_xlabel('$t$, ms')
+    if idx == 0:
+        axs[0, idx].set_ylabel('$V$, mV')
+        axs[1, idx].set_ylabel('$I$, nA')
+handles, labels = [(h + l) for h, l
+                   in zip(axs[0, idx].get_legend_handles_labels(),
+                   axs[1, idx].get_legend_handles_labels())]
+fig.legend(handles, labels)
+tight_layout()
+show()
diff --git a/examples/hh_sbi_synthetic_traces.py b/examples/hh_sbi_synthetic_traces.py
new file mode 100644
index 00000000..f24f5945
--- /dev/null
+++ b/examples/hh_sbi_synthetic_traces.py
@@ -0,0 +1,119 @@
+from brian2modelfitting import *
+from brian2 import *
+
+
+# Input and output data traces simulation
+dt = 0.05*ms
+t_on = 50*ms
+t_total = 350*ms
+t_off = t_total - t_on
+I = 1*nA
+gleak = 10*nS
+Eleak = -70*mV
+VT = -60*mV
+C = 200*pF
+ENa = 53*mV
+EK = -107*mV
+ground_truth_params = {'gNa': 32*uS, 'gK': 1*uS}
+eqs = '''
+    dVm/dt = -(gNa*m**3*h*(Vm - ENa) + gK*n**4*(Vm - EK) + gleak*(Vm - Eleak) - I_inj) / C : volt
+    I_inj = int(t >= t_on and t < t_off)*I : amp (shared)
+    dm/dt = alpham*(1-m) - betam*m : 1
+    dn/dt = alphan*(1-n) - betan*n : 1
+    dh/dt = alphah*(1-h) - betah*h : 1
+    alpham = (-0.32/mV) * (Vm - VT - 13.*mV) / (exp((-(Vm - VT - 13.*mV))/(4.*mV)) - 1)/ms : Hz
+    betam = (0.28/mV) * (Vm - VT - 40.*mV) / (exp((Vm - VT - 40.*mV)/(5.*mV)) - 1)/ms : Hz
+    alphah = 0.128 * exp(-(Vm - VT - 17.*mV) / (18.*mV))/ms : Hz
+    betah = 4/(1 + exp((-(Vm - VT - 40.*mV)) / (5.*mV)))/ms : Hz
+    alphan = (-0.032/mV) * (Vm - VT - 15.*mV) / (exp((-(Vm - VT - 15.*mV)) / (5.*mV)) - 1)/ms : Hz
+    betan = 0.5*exp(-(Vm - VT - 10.*mV) / (40.*mV))/ms : Hz
+
+    # parameters
+    gNa : siemens (constant)
+    gK : siemens (constant)
+    '''
+neurons = NeuronGroup(1, eqs, dt=dt,
+                      threshold='m>0.5', refractory='m>0.5',
+                      method='exponential_euler')
+state_monitor = StateMonitor(neurons, 'Vm', record=True)
+spike_monitor = SpikeMonitor(neurons, record=False)
+neurons.gNa = ground_truth_params['gNa']
+neurons.gK = ground_truth_params['gK']
+neurons.Vm = 'Eleak'
+neurons.m = '1/(1 + betam/alpham)'
+neurons.h = '1/(1 + betah/alphah)'
+neurons.n = '1/(1 + betan/alphan)'
+run(t_total)
+I_inj = ((state_monitor.t >= t_on) & (state_monitor.t < t_off)) * I
+t = state_monitor.t
+out_trace = state_monitor.Vm
+
+# Inference
+start_scope()
+
+# A set of equations for the Inferencer class
+eqs = '''
+    dVm/dt = -(gNa*m**3*h*(Vm - ENa) + gK*n**4*(Vm - EK) + gleak*(Vm - Eleak) - I_inj) / C : volt
+    dm/dt = alpham*(1-m) - betam*m : 1
+    dn/dt = alphan*(1-n) - betan*n : 1
+    dh/dt = alphah*(1-h) - betah*h : 1
+    alpham = (-0.32/mV) * (Vm - VT - 13.*mV) / (exp((-(Vm - VT - 13.*mV))/(4.*mV)) - 1)/ms : Hz
+    betam = (0.28/mV) * (Vm - VT - 40.*mV) / (exp((Vm - VT - 40.*mV)/(5.*mV)) - 1)/ms : Hz
+    alphah = 0.128 * exp(-(Vm - VT - 17.*mV) / (18.*mV))/ms : Hz
+    betah = 4/(1 + exp((-(Vm - VT - 40.*mV)) / (5.*mV)))/ms : Hz
+    alphan = (-0.032/mV) * (Vm - VT - 15.*mV) / (exp((-(Vm - VT - 15.*mV)) / (5.*mV)) - 1)/ms : Hz
+    betan = 0.5*exp(-(Vm - VT - 10.*mV) / (40.*mV))/ms : Hz
+
+    # parameters
+    gNa : siemens (constant)
+    gK : siemens (constant)
+    '''
+
+# Instantiating the inferencer object with input and output data traces
+inferencer = Inferencer(dt=dt, model=eqs,
+                        input={'I_inj': I_inj.reshape(1, -1)},
+                        output={'Vm': out_trace},
+                        method='exponential_euler',
+                        threshold='m>0.5', refractory='m>0.5',
+                        param_init={'Vm': 'Eleak',
+                                    'm': '1/(1 + betam/alpham)',
+                                    'h': '1/(1 + betah/alphah)',
+                                    'n': '1/(1 + betan/alphan)'})
+
+# Amortized inference performed with automatic feature extraction
+inferencer.infer(n_samples=10_000,
+                 inference_method='SNPE',
+                 density_estimator_model='mdn',
+                 gNa=[.5*uS, 80.*uS],
+                 gK=[1e-4*uS, 15.*uS])
+
+# Sample from posterior
+inferencer.sample((10_000, ))
+
+# Visualize estimated posterior distribution
+limits = {'gNa': [.5*uS, 80.*uS],
+          'gK': [1e-4*uS, 15.*uS]}
+labels = {'gNa': r'$\overline{g}_\mathrm{Na}$',
+          'gK': r'$\overline{g}_\mathrm{K}$'}
+inferencer.pairplot(limits=limits,
+                    ticks=limits,
+                    labels=labels,
+                    points=ground_truth_params,
+                    points_offdiag={'markersize': 6},
+                    points_colors=['r'],
+                    figsize=(6, 6))
+
+# Visualize sampled traces
+inf_trace = inferencer.generate_traces()
+
+fig, axs = subplots(2, 1, sharex=True,
+                    gridspec_kw={'height_ratios': [3, 1]}, figsize=(7, 3))
+axs[0].plot(t, out_trace.ravel()/mV, 'C3-', lw=3, label='simulated recordings')
+axs[0].plot(t, inf_trace.ravel()/mV, 'k--', lw=2, label='posterior sample')
+axs[0].set(ylabel='$V_m$, mV')
+axs[0].legend(loc='upper right')
+axs[1].plot(t, I_inj.ravel()/nA, 'k-', lw=3, label='stimulus')
+axs[1].set(xlabel='$t$, s', ylabel='I, nA')
+axs[1].legend(loc='upper right')
+tight_layout()
+show()
diff --git a/examples/tutorial_sbi_helpers.py b/examples/tutorial_sbi_helpers.py
new file mode 100644
index 00000000..422a4c29
--- /dev/null
+++ b/examples/tutorial_sbi_helpers.py
@@ -0,0 +1,107 @@
+from brian2 import *
+
+
+def spike_times(t, x):
+    """Return time points in seconds in which a spike occurs.
+
+    Parameters
+    ----------
+    t : numpy.ndarray
+        Time steps in ms.
+    x : numpy.ndarray
+        Voltage trace in mV.
+
+    Returns
+    -------
+    numpy.ndarray
+        Spike times in seconds.
+    """
+    x = x.copy()
+
+    ind = where(x < -0.04)
+    x[ind] = -0.04
+    ind = where(diff(x) < 0)
+    x[ind] = -0.04
+
+    ind = where(diff(x) < 0)
+    spike_times = array(t)[ind]
+
+    if spike_times.shape[0] > 0:
+        spike_times = spike_times[append(1, diff(spike_times)) > 0.5]
+    return spike_times / 1000
+
+
+def plot_traces(t, inp_traces, out_traces, spike_times=None, inf_traces=None):
+    """Visualize input current and output voltage traces.
+
+    Parameters
+    ----------
+    t : numpy.ndarray
+        Time steps in ms.
+    inp_traces : numpy.ndarray
+        Synaptic stimulus data in Amperes.
+    inp_traces : numpy.ndarray
+        Recorded voltage traces in mV.
+    spike_times : numpy.ndarray
+        Spike times in s.
+    inf_traces : numpy.ndarray, optional
+        Sampled voltage traces from the approximated posterior in mV.
+
+    Returns
+    -------
+    tuple
+        Figure and axes.
+    """
+    if (inp_traces.ndim, out_traces.ndim) != (2, 2):
+        inp_traces = inp_traces.reshape(1, -1)
+        out_traces = out_traces.reshape(1, -1)
+    if inf_traces is not None:
+        if inf_traces.ndim != 2:
+            inf_traces = inf_traces.reshape(1, -1)
+    if spike_times is not None:
+        spike_times = spike_times.ravel()
+        spike_i = []
+        for spike_time in spike_times:
+            spike_i.append(where(isclose(spike_time * 1000, t))[0].item())
+        spike_v = (out_traces.min(), out_traces.max())
+
+    fig, ax = subplots(2, 1, sharex=True,
+                       gridspec_kw={'height_ratios': [3, 1]}, figsize=(6, 3))
+
+    # voltage traces
+    ax[0].plot(t, out_traces.T, 'C0', lw=3, label='recordings', zorder=1)
+    if inf_traces is not None:
+        ax[0].plot(t, inf_traces.T, 'C1--', lw=2, label='sampled traces',
+                   zorder=2)
+    if spike_times is not None:
+        ax[0].vlines(t[spike_i], *spike_v, lw=2, color='C3', label='spikes',
+                     zorder=3)
+    ax[0].set(ylabel='$V$ [mV]')
+    ax[0].legend(loc='upper right')
+
+    # input stimulus current
+    ax[1].plot(t, inp_traces.T * 1e9, lw=2, c='k', label='stimulus')
+    ax[1].set(xlabel='$t$ [ms]', ylabel='$I$, nA')
+    ax[1].legend(loc='upper right')
+
+    tight_layout()
+    return fig, ax
+
+
+def plot_cond_coeff_mat(cond_coeff_mat):
+    """Visualize conditional coerrelation matrix.
+
+    Parameters
+    ----------
+    cond_coeff_mat : numpy.ndarray
+        Average conditional correlation matrix.
+
+    Returnss
+    -------
+    tuple
+        Figure and axes.
+    """
+    fig, ax = subplots(1, 1, figsize=(4, 4))
+    im = imshow(cond_coeff_mat, clim=[-1, 1])
+    _ = fig.colorbar(im)
+    return fig, ax
diff --git a/requirements.txt b/requirements.txt
index 83b5aee5..2d055041 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,3 +5,4 @@ scikit-learn!=0.23.0
 lmfit
 tqdm
 efel
+sbi
diff --git a/setup.py b/setup.py
index bd9d745e..ba3db004 100644
--- a/setup.py
+++ b/setup.py
@@ -1,27 +1,22 @@
 #! /usr/bin/env python
-'''
-brian2modelfitting setup script
-'''
+"""brian2modelfitting setup script"""
+
 import os
-import sys
 
 from setuptools import setup, find_packages
 
-def readme():
-    with open('README.md') as f:
-        return f.read()
 
 version = {}
 with open(os.path.join('brian2modelfitting', 'version.py')) as fp:
     exec(fp.read(), version)
 
-# Note that this does not set a version number explicitly, but automatically
-# figures out a version based on git tags
+with open('README.md') as f:
+    long_description = f.read()
+
 setup(name='brian2modelfitting',
       url='https://github.com/brian-team/brian2modelfitting',
       version=version['version'],
       packages=find_packages(),
-      # package_data={},
       install_requires=['numpy',
                         'brian2>=2.2',
                         'setuptools',
@@ -32,25 +27,26 @@ def readme():
       provides=['brian2modelfitting'],
       extras_require={'test': ['pytest'],
                       'docs': ['sphinx>=1.8'],
-                      'full': ['efel',
-                               'lmfit']},
+                      'full': ['efel', 'lmfit', 'sbi>=0.16.0'],
+                      },
       python_requires='>=3.6',
       use_2to3=False,
       zip_safe=False,
       description='Modelfitting Toolbox for the Brian 2 simulator',
-      long_description=readme(),
+      long_description=long_description,
       author='Aleksandra Teska, Marcel Stimberg, Romain Brette, Dan Goodman',
       author_email='team@briansimulator.org',
       license='CeCILL-2.1',
       classifiers=[
           'Development Status :: 4 - Beta',
           'Intended Audience :: Science/Research',
-          'License :: OSI Approved :: CEA CNRS Inria Logiciel Libre License, version 2.1 (CeCILL-2.1)',
+          'License :: OSI Approved :: CEA CNRS Inria Logiciel Libre License, '
+          'version 2.1 (CeCILL-2.1)',
           'Natural Language :: English',
           'Operating System :: OS Independent',
           'Programming Language :: Python',
           'Programming Language :: Python :: 3',
           'Topic :: Scientific/Engineering :: Bio-Informatics'
-      ],
-      keywords='model fitting computational neuroscience'
+          ],
+      keywords='model fitting computational neuroscience',
       )


Test patch:
diff --git a/brian2modelfitting/tests/test_inferencer.py b/brian2modelfitting/tests/test_inferencer.py
new file mode 100644
index 00000000..819f8a29
--- /dev/null
+++ b/brian2modelfitting/tests/test_inferencer.py
@@ -0,0 +1,351 @@
+from brian2 import nS, siemens, mV, volt, ms, TimedArray, Equations
+from brian2.devices.device import reinit_devices
+from brian2modelfitting.inferencer import (Inferencer,
+                                           get_param_dict,
+                                           calc_prior)
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal
+import pytest
+
+
+def test_get_param_dict():
+    param_names = ['a', 'b']
+    n_samples = 3
+    n_traces = 3
+    n_neurons = n_samples * n_traces
+
+    # check theta
+    theta = np.array([[1, 2],
+                      [3, 4],
+                      [5, 6]])
+    _theta = np.array([[1, 2],
+                       [1, 2],
+                       [1, 2],
+                       [3, 4],
+                       [3, 4],
+                       [3, 4],
+                       [5, 6],
+                       [5, 6],
+                       [5, 6]])
+    assert_equal(_theta, np.repeat(theta, repeats=n_traces, axis=0))
+
+    # check dictionary type and values
+    theta = np.repeat(theta, repeats=n_traces, axis=0)
+    d_param = {'a': np.array([1., 1., 1., 3., 3., 3., 5., 5., 5.]),
+               'b': np.array([2., 2., 2., 4., 4., 4., 6., 6., 6.])}
+    d_param_ret = get_param_dict(param_names=param_names,
+                                 param_values=_theta,
+                                 n_values=n_neurons)
+    assert isinstance(d_param_ret, dict)
+    assert_equal(d_param, d_param_ret)
+
+
+def test_calc_prior():
+    n_samples = 10
+    param_names = ['a', 'b']
+    a_lower_bound = 1
+    a_upper_bound = 10
+    a_mean = (a_lower_bound + a_upper_bound) / 2.
+    b_lower_bound = 0
+    b_upper_bound = 1
+    b_mean = (b_lower_bound + b_upper_bound) / 2.
+    prior_mean = np.array([a_mean, b_mean], dtype=np.float32)
+    a_std = (a_upper_bound - a_lower_bound) / np.sqrt(12)
+    b_std = (b_upper_bound - b_lower_bound) / np.sqrt(12)
+    prior_std = np.array([a_std, b_std], dtype=np.float32)
+
+    prior_sbi = calc_prior(param_names=param_names,
+                           a=[a_lower_bound*volt, a_upper_bound*volt],
+                           b=[b_lower_bound*volt, b_upper_bound*volt])
+
+    # check mean value
+    prior_sbi_mean = prior_sbi.mean.numpy()
+    assert_equal(prior_mean, prior_sbi_mean)
+
+    # check standard deviation
+    prior_sbi_std = prior_sbi.stddev.numpy()
+    assert_equal(prior_std, prior_sbi_std)
+
+    # check the type and shape of the prior
+    prior_sbi_sampled = prior_sbi.sample((n_samples, )).numpy()
+    assert isinstance(prior_sbi_sampled[0, 0], np.float32)
+    assert_equal(prior_sbi_sampled.shape, (n_samples, len(param_names)))
+
+
+E = 40*mV
+input_traces = np.zeros((10, 5))*volt
+for i in range(5):
+    input_traces[5:, i] = i * 10*mV
+input_traces = input_traces.T
+output_traces = input_traces * 10*nS
+output_traces = output_traces.T
+model = '''
+    I = g * (v - E) : amp
+    g : siemens (constant)
+    '''
+model_full = '''
+    I = g * (v - 40*mV) : amp
+    g : siemens (constant)
+    '''
+
+
+@pytest.fixture
+def setup(request):
+    dt = 0.01 * ms
+    inferencer = Inferencer(dt=dt,
+                            model=model,
+                            input={'v': input_traces},
+                            output={'I': output_traces})
+
+    def fin():
+        reinit_devices()
+    request.addfinalizer(fin)
+
+    return dt, inferencer
+
+
+@pytest.fixture
+def setup_features(request):
+    dt = 0.01 * ms
+    inferencer = Inferencer(dt=dt,
+                            model=model,
+                            input={'v': input_traces},
+                            output={'I': output_traces},
+                            features={'I': [lambda x: x.mean()]})
+
+    def fin():
+        reinit_devices()
+    request.addfinalizer(fin)
+
+    return dt, inferencer
+
+
+@pytest.fixture
+def setup_full(request):
+    dt = 0.01 * ms
+    inferencer = Inferencer(dt=dt,
+                            model=model_full,
+                            input={'v': input_traces},
+                            output={'I': output_traces},
+                            features={'I': [lambda x: x.mean()]})
+
+    def fin():
+        reinit_devices()
+    request.addfinalizer(fin)
+
+    return dt, inferencer
+
+
+def test_inferencer_init(setup):
+    _, inferencer = setup
+    attr_inferencer = ['dt', 'output_var', 'output', 'param_names',
+                       'n_traces', 'sim_time', 'output_dim', 'model',
+                       'input_traces', 'param_init', 'method',
+                       'threshold', 'reset', 'refractory', 'x_o',
+                       'features', 'params', 'n_samples', 'samples',
+                       'posterior', 'theta', 'x', 'sbi_device']
+    for attr in attr_inferencer:
+        assert hasattr(inferencer, attr)
+
+    assert inferencer.params is None
+    assert inferencer.n_samples is None
+    assert inferencer.samples is None
+    assert inferencer.posterior is None
+    assert inferencer.theta is None
+    assert inferencer.x is None
+
+    assert isinstance(inferencer.input_traces, TimedArray)
+    assert isinstance(inferencer.model, Equations)
+
+    target_var = f'{inferencer.output_var[0]}_target'
+    assert target_var in inferencer.model
+    assert inferencer.model[target_var].dim is inferencer.output_dim[0]
+
+    assert inferencer.features is None
+
+    observations = []
+    for o in output_traces:
+        o = np.array(o)
+        observations.append(o.ravel().astype(np.float32))
+    x_o = np.concatenate(observations)
+    assert_almost_equal(inferencer.x_o, x_o)
+
+
+def test_inferencer_init_features(setup_features):
+    _, inferencer = setup_features
+    attr_inferencer = ['dt', 'output_var', 'output', 'param_names',
+                       'n_traces', 'sim_time', 'output_dim', 'model',
+                       'input_traces', 'param_init', 'method',
+                       'threshold', 'reset', 'refractory', 'x_o',
+                       'features', 'params', 'n_samples', 'samples',
+                       'posterior', 'theta', 'x', 'sbi_device']
+    for attr in attr_inferencer:
+        assert hasattr(inferencer, attr)
+
+    assert inferencer.params is None
+    assert inferencer.n_samples is None
+    assert inferencer.samples is None
+    assert inferencer.posterior is None
+    assert inferencer.theta is None
+    assert inferencer.x is None
+
+    assert isinstance(inferencer.input_traces, TimedArray)
+    assert isinstance(inferencer.model, Equations)
+
+    target_var = f'{inferencer.output_var[0]}_target'
+    assert target_var in inferencer.model
+    assert inferencer.model[target_var].dim is inferencer.output_dim[0]
+
+    assert callable(inferencer.features['I'][0])
+
+    observations = []
+    for o in output_traces:
+        o = np.array(o)
+        print(o)
+        observations.append(o.mean())
+    x_o = np.array(observations, dtype=np.float32)
+    print(x_o)
+    print(inferencer.x_o)
+    assert_almost_equal(inferencer.x_o, x_o)
+
+
+def test_init_prior(setup):
+    _, inferencer = setup
+    lower_bound = 1e-9
+    upper_bound = 100e-9
+    g_mean = np.array((lower_bound + upper_bound) / 2, dtype=np.float32)
+    g_std = np.array((upper_bound - lower_bound) / np.sqrt(12),
+                     dtype=np.float32)
+    g_std = g_std.astype(np.float32)
+    prior = inferencer.init_prior(g=[lower_bound*siemens, upper_bound*siemens])
+
+    # check mean value
+    prior_mean = prior.mean.numpy()
+    assert_almost_equal(g_mean, prior_mean)
+
+    # check standard deviation
+    prior_std = prior.stddev.numpy()
+    assert_almost_equal(g_std, prior_std)
+
+
+def test_generate_training_data(setup):
+    _, inferencer = setup
+    lower_bound = 1*nS
+    upper_bound = 100*nS
+    n_samples = 100
+    prior = inferencer.init_prior(g=[lower_bound, upper_bound])
+    theta = inferencer.generate_training_data(n_samples, prior)
+    assert theta.shape == (100, 1)
+
+
+def test_extract_summary_statistics(setup):
+    _, inferencer = setup
+    lower_bound = 1*nS
+    upper_bound = 100*nS
+    n_samples = 100
+    n_traces = 5
+    prior = inferencer.init_prior(g=[lower_bound, upper_bound])
+    theta = inferencer.generate_training_data(n_samples, prior)
+    x = inferencer.extract_summary_statistics(theta)
+    assert x.shape == (100, n_traces * 10)
+
+
+def test_extract_summary_statistics_errors(setup):
+    _, inferencer = setup
+    lower_bound = 1*nS
+    upper_bound = 100*nS
+    n_samples = 100
+    prior = inferencer.init_prior(g=[lower_bound, upper_bound])
+    theta = inferencer.generate_training_data(n_samples, prior)
+    with pytest.raises(Exception):
+        _ = inferencer.extract_summary_statistics(theta, level=2)
+
+
+def test_save_summary_statistics_errors(setup):
+    _, inferencer = setup
+    with pytest.raises(AttributeError):
+        inferencer.save_summary_statistics('summary_stats')
+
+    theta = np.array([[1, 2], [3, 4]])
+    x = np.arange(2).reshape(2, 1)
+    with pytest.raises(TypeError):
+        inferencer.save_summary_statistics(1, theta=theta, x=x)
+        inferencer.save_summary_statistics(1, theta='theta', x=x)
+        inferencer.save_summary_statistics(1, theta=theta, x='x')
+        inferencer.save_summary_statistics(1, theta='theta', x='x')
+
+
+def test_load_summary_statistics(setup):
+    _, inferencer = setup
+    theta = np.array([[1, 2], [3, 4]])
+    x = np.arange(2).reshape(2, 1)
+    inferencer.save_summary_statistics('summ_stats.npz', theta=theta, x=x)
+    theta_load, x_load = inferencer.load_summary_statistics('summ_stats.npz')
+    assert_equal(theta_load, theta)
+    assert_equal(x_load, x)
+
+
+def test_init_inference_errors(setup):
+    _, inferencer = setup
+    prior = inferencer.init_prior(g=[1*nS, 100*nS])
+    with pytest.raises(NameError):
+        _ = inferencer.init_inference(inference_method='test',
+                                      density_estimator_model='mdn',
+                                      prior=prior,
+                                      sbi_device='cpu')
+        _ = inferencer.init_inference(inference_method='SNPE',
+                                      density_estimator_model='test',
+                                      prior=prior,
+                                      sbi_device='cpu')
+
+    # sbi_device should be automatically set to 'cpu'
+    _ = inferencer.init_inference(inference_method='SNPE',
+                                  density_estimator_model='mdn',
+                                  prior=prior,
+                                  sbi_device='test')
+    assert inferencer.sbi_device == 'cpu'
+
+
+def test_infer_step(setup_full):
+    _, inferencer = setup_full
+    prior = inferencer.init_prior(g=[1*nS, 100*nS])
+    inference = inferencer.init_inference(inference_method='SNPE',
+                                          density_estimator_model='mdn',
+                                          prior=prior,
+                                          sbi_device='cpu')
+    posterior = inferencer.infer_step(proposal=prior,
+                                      n_samples=10,
+                                      inference=inference)
+    assert posterior._method_family == 'snpe'
+    assert_equal(np.array(posterior._x_shape), np.array([1, 5]))
+
+
+def test_infer_step_errors(setup_full):
+    _, inferencer = setup_full
+    prior = inferencer.init_prior(g=[1*nS, 100*nS])
+    inference = inferencer.init_inference(inference_method='SNPE',
+                                          density_estimator_model='mdn',
+                                          prior=prior,
+                                          sbi_device='cpu')
+    with pytest.raises(ValueError):
+        _ = inferencer.infer_step(proposal=prior,
+                                  inference=inference)
+
+
+def test_infer(setup):
+    _, inferencer = setup
+    posterior = inferencer.infer(n_samples=10, g=[1*nS, 100*nS])
+    assert posterior == inferencer.posterior
+
+    posterior_multi_round = inferencer.infer(n_samples=10, n_rounds=2,
+                                             restart=True, g=[1*nS, 100*nS])
+    assert posterior_multi_round != posterior
+    assert posterior_multi_round == inferencer.posterior
+
+
+def test_load_posterior(setup):
+    _, inferencer = setup
+    _ = inferencer.infer(n_samples=10, g=[1*nS, 100*nS])
+    inferencer.save_posterior('posterior.pth')
+    posterior_load = inferencer.load_posterior('posterior.pth')
+    assert inferencer.posterior.net == posterior_load.net