Implement proper truncation for prior distributions

doc/conf.py

@@ -87,6 +87,7 @@
 nb_execution_mode = "force"
 nb_execution_raise_on_error = True
 nb_execution_show_tb = True
+nb_execution_timeout = 90  # max. seconds/cell
 
 source_suffix = {
     ".rst": "restructuredtext",

doc/example/distributions.ipynb

@@ -23,59 +23,99 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "initial_id",
    "metadata": {
     "collapsed": true
    },
-   "outputs": [],
    "source": [
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "import seaborn as sns\n",
     "\n",
     "from petab.v1.C import *\n",
+    "from petab.v1.parameters import unscale\n",
     "from petab.v1.priors import Prior\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",
+    "\n",
+    "def plot_single(\n",
+    "    prior: Prior, scaled: bool = False, ax=None, sample: np.array = None\n",
+    "):\n",
+    "    fig = None\n",
     "    if ax is None:\n",
     "        fig, ax = plt.subplots()\n",
     "\n",
-    "    sample = prior.sample(10000)\n",
+    "    if sample is None:\n",
+    "        sample = prior.sample(20_000)\n",
     "\n",
-    "    # pdf\n",
+    "    # assuming scaled sample\n",
+    "    if not scaled:\n",
+    "        sample = unscale(sample, prior.transformation)\n",
+    "        bounds = prior.bounds\n",
+    "    else:\n",
+    "        bounds = (\n",
+    "            (prior.lb_scaled, prior.ub_scaled)\n",
+    "            if prior.bounds is not None\n",
+    "            else None\n",
+    "        )\n",
+    "\n",
+    "    # plot pdf\n",
     "    xmin = min(\n",
-    "        sample.min(),\n",
-    "        prior.lb_scaled if prior.bounds is not None else sample.min(),\n",
+    "        sample.min(), bounds[0] if prior.bounds is not None else sample.min()\n",
     "    )\n",
     "    xmax = max(\n",
-    "        sample.max(),\n",
-    "        prior.ub_scaled if prior.bounds is not None else sample.max(),\n",
+    "        sample.max(), bounds[1] if prior.bounds is not None else sample.max()\n",
     "    )\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(\n",
+    "            f\"Parameter value on parameter scale ({prior.transformation})\"\n",
+    "        )\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, strict=False))\n",
     "    ax.legend(unique_labels.values(), unique_labels.keys())\n",
-    "    plt.show()"
-   ]
+    "\n",
+    "    if ax is None:\n",
+    "        plt.show()"
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -85,38 +125,38 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "4f09e50a3db06d9f",
    "metadata": {},
-   "outputs": [],
    "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)))"
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
    "id": "dab4b2d1e0f312d8",
    "metadata": {},
-   "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`), the distribution parameters are used as is without applying the `parameterScale` to them. The exception are the `parameterScale*`-type distributions, as explained below. 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"
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "f6192c226f179ef9",
    "metadata": {},
-   "outputs": [],
    "source": [
     "plot(Prior(NORMAL, (10, 2), transformation=LIN))\n",
     "plot(Prior(NORMAL, (10, 2), transformation=LOG))\n",
     "\n",
     "# Note that the log-normal distribution is different\n",
     "#  from a log-transformed normal distribution:\n",
     "plot(Prior(LOG_NORMAL, (10, 2), transformation=LIN))"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
@@ -126,53 +166,69 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "34c95268e8921070",
    "metadata": {},
-   "outputs": [],
    "source": [
     "plot(Prior(LOG_NORMAL, (10, 2), transformation=LOG))\n",
     "plot(Prior(PARAMETER_SCALE_NORMAL, (10, 2)))"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
    "id": "263c9fd31156a4d5",
    "metadata": {},
-   "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 parameters 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`."
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "5ca940bc24312fc6",
    "metadata": {},
-   "outputs": [],
    "source": [
+    "# different, because transformation!=LIN\n",
     "plot(Prior(UNIFORM, (0.01, 2), transformation=LOG10))\n",
     "plot(Prior(PARAMETER_SCALE_UNIFORM, (0.01, 2), transformation=LOG10))\n",
     "\n",
+    "# same, because transformation=LIN\n",
     "plot(Prior(UNIFORM, (0.01, 2), transformation=LIN))\n",
     "plot(Prior(PARAMETER_SCALE_UNIFORM, (0.01, 2), transformation=LIN))"
-   ]
+   ],
+   "outputs": [],
+   "execution_count": null
   },
   {
    "cell_type": "markdown",
    "id": "b1a8b17d765db826",
    "metadata": {},
-   "source": "To prevent the sampled parameters from exceeding the bounds, the sampled parameters are clipped to the bounds. The bounds are defined in the parameter table. Note that the current implementation does not support sampling from a truncated distribution. Instead, the samples are clipped to the bounds. This may introduce unwanted bias, and thus, should only be used with caution (i.e., the bounds should be chosen wide enough):"
+   "source": "The given distributions are truncated at the bounds defined in the parameter table:"
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "4ac42b1eed759bdd",
    "metadata": {},
-   "outputs": [],
    "source": [
+    "plot(Prior(NORMAL, (0, 1), bounds=(-2, 2)))\n",
+    "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(\n",
-    "    Prior(NORMAL, (0, 1), bounds=(-4, 4))\n",
-    ")  # negligible clipping-bias at 4 sigma\n",
-    "plot(Prior(UNIFORM, (0, 1), bounds=(0.1, 0.9)))  # significant clipping-bias"
-   ]
+    "    Prior(\n",
+    "        PARAMETER_SCALE_UNIFORM,\n",
+    "        (-3, 1),\n",
+    "        bounds=(1e-2, 1),\n",
+    "        transformation=LOG10,\n",
+    "    )\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "This results in a constant shift in the probability density, compared to the non-truncated version (https://en.wikipedia.org/wiki/Truncated_distribution), such that the probability density still sums to 1.",
+   "id": "67de0cace55617a2"
   },
   {
    "cell_type": "markdown",
@@ -182,22 +238,24 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
    "id": "581e1ac431860419",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "plot(Prior(NORMAL, (10, 1), bounds=(6, 14), transformation=\"log10\"))\n",
+    "plot(Prior(NORMAL, (10, 1), bounds=(6, 11), transformation=\"log10\"))\n",
     "plot(\n",
     "    Prior(\n",
     "        PARAMETER_SCALE_NORMAL,\n",
-    "        (10, 1),\n",
-    "        bounds=(10**6, 10**14),\n",
+    "        (2, 1),\n",
+    "        bounds=(10**0, 10**3),\n",
     "        transformation=\"log10\",\n",
     "    )\n",
     ")\n",
-    "plot(Prior(LAPLACE, (10, 2), bounds=(6, 14)))"
-   ]
+    "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)))"
+   ],
+   "outputs": [],
+   "execution_count": null
   }
  ],
  "metadata": {

petab/v1/C.py

@@ -208,6 +208,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]