prior always on linear

Author: dweindl

doc/example/distributions.ipynb

@@ -33,42 +33,73 @@
     "\n",
     "from petab.v1.C import *\n",
     "from petab.v1.priors import Prior\n",
+    "from petab.v1.parameters import scale, unscale\n",
+    "\n",
     "\n",
     "sns.set_style(None)\n",
     "\n",
     "\n",
-    "def plot(prior: Prior, ax=None):\n",
+    "def plot(prior: Prior):\n",
     "    \"\"\"Visualize a distribution.\"\"\"\n",
+    "    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "    sample = prior.sample(20_000, x_scaled=True)\n",
+    "\n",
+    "    fig.suptitle(str(prior))\n",
+    "\n",
+    "    plot_single(prior, ax=ax1, sample=sample, scaled=False)\n",
+    "    plot_single(prior, ax=ax2, sample=sample, scaled=True)\n",
+    "    plt.tight_layout()\n",
+    "    plt.show()\n",
+    "\n",
+    "def plot_single(prior: Prior, scaled: bool = False, ax=None, sample: np.array = None):\n",
+    "    fig = None\n",
     "    if ax is None:\n",
     "        fig, ax = plt.subplots()\n",
     "\n",
-    "    sample = prior.sample(20_000)\n",
+    "    if sample is None:\n",
+    "        sample = prior.sample(20_000)\n",
+    "\n",
+    "    # assuming scaled sample\n",
+    "    if not scaled:\n",
+    "        sample = unscale(sample, prior.transformation)\n",
+    "        bounds = prior.bounds\n",
+    "    else:\n",
+    "        bounds = (prior.lb_scaled, prior.ub_scaled) if prior.bounds is not None else None\n",
     "\n",
-    "    # pdf\n",
-    "    xmin = min(sample.min(), prior.lb_scaled if prior.bounds is not None else sample.min())\n",
-    "    xmax = max(sample.max(), prior.ub_scaled if prior.bounds is not None else sample.max())\n",
+    "    # plot pdf\n",
+    "    xmin = min(sample.min(), bounds[0] if prior.bounds is not None else sample.min())\n",
+    "    xmax = max(sample.max(), bounds[1] if prior.bounds is not None else sample.max())\n",
     "    padding = 0.1 * (xmax - xmin)\n",
     "    xmin -= padding\n",
     "    xmax += padding\n",
     "    x = np.linspace(xmin, xmax, 500)\n",
-    "    y = prior.pdf(x)\n",
+    "    y = prior.pdf(x, x_scaled=scaled, rescale=scaled)\n",
     "    ax.plot(x, y, color='red', label='pdf')\n",
     "\n",
     "    sns.histplot(sample, stat='density', ax=ax, label=\"sample\")\n",
     "\n",
-    "    # bounds\n",
+    "    # plot bounds\n",
     "    if prior.bounds is not None:\n",
-    "        for bound in (prior.lb_scaled, prior.ub_scaled):\n",
+    "        for bound in bounds:\n",
     "            if bound is not None and np.isfinite(bound):\n",
     "                ax.axvline(bound, color='black', linestyle='--', label='bound')\n",
     "\n",
-    "    ax.set_title(str(prior))\n",
-    "    ax.set_xlabel('Parameter value on the parameter scale')\n",
+    "    if fig is not None:\n",
+    "        ax.set_title(str(prior))\n",
+    "\n",
+    "    if scaled:\n",
+    "        ax.set_xlabel(f'Parameter value on parameter scale ({prior.transformation})')\n",
+    "        ax.set_ylabel(\"Rescaled density\")\n",
+    "    else:\n",
+    "        ax.set_xlabel('Parameter value')\n",
+    "\n",
     "    ax.grid(False)\n",
     "    handles, labels = ax.get_legend_handles_labels()\n",
     "    unique_labels = dict(zip(labels, handles))\n",
     "    ax.legend(unique_labels.values(), unique_labels.keys())\n",
-    "    plt.show()"
+    "\n",
+    "    if ax is None:\n",
+    "        plt.show()\n"
    ],
    "id": "initial_id",
    "outputs": [],
@@ -84,11 +115,11 @@
    "metadata": {},
    "cell_type": "code",
    "source": [
-    "plot(Prior(UNIFORM, (0, 1)))\n",
-    "plot(Prior(NORMAL, (0, 1)))\n",
-    "plot(Prior(LAPLACE, (0, 1)))\n",
-    "plot(Prior(LOG_NORMAL, (0, 1)))\n",
-    "plot(Prior(LOG_LAPLACE, (1, 0.5)))"
+    "plot_single(Prior(UNIFORM, (0, 1)))\n",
+    "plot_single(Prior(NORMAL, (0, 1)))\n",
+    "plot_single(Prior(LAPLACE, (0, 1)))\n",
+    "plot_single(Prior(LOG_NORMAL, (0, 1)))\n",
+    "plot_single(Prior(LOG_LAPLACE, (1, 0.5)))"
    ],
    "id": "4f09e50a3db06d9f",
    "outputs": [],
@@ -97,7 +128,7 @@
   {
    "metadata": {},
    "cell_type": "markdown",
-   "source": "If a parameter scale is specified (`parameterScale=lin|log|log10` not a `parameterScale*`-type distribution), the sample is transformed accordingly (but not the distribution parameters):\n",
+   "source": "If a parameter scale is specified (`parameterScale=lin|log|log10`) and the chosen distribution is not a `parameterScale*`-type distribution, then the distribution parameters are taken as is, i.e., the `parameterScale` is not applied to the distribution parameters. In the context of PEtab prior distributions, `parameterScale` will only be used for the start point sampling for optimization, where the sample will be transformed accordingly. This is demonstrated below. The left plot always shows the prior distribution for unscaled parameter values, and the right plot shows the prior distribution for scaled parameter values. Note that in the objective function, the prior is always on the unscaled parameters.\n",
    "id": "dab4b2d1e0f312d8"
   },
   {
@@ -134,7 +165,7 @@
   {
    "metadata": {},
    "cell_type": "markdown",
-   "source": "Prior distributions can also be defined on the parameter scale by using the types `parameterScaleUniform`, `parameterScaleNormal` or `parameterScaleLaplace`. In these cases, 1) the distribution parameter are interpreted on the transformed parameter scale, and 2) a sample from the given distribution is used directly, without applying any transformation according to `parameterScale` (this implies, that for `parameterScale=lin`, there is no difference between `parameterScaleUniform` and `uniform`):",
+   "source": "Prior distributions can also be defined on the scaled parameters (i.e., transformed according to `parameterScale`) by using the types `parameterScaleUniform`, `parameterScaleNormal` or `parameterScaleLaplace`. In these cases, the distribution parameter are interpreted on the transformed parameter scale (but not the parameter bounds, see below). This implies, that for `parameterScale=lin`, there is no difference between `parameterScaleUniform` and `uniform`.",
    "id": "263c9fd31156a4d5"
   },
   {
@@ -167,7 +198,7 @@
     "plot(Prior(UNIFORM, (0, 1), bounds=(0.1, 0.9)))\n",
     "plot(Prior(UNIFORM, (1e-8, 1), bounds=(0.1, 0.9), transformation=LOG10))\n",
     "plot(Prior(LAPLACE, (0, 1), bounds=(-0.5, 0.5)))\n",
-    "plot(Prior(PARAMETER_SCALE_UNIFORM, (-3, 1), bounds=(1e-2, 1), transformation=LOG10))\n"
+    "plot(Prior(PARAMETER_SCALE_UNIFORM, (-3, 1), bounds=(1e-2, 1), transformation=LOG10))"
    ],
    "id": "4ac42b1eed759bdd",
    "outputs": [],
@@ -184,7 +215,7 @@
    "cell_type": "code",
    "source": [
     "plot(Prior(NORMAL, (10, 1), bounds=(6, 11), transformation=\"log10\"))\n",
-    "plot(Prior(PARAMETER_SCALE_NORMAL, (10, 1), bounds=(10**9, 10**14), transformation=\"log10\"))\n",
+    "plot(Prior(PARAMETER_SCALE_NORMAL, (2, 1), bounds=(10**0, 10**3), transformation=\"log10\"))\n",
     "plot(Prior(LAPLACE, (10, 2), bounds=(6, 14)))\n",
     "plot(Prior(LOG_LAPLACE, (1, 0.5), bounds=(0.5, 8)))\n",
     "plot(Prior(LOG_NORMAL, (2, 1), bounds=(0.5, 8)))"

petab/v1/C.py

@@ -207,6 +207,13 @@
     PARAMETER_SCALE_LAPLACE,
 ]
 
+#: parameterScale*-type prior distributions
+PARAMETER_SCALE_PRIOR_TYPES = [
+    PARAMETER_SCALE_UNIFORM,
+    PARAMETER_SCALE_NORMAL,
+    PARAMETER_SCALE_LAPLACE,
+]
+
 #: Supported noise distributions
 NOISE_MODELS = [NORMAL, LAPLACE]
 

petab/v1/distributions.py

@@ -168,11 +168,11 @@ def _pdf_transformed_untruncated(self, x) -> np.ndarray | float:
 
         # 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._logbase))) if self._logbase else 1
-        )
+        with np.errstate(invalid="ignore", divide="ignore"):
+            chain_rule_factor = (
+                (1 / (x * np.log(self._logbase))) if self._logbase else 1
+            )
 
-        with np.errstate(invalid="ignore"):
             return np.where(
                 x > 0,
                 self._pdf_untransformed_untruncated(self._log(x))
@@ -242,6 +242,19 @@ def _ppf_transformed_untruncated(self, q) -> np.ndarray | float:
         """
         return self._exp(self._ppf_untransformed_untruncated(q))
 
+    def ppf(self, q) -> np.ndarray | float:
+        """Percent point function at q.
+
+        :param q: The quantile at which to evaluate the PPF.
+        :return: The value of the PPF at ``q``.
+        """
+        if self._trunc is None:
+            return self._ppf_transformed_untruncated(q)
+
+        # Adjust quantiles to account for truncation
+        adjusted_q = self._cd_low + q * (self._cd_high - self._cd_low)
+        return self._ppf_transformed_untruncated(adjusted_q)
+
     def _inverse_transform_sample(self, shape) -> np.ndarray | float:
         """Generate an inverse transform sample from the transformed and
         truncated distribution.

petab/v1/priors.py

@@ -161,24 +161,31 @@ def type(self) -> str:
 
     @property
     def parameters(self) -> tuple:
+        """The parameters of the distribution."""
         return self._parameters
 
     @property
     def bounds(self) -> tuple[float, float] | None:
+        """The non-scaled bounds of the distribution."""
         return self._bounds
 
     @property
     def transformation(self) -> str:
+        """The `parameterScale`."""
         return self._transformation
 
-    def sample(self, shape=None) -> np.ndarray:
+    def sample(self, shape=None, x_scaled=False) -> np.ndarray | float:
         """Sample from the distribution.
 
         :param shape: The shape of the sample.
+        :param x_scaled: Whether the sample should be on the parameter scale.
         :return: A sample from the distribution.
         """
         raw_sample = self.distribution.sample(shape)
-        return self._scale_sample(raw_sample)
+        if x_scaled:
+            return self._scale_sample(raw_sample)
+        else:
+            return raw_sample
 
     def _scale_sample(self, sample):
         """Scale the sample to the parameter space"""
@@ -196,14 +203,8 @@ def ub_scaled(self) -> float:
         """The upper bound on the parameter scale."""
         return scale(self.bounds[1], self.transformation)
 
-    def pdf(self, x) -> np.ndarray | float:
-        """Probability density function at x.
-
-        :param x: The value at which to evaluate the PDF.
-            ``x`` is assumed to be on the parameter scale.
-        :return: The value of the PDF at ``x``. ``x`` is assumed to be on the
-            parameter scale.
-        """
+    def _chain_rule_coeff(self, x) -> np.ndarray | float:
+        """The chain rule coefficient for the transformation at x."""
         x = unscale(x, self.transformation)
 
         # scale the PDF to the parameter scale
@@ -216,36 +217,50 @@ def pdf(self, x) -> np.ndarray | float:
         else:
             raise ValueError(f"Unknown transformation: {self.transformation}")
 
-        return self.distribution.pdf(x) * coeff
+        return coeff
+
+    def pdf(
+        self, x, x_scaled: bool = False, rescale=False
+    ) -> np.ndarray | float:
+        """Probability density function at x.
+
+        This accounts for truncation, independent of the `bounds_truncate`
+        parameter.
 
-    def neglogprior(self, x) -> np.ndarray | float:
+        :param x: The value at which to evaluate the PDF.
+            ``x`` is assumed to be on the parameter scale.
+        :param x_scaled: Whether ``x`` is on the parameter scale.
+        :param rescale: Whether to rescale the PDF to integrate to 1 on the
+            parameter scale. Only used if ``x_scaled`` is ``True``.
+        :return: The value of the PDF at ``x``.
+        """
+        if x_scaled:
+            coeff = self._chain_rule_coeff(x) if rescale else 1
+            x = unscale(x, self.transformation)
+            return self.distribution.pdf(x) * coeff
+
+        return self.distribution.pdf(x)
+
+    def neglogprior(
+        self, x: np.array | float, x_scaled: bool = False
+    ) -> np.ndarray | float:
         """Negative log-prior at x.
 
         :param x: The value at which to evaluate the negative log-prior.
-            ``x`` is assumed to be on the parameter scale.
+        :param x_scaled: Whether ``x`` is on the parameter scale.
+            Note that the prior is always evaluated on the non-scaled
+            parameters.
         :return: The negative log-prior at ``x``.
         """
-        # FIXME: the prior is always defined on linear scale
         if self._bounds_truncate:
             # the truncation is handled by the distribution
-            return -np.log(self.pdf(x))
+            # the prior is always evaluated on the non-scaled parameters
+            return -np.log(self.pdf(x, x_scaled=x_scaled, rescale=False))
 
         # we want to evaluate the prior on the untruncated distribution
-        x = unscale(x, self.transformation)
-
-        # scale the PDF to the parameter scale
-        if self.transformation == C.LIN:
-            coeff = 1
-        elif self.transformation == C.LOG10:
-            coeff = x * np.log(10)
-        elif self.transformation == C.LOG:
-            coeff = x
-        else:
-            raise ValueError(f"Unknown transformation: {self.transformation}")
-
-        return -np.log(
-            self.distribution._pdf_transformed_untruncated(x) * coeff
-        )
+        if x_scaled:
+            x = unscale(x, self.transformation)
+        return -np.log(self.distribution._pdf_transformed_untruncated(x))
 
     @staticmethod
     def from_par_dict(
@@ -339,6 +354,7 @@ def priors_to_measurements(problem: Problem):
         return new_problem
 
     def scaled_observable_formula(parameter_id, parameter_scale):
+        # The location parameter of the prior
         if parameter_scale == LIN:
             return parameter_id
         if parameter_scale == LOG:
@@ -367,6 +383,12 @@ def scaled_observable_formula(parameter_id, parameter_scale):
             #  offset
             raise NotImplementedError("Uniform priors are not supported.")
 
+        if prior_type not in (C.NORMAL, C.LAPLACE):
+            # we can't (easily) handle parameterScale* priors or log*-priors
+            raise NotImplementedError(
+                f"Objective prior type {prior_type} is not implemented."
+            )
+
         parameter_id = row.name
         prior_parameters = tuple(
             map(
@@ -391,7 +413,9 @@ def scaled_observable_formula(parameter_id, parameter_scale):
             OBSERVABLE_ID: new_obs_id,
             OBSERVABLE_FORMULA: scaled_observable_formula(
                 parameter_id,
-                parameter_scale if "parameterScale" in prior_type else LIN,
+                parameter_scale
+                if prior_type in C.PARAMETER_SCALE_PRIOR_TYPES
+                else LIN,
             ),
             NOISE_FORMULA: f"noiseParameter1_{new_obs_id}",
         }
@@ -400,12 +424,13 @@ def scaled_observable_formula(parameter_id, parameter_scale):
         elif OBSERVABLE_TRANSFORMATION in new_problem.observable_df:
             # only set default if the column is already present
             new_observable[OBSERVABLE_TRANSFORMATION] = LIN
-
+        # type of the underlying distribution
         if prior_type in (NORMAL, PARAMETER_SCALE_NORMAL, LOG_NORMAL):
             new_observable[NOISE_DISTRIBUTION] = NORMAL
         elif prior_type in (LAPLACE, PARAMETER_SCALE_LAPLACE, LOG_LAPLACE):
             new_observable[NOISE_DISTRIBUTION] = LAPLACE
         else:
+            # we can't (easily) handle uniform priors in PEtab v1
             raise NotImplementedError(
                 f"Objective prior type {prior_type} is not implemented."
             )

petab/v1/sampling.py

@@ -80,7 +80,7 @@ def sample_parameter_startpoints(
         [
             Prior.from_par_dict(
                 row, type_="initialization", bounds_truncate=True
-            ).sample(n_starts)
+            ).sample(n_starts, x_scaled=True)
             for row in par_to_estimate.to_dict("records")
         ]
     ).T