diff --git a/doc/example/distributions.ipynb b/doc/example/distributions.ipynb index 6442a548..86235fe1 100644 --- a/doc/example/distributions.ipynb +++ b/doc/example/distributions.ipynb @@ -177,20 +177,11 @@ "source": [ "plot(Prior(NORMAL, (10, 1), bounds=(6, 14), transformation=\"log10\"))\n", "plot(Prior(PARAMETER_SCALE_NORMAL, (10, 1), bounds=(10**6, 10**14), transformation=\"log10\"))\n", - "plot(Prior(LAPLACE, (10, 2), bounds=(6, 14)))\n", - "\n" + "plot(Prior(LAPLACE, (10, 2), bounds=(6, 14)))" ], "id": "581e1ac431860419", "outputs": [], "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": "", - "id": "633733651bbc3ef0", - "outputs": [], - "execution_count": null } ], "metadata": { diff --git a/petab/v1/distributions.py b/petab/v1/distributions.py index be7bad89..401dd0a9 100644 --- a/petab/v1/distributions.py +++ b/petab/v1/distributions.py @@ -4,35 +4,91 @@ import abc import numpy as np -from scipy.stats import laplace, lognorm, loguniform, norm, uniform +from scipy.stats import laplace, norm, uniform __all__ = [ "Distribution", "Normal", - "LogNormal", "Uniform", - "LogUniform", "Laplace", - "LogLaplace", ] class Distribution(abc.ABC): - """A univariate probability distribution.""" + """A univariate probability distribution. + + This class provides a common interface for sampling from and evaluating + the probability density function of a univariate probability distribution. + + The distribution can be transformed by applying a logarithm to the samples + and the PDF. This is useful, e.g., for log-normal distributions. + + :param log: If ``True``, the distribution is transformed to its + corresponding log distribution (e.g., Normal -> LogNormal). + If a float, the distribution is transformed to its corresponding + log distribution with the given base (e.g., Normal -> Log10Normal). + If ``False``, no transformation is applied. + """ + + def __init__(self, log: bool | float = False): + if log is True: + log = np.exp(1) + self._log = log + + def _undo_log(self, x: np.ndarray | float) -> np.ndarray | float: + """Undo the log transformation. + + :param x: The sample to transform. + :return: The transformed sample + """ + if self._log is False: + return x + return self._log**x + + def _apply_log(self, x: np.ndarray | float) -> np.ndarray | float: + """Apply the log transformation. + + :param x: The value to transform. + :return: The transformed value. + """ + if self._log is False: + return x + return np.log(x) / np.log(self._log) - @abc.abstractmethod def sample(self, shape=None) -> np.ndarray: """Sample from the distribution. :param shape: The shape of the sample. :return: A sample from the distribution. """ - ... + sample = self._sample(shape) + return self._undo_log(sample) @abc.abstractmethod + def _sample(self, shape=None) -> np.ndarray: + """Sample from the underlying distribution. + + :param shape: The shape of the sample. + :return: A sample from the underlying distribution, + before applying, e.g., the log transformation. + """ + ... + def pdf(self, x): """Probability density function at x. + :param x: The value at which to evaluate the PDF. + :return: The value of the PDF at ``x``. + """ + # handle the log transformation; see also: + # https://en.wikipedia.org/wiki/Probability_density_function#Scalar_to_scalar + chain_rule_factor = (1 / (x * np.log(self._log))) if self._log else 1 + return self._pdf(self._apply_log(x)) * chain_rule_factor + + @abc.abstractmethod + def _pdf(self, x): + """Probability density function of the underlying distribution at x. + :param x: The value at which to evaluate the PDF. :return: The value of the PDF at ``x``. """ @@ -40,203 +96,138 @@ def pdf(self, x): class Normal(Distribution): - """A normal distribution.""" + """A (log-)normal distribution. + + :param loc: The location parameter of the distribution. + :param scale: The scale parameter of the distribution. + :param truncation: The truncation limits of the distribution. + :param log: If ``True``, the distribution is transformed to a log-normal + distribution. If a float, the distribution is transformed to a + log-normal distribution with the given base. + If ``False``, no transformation is applied. + If a transformation is applied, the location and scale parameters + and the truncation limits are the location, scale and truncation limits + of the underlying normal distribution. + """ def __init__( self, - mean: float, - std: float, + loc: float, + scale: float, truncation: tuple[float, float] | None = None, + log: bool | float = False, ): - super().__init__() - self._mean = mean - self._std = std + super().__init__(log=log) + self._loc = loc + self._scale = scale self._truncation = truncation if truncation is not None: raise NotImplementedError("Truncation is not yet implemented.") def __repr__(self): - return ( - f"Normal(mean={self._mean}, std={self._std}, " - f"truncation={self._truncation})" - ) - - def sample(self, shape=None): - return np.random.normal(loc=self._mean, scale=self._std, size=shape) - - def pdf(self, x): - return norm.pdf(x, loc=self._mean, scale=self._std) - - -class LogNormal(Distribution): - """A log-normal distribution. - - :param mean: The mean of the underlying normal distribution. - :param std: The standard deviation of the underlying normal distribution. - - """ - - def __init__( - self, - mean: float, - std: float, - truncation: tuple[float, float] | None = None, - base: float = np.exp(1), - ): - super().__init__() - self._mean = mean - self._std = std - self._truncation = truncation - self._base = base - - if truncation is not None: - raise NotImplementedError("Truncation is not yet implemented.") + trunc = f", truncation={self._truncation}" if self._truncation else "" + log = f", log={self._log}" if self._log else "" + return f"Normal(loc={self._loc}, scale={self._scale}{trunc}{log})" - if base != np.exp(1): - raise NotImplementedError("Only base e is supported.") + def _sample(self, shape=None): + return np.random.normal(loc=self._loc, scale=self._scale, size=shape) - def __repr__(self): - return ( - f"LogNormal(mean={self._mean}, std={self._std}, " - f"base={self._base}, truncation={self._truncation})" - ) + def _pdf(self, x): + return norm.pdf(x, loc=self._loc, scale=self._scale) - def sample(self, shape=None): - return np.random.lognormal( - mean=self._mean, sigma=self._std, size=shape - ) + @property + def loc(self): + """The location parameter of the underlying distribution.""" + return self._loc - def pdf(self, x): - return lognorm.pdf(x, scale=np.exp(self._mean), s=self._std) + @property + def scale(self): + """The scale parameter of the underlying distribution.""" + return self._scale class Uniform(Distribution): - """A uniform distribution.""" + """A (log-)uniform distribution. + + :param low: The lower bound of the distribution. + :param high: The upper bound of the distribution. + :param log: If ``True``, the distribution is transformed to a log-uniform + distribution. If a float, the distribution is transformed to a + log-uniform distribution with the given base. + If ``False``, no transformation is applied. + If a transformation is applied, the lower and upper bounds are the + lower and upper bounds of the underlying uniform distribution. + """ def __init__( self, low: float, high: float, + *, + log: bool | float = False, ): - super().__init__() + super().__init__(log=log) self._low = low self._high = high def __repr__(self): - return f"Uniform(low={self._low}, high={self._high})" + log = f", log={self._log}" if self._log else "" + return f"Uniform(low={self._low}, high={self._high}{log})" - def sample(self, shape=None): + def _sample(self, shape=None): return np.random.uniform(low=self._low, high=self._high, size=shape) - def pdf(self, x): + def _pdf(self, x): return uniform.pdf(x, loc=self._low, scale=self._high - self._low) -class LogUniform(Distribution): - """A log-uniform distribution. - - :param low: The lower bound of the underlying normal distribution. - :param high: The upper bound of the underlying normal distribution. - """ - - def __init__( - self, - low: float, - high: float, - base: float = np.exp(1), - ): - super().__init__() - self._low = low - self._high = high - self._base = base - # re-scaled distribution parameters as required by - # scipy.stats.loguniform - self._low_internal = np.exp(np.log(base) * low) - self._high_internal = np.exp(np.log(base) * high) - - def __repr__(self): - return ( - f"LogUniform(low={self._low}, high={self._high}, " - f"base={self._base})" - ) - - def sample(self, shape=None): - return loguniform.rvs( - self._low_internal, self._high_internal, size=shape - ) - - def pdf(self, x): - return loguniform.pdf(x, self._low_internal, self._high_internal) - - class Laplace(Distribution): - """A Laplace distribution.""" + """A (log-)Laplace distribution. + + :param loc: The location parameter of the distribution. + :param scale: The scale parameter of the distribution. + :param truncation: The truncation limits of the distribution. + :param log: If ``True``, the distribution is transformed to a log-Laplace + distribution. If a float, the distribution is transformed to a + log-Laplace distribution with the given base. + If ``False``, no transformation is applied. + If a transformation is applied, the location and scale parameters + and the truncation limits are the location, scale and truncation limits + of the underlying Laplace distribution. + """ def __init__( self, loc: float, scale: float, truncation: tuple[float, float] | None = None, + log: bool | float = False, ): - super().__init__() + super().__init__(log=log) self._loc = loc self._scale = scale self._truncation = truncation if truncation is not None: raise NotImplementedError("Truncation is not yet implemented.") - def sample(self, shape=None): + def __repr__(self): + trunc = f", truncation={self._truncation}" if self._truncation else "" + log = f", log={self._log}" if self._log else "" + return f"Laplace(loc={self._loc}, scale={self._scale}{trunc}{log})" + + def _sample(self, shape=None): return np.random.laplace(loc=self._loc, scale=self._scale, size=shape) - def pdf(self, x): + def _pdf(self, x): return laplace.pdf(x, loc=self._loc, scale=self._scale) - -class LogLaplace(Distribution): - """A log-Laplace distribution.""" - - def __init__( - self, - loc: float, - scale: float, - truncation: tuple[float, float] | None = None, - base: float = np.exp(1), - ): - super().__init__() - self._loc = loc - self._scale = scale - self._truncation = truncation - self._base = base - if truncation is not None: - raise NotImplementedError("Truncation is not yet implemented.") - if base != np.exp(1): - raise NotImplementedError("Only base e is supported.") - - def __repr__(self): - return ( - f"LogLaplace(loc={self._loc}, scale={self._scale}, " - f"base={self._base}, truncation={self._truncation})" - ) - @property def loc(self): - """The mean of the underlying Laplace distribution.""" + """The location parameter of the underlying distribution.""" return self._loc @property def scale(self): - """The scale of the underlying Laplace distribution.""" + """The scale parameter of the underlying distribution.""" return self._scale - - def sample(self, shape=None): - return np.exp( - np.random.laplace(loc=self._loc, scale=self._scale, size=shape) - ) - - def pdf(self, x): - return ( - 1 - / (2 * self.scale * x) - * np.exp(-np.abs(np.log(x) - self._loc) / self._scale) - ) diff --git a/petab/v1/priors.py b/petab/v1/priors.py index 554a22a7..f0f37f75 100644 --- a/petab/v1/priors.py +++ b/petab/v1/priors.py @@ -97,21 +97,17 @@ def __init__( case (C.LAPLACE, _) | (C.PARAMETER_SCALE_LAPLACE, C.LIN): self.distribution = Laplace(*parameters) case (C.PARAMETER_SCALE_UNIFORM, C.LOG): - self.distribution = LogUniform(*parameters) + self.distribution = Uniform(*parameters, log=True) case (C.LOG_NORMAL, _) | (C.PARAMETER_SCALE_NORMAL, C.LOG): - self.distribution = LogNormal(*parameters) + self.distribution = Normal(*parameters, log=True) case (C.LOG_LAPLACE, _) | (C.PARAMETER_SCALE_LAPLACE, C.LOG): - self.distribution = LogLaplace(*parameters) + self.distribution = Laplace(*parameters, log=True) case (C.PARAMETER_SCALE_UNIFORM, C.LOG10): - self.distribution = LogUniform(*parameters, base=10) + self.distribution = Uniform(*parameters, log=10) case (C.PARAMETER_SCALE_NORMAL, C.LOG10): - self.distribution = LogNormal( - np.log(10) * parameters[0], np.log(10) * parameters[1] - ) + self.distribution = Normal(*parameters, log=10) case (C.PARAMETER_SCALE_LAPLACE, C.LOG10): - self.distribution = LogLaplace( - np.log(10) * parameters[0], np.log(10) * parameters[1] - ) + self.distribution = Laplace(*parameters, log=10) case _: raise ValueError( "Unsupported distribution type / transformation: " diff --git a/tests/v1/test_distributions.py b/tests/v1/test_distributions.py index c37ef2b0..6e91698b 100644 --- a/tests/v1/test_distributions.py +++ b/tests/v1/test_distributions.py @@ -14,12 +14,14 @@ list( product( [ - Normal(1, 1), - LogNormal(2, 1), + Normal(2, 1), + Normal(2, 1, log=True), + Normal(2, 1, log=10), Uniform(2, 4), - LogUniform(1, 2), + Uniform(-2, 4, log=True), + Uniform(2, 4, log=10), Laplace(1, 2), - LogLaplace(1, 0.5), + Laplace(1, 0.5, log=True), ], [LIN, LOG, LOG10], )