Source code for elfi.methods.post_processing

"""Post-processing for posterior samples from other ABC algorithms.

References
----------
Fundamentals and Recent Developments in Approximate Bayesian Computation
Lintusaari et. al
Syst Biol (2017) 66 (1): e66-e82.
https://doi.org/10.1093/sysbio/syw077

"""
import warnings

import numpy as np
from sklearn.linear_model import LinearRegression

from . import results

__all__ = ('LinearAdjustment', 'adjust_posterior')


class RegressionAdjustment(object):
    """Base class for regression adjustments.

    Each parameter is assumed to be a scalar. A local regression is
    fitted for each parameter individually using the values of the
    summary statistics as the regressors.  The regression model can be
    any object implementing a `fit()` method. All keyword arguments
    given to the constructor are passed to the regression model.

    Subclasses need to implement the methods `_adjust` and
    `_input_variables`.  They must also specify the class variables
    `_regression_model` and `_name`.  See the individual documentation
    and the `LinearAdjustment` class for further detail.

    Parameters
    ----------
    kwargs**
      keyword arguments to pass to the regression model

    Attributes
    ----------
    regression_models
      a list of fitted regression model instances
    parameter_names
      a list of parameter names
    sample
      the sample object from an ABC algorithm
    X
      the regressors for the regression model

    """

    _regression_model = None
    _name = 'RegressionAdjustment'

    def __init__(self, **kwargs):
        self._model_kwargs = kwargs
        self._fitted = False
        self.regression_models = []
        self._X = None
        self._sample = None
        self._parameter_names = None
        self._finite = []

    @property
    def parameter_names(self):
        self._check_fitted()
        return self._parameter_names

    @property
    def sample(self):
        self._check_fitted()
        return self._sample

    @property
    def X(self):
        self._check_fitted()
        return self._X

    def _check_fitted(self):
        if not self._fitted:
            raise ValueError("The regression model must be fitted first. " "Use the fit() method.")

    def fit(self, sample, model, summary_names, parameter_names=None):
        """Fit a regression adjustment model to the posterior sample.

        Non-finite values in the summary statistics and parameters
        will be omitted.

        Parameters
        ----------
        sample : elfi.methods.Sample
          a sample object from an ABC method
        model : elfi.ElfiModel
          the inference model
        summary_names : list[str]
          a list of names for the summary nodes
        parameter_names : list[str] (optional)
          a list of parameter names

        """
        self._X = self._input_variables(model, sample, summary_names)
        self._sample = sample
        self._parameter_names = parameter_names or sample.parameter_names
        self._get_finite()

        for pair in self._pairs():
            self.regression_models.append(self._fit1(*pair))

        self._fitted = True

    def _fit1(self, X, y):
        return self._regression_model(**self._model_kwargs).fit(X, y)

    def _pairs(self):
        # TODO: Access the variables through the getters
        for (i, name) in enumerate(self._parameter_names):
            X = self._X[self._finite[i], :]
            p = self._sample.outputs[name][self._finite[i]]
            yield X, p

    def adjust(self):
        """Adjust the posterior.

        Only the non-finite values used to fit the regression model
        will be adjusted.

        Returns
        -------
          a Sample object containing the adjusted posterior

        """
        outputs = {}
        for (i, name) in enumerate(self.parameter_names):
            theta_i = self.sample.outputs[name][self._finite[i]]
            adjusted = self._adjust(i, theta_i, self.regression_models[i])
            outputs[name] = adjusted

        res = results.Sample(
            method_name=self._name, outputs=outputs, parameter_names=self._parameter_names)
        return res

    def _adjust(self, i, theta_i, regression_model):
        """Adjust a single parameter using a fitted regression model.

        Parameters
        ----------
        i : int
          the index of the parameter
        theta_i : np.ndarray
          a vector of parameter values to adjust
        regression_model
          a fitted regression model

        Returns
        -------
        adjusted_theta_i : np.ndarray
          an adjusted version of the parameter values

        """
        raise NotImplementedError

    def _input_variables(self, model, sample, summary_names):
        """Construct a matrix of regressors.

        Parameters
        ----------
        model : elfi.ElfiModel
          the inference model
        sample
          a sample object from an ABC algorithm
        summary_names : list[str]
          names of the summary nodes

        Returns
        -------
        X
          a numpy array of regressors

        """
        raise NotImplementedError

    def _get_finite(self):
        # TODO: Access the variables through the getters
        finite_inputs = np.isfinite(self._X).all(axis=1)
        finite = [
            finite_inputs & np.isfinite(self._sample.outputs[p]) for p in self._parameter_names
        ]
        all_finite = all(map(all, finite))
        self._finite = finite
        if not (all(finite_inputs) and all_finite):
            warnings.warn("Non-finite inputs and outputs will be omitted.")


[docs]class LinearAdjustment(RegressionAdjustment): """Regression adjustment using a local linear model.""" _regression_model = LinearRegression _name = 'LinearAdjustment' def _adjust(self, i, theta_i, regression_model): b = regression_model.coef_ return theta_i - self.X[self._finite[i], :].dot(b) def _input_variables(self, model, sample, summary_names): """Regress on the differences to the observed summaries.""" observed_summaries = np.stack([model[s].observed for s in summary_names], axis=1) summaries = np.stack([sample.outputs[name] for name in summary_names], axis=1) return summaries - observed_summaries
[docs]def adjust_posterior(sample, model, summary_names, parameter_names=None, adjustment='linear'): """Adjust the posterior using local regression. Note that the summary nodes need to be explicitly included to the sample object with the `output_names` keyword argument when performing the inference. Parameters ---------- sample : elfi.methods.results.Sample a sample object from an ABC algorithm model : elfi.ElfiModel the inference model summary_names : list[str] names of the summary nodes parameter_names : list[str] (optional) names of the parameters adjustment : RegressionAdjustment or string a regression adjustment object or a string specification Accepted values for the string specification: - 'linear' Returns ------- elfi.methods.results.Sample a Sample object with the adjusted posterior Examples -------- >>> import elfi >>> from elfi.examples import gauss >>> m = gauss.get_model() >>> res = elfi.Rejection(m['d'], output_names=['ss_mean', 'ss_var'], ... batch_size=10).sample(500, bar=False) >>> adj = adjust_posterior(res, m, ['ss_mean', 'ss_var'], ['mu'], LinearAdjustment()) """ adjustment = _get_adjustment(adjustment) adjustment.fit( model=model, sample=sample, parameter_names=parameter_names, summary_names=summary_names) return adjustment.adjust()
def _get_adjustment(adjustment): adjustments = {'linear': LinearAdjustment} if isinstance(adjustment, RegressionAdjustment): return adjustment elif isinstance(adjustment, str): try: return adjustments.get(adjustment, None)() except TypeError: raise ValueError("Could not find " "adjustment method:{}".format(adjustment))