diff --git a/calibr8/__init__.py b/calibr8/__init__.py index af6a14a..298e1df 100644 --- a/calibr8/__init__.py +++ b/calibr8/__init__.py @@ -1,3 +1,10 @@ + +from .contrib.noise import ( + NormalNoise, + LaplaceNoise, + LogNormalNoise, + StudentTNoise, +) from .contrib.base import ( BaseAsymmetricLogisticT, BaseLogIndependentAsymmetricLogisticT, @@ -6,6 +13,7 @@ ) from .core import ( CalibrationModel, + DistributionMixin, InferenceResult, UnivariateInferenceResult, NumericPosterior, diff --git a/calibr8/contrib/base.py b/calibr8/contrib/base.py index 4f41501..da8518f 100644 --- a/calibr8/contrib/base.py +++ b/calibr8/contrib/base.py @@ -2,111 +2,14 @@ This module implements generic, reusable calibration models that can be subclassed to implement custom calibration models. """ -import numpy -import scipy import typing -import warnings +from . import noise from .. import core -from .. import utils -try: - try: - import pymc3 as pm - except ModuleNotFoundError: - import pymc as pm -except ModuleNotFoundError: - pm = utils.ImportWarner('pymc3') - -class BaseModelT(core.CalibrationModel): - def loglikelihood(self, *, y, x, name: str=None, replicate_id: str=None, dependent_key: str=None, theta=None, **dist_kwargs): - """ Loglikelihood of observation (dependent variable) given the independent variable. - - If both x and y are 1D-vectors, they must have the same length and the likelihood will be evaluated elementwise. - - For a 2-dimensional `x`, the implementation *should* broadcast and return a result that has - the same length as the first dimension of `x`. - - Parameters - ---------- - y : scalar or array-like - observed measurements (dependent variable) - x : scalar, array-like or TensorVariable - assumed independent variable - name : str - Name for the likelihood variable in a PyMC3 model (tensor mode). - Previously this was `f'{replicate_id}.{dependent_key}'`. - replicate_id : optional, str - Deprecated; pass the `name` kwarg instead. - dependent_key : optional, str - Deprecated; pass the `name` kwarg instead. - theta : optional, array-like - Parameters for the calibration model to use instead of `theta_fitted`. - The vector must have the correct length, but can have numeric and or symbolic entries. - Use this kwarg to run MCMC on calibration model parameters. - **dist_kwargs : dict - Additional keyword arguments are forwarded to the `pm.StudentT` distribution. - Most prominent example: `dims`. - - Returns - ------- - L : float or TensorVariable - sum of log-likelihoods - """ - if theta is None: - if self.theta_fitted is None: - raise Exception('No parameter vector was provided and the model is not fitted with data yet.') - theta = self.theta_fitted - - if not isinstance(x, (list, numpy.ndarray, float, int)) and not utils.istensor(x): - raise ValueError( - f'Input x must be a scalar, TensorVariable or an array-like object, but not {type(x)}' - ) - if not isinstance(y, (list, numpy.ndarray, float, int)) and not utils.istensor(x): - raise ValueError( - f'Input y must be a scalar or an array-like object, but not {type(y)}' - ) - - mu, scale, df = self.predict_dependent(x, theta=theta) - if utils.istensor(x) or utils.istensor(theta): - if pm.Model.get_context(error_if_none=False) is not None: - if replicate_id and dependent_key: - warnings.warn( - "The `replicate_id` and `dependent_key` parameters are deprecated. Use `name` instead.", - DeprecationWarning - ) - name = f'{replicate_id}.{dependent_key}' - if not name: - raise ValueError("A `name` must be specified for the PyMC likelihood.") - rv = pm.StudentT( - name, - mu=mu, - sigma=scale, - nu=df, - observed=y, - **dist_kwargs or {} - ) - else: - rv = pm.StudentT.dist( - mu=mu, - sigma=scale, - nu=df, - **dist_kwargs or {} - ) - # The API to get log-likelihood tensors differs between PyMC versions - if pm.__version__[0] == "3": - if isinstance(rv, pm.model.ObservedRV): - return rv.logpt.sum() - elif isinstance(rv, pm.Distribution): - return rv.logp(y).sum() - else: - return pm.logpt(rv, y, sum=True) - else: - # If `x` is given as a column vector, this model can broadcast automatically. - # This gives considerable performance benefits for the `likelihood(..., scan_x=True)` - # case which is relevant for `infer_independent`. - return numpy.sum(scipy.stats.t.logpdf(x=y, loc=mu, scale=scale, df=df), axis=-1) +class BaseModelT(core.CalibrationModel, noise.StudentTNoise): + pass class BasePolynomialModelT(BaseModelT): diff --git a/calibr8/contrib/noise.py b/calibr8/contrib/noise.py new file mode 100644 index 0000000..9df722a --- /dev/null +++ b/calibr8/contrib/noise.py @@ -0,0 +1,64 @@ +import numpy +import scipy.stats + +from .. core import DistributionMixin +from .. utils import HAS_PYMC, pm + + +class NormalNoise(DistributionMixin): + """Normal noise, predicted in terms of mean and standard deviation.""" + scipy_dist = scipy.stats.norm + pymc_dist = pm.Normal if HAS_PYMC else None + + @staticmethod + def to_scipy(*params): + return dict(loc=params[0], scale=params[1]) + + @staticmethod + def to_pymc(*params): + return dict(mu=params[0], sigma=params[1]) + + +class LaplaceNoise(DistributionMixin): + """Normal noise, predicted in terms of mean and scale.""" + scipy_dist = scipy.stats.laplace + pymc_dist = pm.Laplace if HAS_PYMC else None + + @staticmethod + def to_scipy(*params): + return dict(loc=params[0], scale=params[1]) + + @staticmethod + def to_pymc(*params): + return dict(mu=params[0], b=params[1]) + + +class LogNormalNoise(DistributionMixin): + """Log-Normal noise, predicted in logarithmic mean and standard deviation. + ⚠ This corresponds to the NumPy/Aesara/PyMC parametrization! + """ + scipy_dist = scipy.stats.lognorm + pymc_dist = pm.Lognormal if HAS_PYMC else None + + @staticmethod + def to_scipy(*params): + # SciPy wants linear scale mean and log scale standard deviation! + return dict(scale=numpy.exp(params[0]), s=params[1]) + + @staticmethod + def to_pymc(*params): + return dict(mu=params[0], sigma=params[1]) + + +class StudentTNoise(DistributionMixin): + """Student-t noise, predicted in terms of mean, scale and degree of freedom.""" + scipy_dist = scipy.stats.t + pymc_dist = pm.StudentT if HAS_PYMC else None + + @staticmethod + def to_scipy(*params): + return dict(loc=params[0], scale=params[1], df=params[2]) + + @staticmethod + def to_pymc(*params): + return dict(mu=params[0], sigma=params[1], nu=params[2]) diff --git a/calibr8/core.py b/calibr8/core.py index 0b4d7ad..f9707a1 100644 --- a/calibr8/core.py +++ b/calibr8/core.py @@ -16,6 +16,7 @@ from typing import Union from . import utils +from . utils import pm __version__ = '6.1.3' @@ -369,7 +370,34 @@ def _infer_univariate_independent( ) -class CalibrationModel: +class DistributionMixin: + """Maps the values returned by `CalibrationModel.predict_dependent` + to a SciPy distribution and its parameters, and optionally also to + a PyMC distribution and its parameters. + """ + scipy_dist = None + pymc_dist = None + + def to_scipy(*args) -> dict: + raise NotImplementedError("This model does not implement a mapping to SciPy distribution parameters.") + + def to_pymc(*args) -> dict: + raise NotImplementedError("This model does not implement a mapping to PyMC distribution parameters.") + + +def _inherits_noisemodel(cls): + """Determines if cls is a sub-type of DistributionMixin that's not DistributionMixin or a calibration model.""" + for m in cls.__mro__: + if ( + issubclass(m, DistributionMixin) + and m is not DistributionMixin + and not issubclass(m, CalibrationModel) + ): + return True + return False + + +class CalibrationModel(DistributionMixin): """A parent class providing the general structure of a calibration model.""" def __init__(self, independent_key:str, dependent_key:str, *, theta_names:typing.Tuple[str], ndim=1): @@ -388,6 +416,14 @@ def __init__(self, independent_key:str, dependent_key:str, *, theta_names:typing Most calibrations are univariate (ndim=1). Multivariate calibrations (ndim > 1) may need to override the `infer_independent` implementation. """ + if not _inherits_noisemodel(type(self)): + warnings.warn( + "This model does not implement a noise model yet." + "\nAdd a noise model mixin to your class definition. For example:" + "\n`class MyModel(CalibrationModel, LaplaceNoise)`", + DeprecationWarning, + stacklevel=2 + ) # make sure that the inheriting type has no required constructor (kw)args args, varargs, varkw, defaults, kwonlyargs, kwonlydefaults, annotations = inspect.getfullargspec(type(self).__init__) n_defaults = 0 if not defaults else len(defaults) @@ -523,7 +559,17 @@ def infer_independent( ci_prob=ci_prob, ) - def loglikelihood(self, *, y, x, theta=None): + def loglikelihood( + self, + *, + y, + x, + name: str=None, + replicate_id: str=None, + dependent_key: str=None, + theta=None, + **dist_kwargs + ): """ Loglikelihood of observation (dependent variable) given the independent variable. If both x and y are 1D-vectors, they must have the same length and the likelihood will be evaluated elementwise. @@ -537,15 +583,72 @@ def loglikelihood(self, *, y, x, theta=None): observed measurements (dependent variable) x : scalar, array-like or TensorVariable assumed independent variable + name : str + Name for the likelihood variable in a PyMC model (tensor mode). + Previously this was `f'{replicate_id}.{dependent_key}'`. + replicate_id : optional, str + Deprecated; pass the `name` kwarg instead. + dependent_key : optional, str + Deprecated; pass the `name` kwarg instead. theta : optional, array-like - model parameters + Parameters for the calibration model to use instead of `theta_fitted`. + The vector must have the correct length, but can have numeric and or symbolic entries. + Use this kwarg to run MCMC on calibration model parameters. + **dist_kwargs : dict + Additional keyword arguments are forwarded to the PyMC distribution. + Most prominent example: `dims`. Returns ------- L : float or TensorVariable sum of log-likelihoods """ - raise NotImplementedError('The loglikelihood function should be implemented by the inheriting class.') + if theta is None: + if self.theta_fitted is None: + raise Exception('No parameter vector was provided and the model is not fitted with data yet.') + theta = self.theta_fitted + + if not isinstance(x, (list, numpy.ndarray, float, int)) and not utils.istensor(x): + raise ValueError( + f'Input x must be a scalar, TensorVariable or an array-like object, but not {type(x)}' + ) + if not isinstance(y, (list, numpy.ndarray, float, int)) and not utils.istensor(x): + raise ValueError( + f'Input y must be a scalar or an array-like object, but not {type(y)}' + ) + + params = self.predict_dependent(x, theta=theta) + if utils.istensor(x) or utils.istensor(theta): + if pm.Model.get_context(error_if_none=False) is not None: + if replicate_id and dependent_key: + warnings.warn( + "The `replicate_id` and `dependent_key` parameters are deprecated. Use `name` instead.", + DeprecationWarning + ) + name = f'{replicate_id}.{dependent_key}' + if not name: + raise ValueError("A `name` must be specified for the PyMC likelihood.") + rv = self.pymc_dist( + name, + **self.to_pymc(*params), + observed=y, + **dist_kwargs or {} + ) + else: + rv = self.pymc_dist.dist( + **self.to_pymc(*params), + **dist_kwargs or {} + ) + # The API to get log-likelihood tensors differs between PyMC versions + if pm.__version__[0] == "3": + if isinstance(rv, pm.model.ObservedRV): + return rv.logpt.sum() + elif isinstance(rv, pm.Distribution): + return rv.logp(y).sum() + else: + return pm.logpt(rv, y, sum=True) + else: + return self.scipy_dist.logpdf(x=y, **self.to_scipy(*params)).sum(axis=-1) def likelihood(self, *, y, x, theta=None, scan_x: bool=False): """ Likelihood of observation (dependent variable) given the independent variable. @@ -572,8 +675,8 @@ def likelihood(self, *, y, x, theta=None, scan_x: bool=False): try: # Try to pass `x` as a column vector to benefit from broadcasting # if that's implemented by the underlying model. - result = numpy.exp(self.loglikelihood(y=y, x=x[:, None], theta=theta)) - if not result.shape == x.shape: + result = numpy.exp(self.loglikelihood(y=y, x=x[..., None], theta=theta)) + if not numpy.shape(result) == numpy.shape(x): raise ValueError("The underlying model does not seem to implement broadcasting.") return result except: diff --git a/calibr8/tests.py b/calibr8/tests.py index d27cd69..f500522 100644 --- a/calibr8/tests.py +++ b/calibr8/tests.py @@ -37,13 +37,25 @@ dir_testfiles = pathlib.Path(pathlib.Path(__file__).absolute().parent, 'testfiles') -class _TestModel(calibr8.CalibrationModel): +class _TestModel(calibr8.CalibrationModel, calibr8.NormalNoise): def __init__(self, independent_key=None, dependent_key=None, theta_names=None): if theta_names is None: theta_names = tuple('a,b,c'.split(',')) super().__init__(independent_key='I', dependent_key='D', theta_names=theta_names) +class _TestBivariateLinearModel(calibr8.CalibrationModel, calibr8.NormalNoise): + def __init__(self, independent_key: str=None, dependent_key: str=None, theta_names=None): + super().__init__(independent_key="x1,x2", dependent_key="y", theta_names="i,s1,s2,sd".split(","), ndim=2) + + def predict_dependent(self, x, *, theta=None): + x = numpy.array(x) + if theta is None: + theta = self.theta_fitted + i, s1, s2, sd = theta + return i + s1 * x[..., 0] + s2 * x[..., 1], sd + + class _TestPolynomialModel(calibr8.BasePolynomialModelT): def __init__(self, independent_key=None, dependent_key=None, theta_names=None, *, scale_degree=0, mu_degree=1): super().__init__(independent_key='I', dependent_key='D', mu_degree=mu_degree, scale_degree=scale_degree) @@ -85,12 +97,50 @@ def test_warnings(self): pass -class TestBasicCalibrationModel: +class TestNoiseModels: + def test_base_class(self): + assert calibr8.DistributionMixin.scipy_dist is None + assert calibr8.DistributionMixin.pymc_dist is None + with pytest.raises(NotImplementedError, match="mapping to SciPy"): + assert calibr8.DistributionMixin.to_scipy() + with pytest.raises(NotImplementedError, match="mapping to PyMC"): + assert calibr8.DistributionMixin.to_pymc() + pass + + @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') + @pytest.mark.parametrize( + "cls,a,b,params", [ + (calibr8.NormalNoise, -1, 1, (0.2, 1.2)), + (calibr8.LaplaceNoise, -1, 1, (-0.2, 0.8)), + (calibr8.LogNormalNoise, 0.1, 2, (-0.1, 0.2)), + (calibr8.StudentTNoise, -2, 10, (1.5, 0.3, 4)), + ] + ) + def test_parametrization_equivalence(self, cls: calibr8.DistributionMixin, a, b, params): + """Validates that SciPy and PyMC parametrization give the same distribution.""" + x = numpy.linspace(a, b, 7) + kwargs_scipy = cls.to_scipy(*params) + kwargs_pymc = cls.to_pymc(*params) + + result_scipy = cls.scipy_dist.logpdf(x=x, **kwargs_scipy) + + rv = cls.pymc_dist.dist(**kwargs_pymc) + if not hasattr(pm, "logp"): + # PyMC v3 syntax + result_pymc = rv.logp(x).eval() + else: + result_pymc = [pm.logp(rv, xi).eval() for xi in x] + + # The resulting log-PDF evaluations should be really close + numpy.testing.assert_allclose(result_scipy, result_pymc) + pass + + +class TestBaseCalibrationModel: def test_init(self): em = _TestModel('I', 'D', theta_names=tuple('c,d,e'.split(','))) assert em.independent_key == 'I' assert em.dependent_key == 'D' - print(em.theta_names) assert em.theta_names == ('c', 'd', 'e') assert em.theta_bounds is None assert em.theta_guess is None @@ -99,6 +149,31 @@ def test_init(self): assert em.cal_independent is None assert em.cal_dependent is None pass + + def test_init_requires_noisemodel(self): + assert issubclass(calibr8.CalibrationModel, calibr8.DistributionMixin) + + class _InvalidModel(calibr8.CalibrationModel): + def __init__(self): + super().__init__("A", "B", theta_names=tuple('abc'), ndim=1) + + assert not calibr8.core._inherits_noisemodel(calibr8.DistributionMixin) + assert not calibr8.core._inherits_noisemodel(calibr8.CalibrationModel) + assert calibr8.core._inherits_noisemodel(calibr8.LaplaceNoise) + assert not calibr8.core._inherits_noisemodel(_InvalidModel) + assert calibr8.core._inherits_noisemodel(_TestModel) + + # Check that a warning is raised when the model class does not inherit a DistributionMixin _sub_class + with pytest.warns(DeprecationWarning, match="does not implement a noise model"): + _InvalidModel() + + # Check that no warning is raised when the inheritance is correct + with pytest.warns(None) as record: + _TestModel() + assert len(record) == 0 + + assert issubclass(_TestModel, calibr8.DistributionMixin) + pass def test_constructor_signature_check(self): class EM_OK(calibr8.CalibrationModel): @@ -109,13 +184,13 @@ def __init__(self, arg1=1, *, kwonly=2, kwonlydefault=4): class EM_args(calibr8.CalibrationModel): def __init__(self, arg1): super().__init__('I', 'D', theta_names=tuple('abc')) - with pytest.raises(TypeError): + with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): EM_args(arg1=3) class EM_kwargs(calibr8.CalibrationModel): def __init__(self, *, kwonly, kwonlydefault=4): super().__init__('I', 'D', theta_names=tuple('abc')) - with pytest.raises(TypeError): + with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): EM_kwargs(kwonly=3) pass @@ -128,7 +203,7 @@ def test_exceptions(self): cmodel.predict_dependent(x) with pytest.raises(NotImplementedError, match="predict_independent function"): cmodel.predict_independent(x) - with pytest.raises(NotImplementedError, match="loglikelihood function"): + with pytest.raises(NotImplementedError, match="predict_dependent function"): cmodel.loglikelihood(y=y, x=x, theta=[1, 2, 3]) with pytest.raises(ValueError, match=r"Unexpected `ci_prob`"): cmodel.loglikelihood = lambda x, y, theta: 1 @@ -223,6 +298,44 @@ def test_objective(self): assert -obj_max(theta) == obj_min(theta) pass + def test_likelihood_multivariate(self): + cm = _TestBivariateLinearModel() + cm.theta_fitted = (0.5, 1, 2, 0.4) + + X = [ + [1, 2], + [0.5, 1], + [1.5, 2.5], + ] + Y = [5.5, 3.0, 7.0] + y_obs = [5.6, 3.0, 6.9] + + # Check the expected output for one coordinate + assert cm.predict_dependent(X[0]) == (Y[0], 0.4) + + # And its broadcasting for multiple coordinates + numpy.testing.assert_array_equal(cm.predict_dependent(X)[0], Y) + + # Expected likelihoods of coordinate/observation paris + LL_elemwise = cm.scipy_dist.logpdf(x=y_obs, loc=Y, scale=0.4).sum() + + # Expected likelihoods of all observation at each coordinate + LL_scan = numpy.array([ + cm.scipy_dist.logpdf(x=y_obs, loc=y, scale=0.4).sum() + for y in Y + ]) + + assert LL_elemwise.shape == () + assert LL_scan.shape == (3,) + + # Testing the underlying loglikelihood + numpy.testing.assert_array_equal(cm.loglikelihood(y=y_obs, x=X), LL_elemwise) + + # Now check if likelihood wraps it correctly + numpy.testing.assert_array_equal(cm.likelihood(y=y_obs, x=X, scan_x=False), numpy.exp(LL_elemwise)) + numpy.testing.assert_array_equal(cm.likelihood(y=y_obs, x=X, scan_x=True), numpy.exp(LL_scan)) + pass + class TestUnivariateInferenceHelpers: def test_get_eti(self): @@ -941,6 +1054,73 @@ def test_likelihood(self): assert x_dense[numpy.argmax(actual)] == 2 pass + @pytest.mark.xfail(reason="Draft test case, see https://github.com/JuBiotech/calibr8/issues/15") + def test_likelihood_nobroadcasting_fallback(self): + class _TestSwitchableBroadcastingModel(calibr8.CalibrationModel, calibr8.NormalNoise): + def __init__(self): + self.x_shapes = [] + self.can_broadcast = False + super().__init__(independent_key="I", dependent_key="D", theta_names="i,s,sd".split(","), ndim=1) + + def predict_dependent(self, x, *, theta=None): + if theta is None: + theta = self.theta_fitted + i, s, sd = self.theta_fitted + # This part broadcasts just fine + x = numpy.array(x) + mu = i + x * s + return mu, sd + + def loglikelihood(self, *, y, x, **kwargs): + # This overrides the native CalibrationModel.loglikelihood with one + # that can be externally set to return non-broadcasted results. + LL_broadcasted = super().loglikelihood(y=y, x=x, **kwargs) + if not self.can_broadcast: + return LL_broadcasted.sum() + return LL_broadcasted + + cm = _TestSwitchableBroadcastingModel() + cm.theta_fitted = [0.5, 0.6, 0.7] + + X = numpy.array([1, 2, 3]) + Y = [1.1, 1.7, 2.3] + y_obs = [1.0, 1.7, 2.4] + + # Check the prediction of the test model + for x, y in zip(X, Y): + assert cm.predict_dependent(x) == (y, 0.7) + + # The predict_dependent can broadcast + mu, sd = cm.predict_dependent(X) + numpy.testing.assert_array_equal(mu, Y) + assert sd == 0.7 + + # Test the switching between can_broadcast modes + cm.can_broadcast = True + LL = cm.loglikelihood(x=X[..., None], y=y_obs) + assert numpy.shape(LL) == (3,) + + cm.can_broadcast = False + LL = cm.loglikelihood(x=X[..., None], y=y_obs) + assert numpy.shape(LL) == () + + # The CalibrationModel.likelihood should give the same results either way. + cm.can_broadcast = True + L_broadcasted = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) + cm.can_broadcast = False + L_looped = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) + + numpy.testing.assert_array_equal(L_broadcasted, L_looped) + + # Of course the values should also be correct. + L_expected = numpy.exp([ + cm.scipy_dist.logpdf(x=y_obs, loc=mui, scale=sd).sum() + for mui in mu + ]) + assert numpy.shape(L_expected) == (3,) + numpy.testing.assert_array_equal(L_broadcasted, L_expected) + pass + class TestBaseAsymmetricLogisticModelT: def test_predict_dependent(self): diff --git a/docs/source/calibr8_inheritance.md b/docs/source/calibr8_inheritance.md new file mode 100644 index 0000000..3b64027 --- /dev/null +++ b/docs/source/calibr8_inheritance.md @@ -0,0 +1,84 @@ +# Inheritance and composition of `calibr8` models +Models in `calibr8` are implemented as classes, inheriting from the `CalibrationModel` interface* either directly or via configureable `calibr8.BaseXYZ` model classes. +Generalization across different distributions in the noise model is provided with mixins**. + +The reasoning for this design choice is explained in the following sections. + +*[Further reading on informal interfaces in Python](https://realpython.com/python-interface/).
+**[Further reading on mixins](https://realpython.com/inheritance-composition-python/#mixing-features-with-mixin-classes).
+ +## Inheritance of model classes +The choice of the inheritance architecture is based on an analysis of two important aspects of a calibration model: +1. What kind of independent variable is the model working with? (continuous/discrete) +2. How many dimensions does the independent variable have? (univariate/multivariate) + +Note that the `(log)likelihood` of a model is agnostic to the dependent variable being continuous/discrete, and therefore we do not have to consider it here. + +The following tables summarizes the long-form matrix of the various combinations and corresponding implementations. +Note that since MCMC can be used for any of them, the "inference strategy" column denotes the computationally simplest generic strategy. + +| ndim | independent variable | integration/inference strategy | `InferenceResult` (properties) | model subclass | +|----------------|----------------------|--------------------|--------------------------------|----------------| +| 1 | continuous | `trapz` | `ContinuousUnivariateInference` (ETI, HDI, median) | `ContinuousUnivariateModel` | +| >1 | continuous | MCMC | `ContinuousMultivariateInference` (idata, ETI, HDI) | `ContinuousMultivariateModel` | +| 1 | discrete | summation & division of likelihoods | *, `DiscreteUnivariateInference` (probability vector, mode) | *, `DiscreteUnivariateModel` | +| >1 | discrete | summation & division of likelihoods | *, `DiscreteAnyvariateInference` (probability ndarray) | *, `DiscreteMultivariateModel` | +| >1 | **mix of continuous & discrete 🀯 | MCMC | *, `MixedMultivariateInference` (idata) | *, `MixedMultivariateModel` | + +*Not implemented.
+**Needs custom SciPy and PyMC distribution. + +## Composition of distribution mixins +The `CalibrationModel`, and specifically it's `loglikelihood` implementation can generalize across distributions (noise models) for the dependent variable based on a mapping of predicted distribution parameters (`predict_dependent`) to the corresponding SciPy, and optionally a PyMC distribution. + +This mapping is managed via class attributes and staticmethods and provided by subclassing a `DistributionMixin`. +For example, a user-implemented model may subclass the `calibr8.LogNormalNoise` and implement a `predict_dependent` method to return a tuple of the two distribution parameters. + +```python +class LogNormalNoise(DistributionMixin): + """Log-Normal noise, predicted in logarithmic mean and standard deviation. + ⚠ This corresponds to the NumPy/Aesara/PyMC parametrization! + """ + scipy_dist = scipy.stats.lognorm + pymc_dist = pm.Lognormal if HAS_PYMC else None + + @staticmethod + def to_scipy(*params): + # SciPy wants linear scale mean and log scale standard deviation! + return dict(scale=numpy.exp(params[0]), s=params[1]) + + @staticmethod + def to_pymc(*params): + return dict(mu=params[0], sigma=params[1]) +``` + +## Localization of implementations +Most of the source codeβ€”and certainly the most intricate codeβ€”of a `calibr8` model is located in the classes provided by the library. +For most applications the user may directly use a common model class from `calibr8`, or implement their own class with nothing more than the `__init__` method. + +If the desired model structure is not already provided by `calibr8.BaseXYZ` types, a custom `predict_dependent` can be implemented by subclassing from a model subclass. +Examples are the use of uncommon noise models, or multivariate calibration models. + +| name β•² class | `CalibrationModel` | model subclass | `CustomModel` | `BaseXYZ` | `UserModel` | +|--------------|--------------------|----------------|---------------|-----------|-------------| +| implemented by | `calibr8` | `calibr8` | user | `calibr8` | user +| inherits from | `object` | `CalibrationModel` | model subclass | model subclass | `BaseXYZ` | +| adds mixin | `DistributionMixin` | | user-specified | `NormalNoise` or `StudentTNoise` | | +| generalizes for | all models | ndim and type of
independent variable | | common model structures
(e.g. univariate polynomial) +| `__init__` | βœ” | βœ” | βœ” | (βœ”) | βœ” | +| `theta_*` | βœ” | ← | ← | ← | ← | +| `save` | βœ” | ← | ← | ← | ← | +| `load` | βœ” | ← | ← | ← | ← | +| `objective` | βœ” | ← | ← | ← | ← | +| `loglikelihood` | βœ” | ← | ← | ← | ← | +| `likelihood` | βœ” | ← | ← | ← | ← | +| `infer_independent` | ❌ | βœ” | ← | ← | ← | +| `predict_dependent` | ❌ | ❌ | βœ” | βœ” | ← | +| `predict_independent` | ❌ | ❌ | (βœ”) | (βœ”) | ← | + + +Symbols +* ❌ `NotImplementedError` at this level +* βœ” implemented at this class/mixin level +* (βœ”) necessity/feasability depends on the type of model +* ← inherits the implementation diff --git a/docs/source/index.rst b/docs/source/index.rst index 80fa90a..b9a143a 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -39,6 +39,7 @@ You can find autogenerated API documentation below: :maxdepth: 2 :caption: API Reference + calibr8_inheritance.md calibr8_core calibr8_contrib_base calibr8_optimization