Skip to content

Commit d9e8193

Browse files
dweindldweindl
committed
Refactor priors
Introduce `Distribution` classes for handling prior distributions and sampling from them. Later on this can be extended to noise models for measurements.
1 parent 9a4efb4 commit d9e8193

File tree

5 files changed

+540
-85
lines changed

5 files changed

+540
-85
lines changed

‎doc/example.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ The following examples should help to get a better idea of how to use the PEtab
1010

1111
example/example_petablint.ipynb
1212
example/example_visualization.ipynb
13+
example/distributions.ipynb
1314

1415
Examples of systems biology parameter estimation problems specified in PEtab
1516
can be found in the `systems biology benchmark model collection <https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab>`_.

‎doc/example/distributions.ipynb‎

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
{
2+
"cells": [
3+
{
4+
"metadata": {},
5+
"cell_type": "markdown",
6+
"source": [
7+
"# Prior distributions in PEtab\n",
8+
"\n",
9+
"This notebook gives a brief overview of the prior distributions in PEtab and how they are represented in the PEtab library.\n",
10+
"\n",
11+
"Prior distributions are used to specify the prior knowledge about the parameters.\n",
12+
"Parameter priors are specified in the parameter table. A prior is defined by its type and its parameters.\n",
13+
"Each prior type has a specific set of parameters. For example, the normal distribution has two parameters: the mean and the standard deviation.\n",
14+
"\n",
15+
"There are two types of priors in PEtab - objective priors and initialization priors:\n",
16+
"\n",
17+
"* *Objective priors* are used to specify the prior knowledge about the parameters that are to be estimated. They will enter the objective function of the optimization problem. They are specified in the `objectivePriorType` and `objectivePriorParameters` columns of the parameter table.\n",
18+
"* *Initialization priors* can be used as a hint for the optimization algorithm. They will not enter the objective function. They are specified in the `initializationPriorType` and `initializationPriorParameters` columns of the parameter table.\n",
19+
"\n",
20+
"\n"
21+
],
22+
"id": "372289411a2aa7b3"
23+
},
24+
{
25+
"cell_type": "code",
26+
"id": "initial_id",
27+
"metadata": {
28+
"collapsed": true
29+
},
30+
"source": [
31+
"from petab.v1.distributions import *\n",
32+
"import seaborn as sns\n",
33+
"import matplotlib.pyplot as plt\n",
34+
"import numpy as np\n",
35+
"from petab.v1.C import *\n",
36+
"\n",
37+
"sns.set_style(None)\n",
38+
"\n",
39+
"\n",
40+
"def plot(distr: Distribution, ax=None):\n",
41+
" \"\"\"Visualize a distribution.\"\"\"\n",
42+
" if ax is None:\n",
43+
" fig, ax = plt.subplots()\n",
44+
"\n",
45+
" sample = distr.sample(10000)\n",
46+
"\n",
47+
" # pdf\n",
48+
" xmin = min(sample.min(), distr.bounds[0] if distr.bounds is not None else sample.min())\n",
49+
" xmax = max(sample.max(), distr.bounds[1] if distr.bounds is not None else sample.max())\n",
50+
" x = np.linspace(xmin, xmax, 500)\n",
51+
" y = distr.pdf(x)\n",
52+
" ax.plot(x, y, color='red', label='pdf')\n",
53+
"\n",
54+
" sns.histplot(sample, stat='density', ax=ax, label=\"sample\")\n",
55+
"\n",
56+
" # bounds\n",
57+
" for bound in distr.bounds or []:\n",
58+
" if bound is not None and np.isfinite(bound):\n",
59+
" ax.axvline(bound, color='black', linestyle='--', label='bound')\n",
60+
"\n",
61+
" ax.set_title(str(distr))\n",
62+
" ax.set_xlabel('Parameter value on the parameter scale')\n",
63+
" ax.grid(False)\n",
64+
" handles, labels = ax.get_legend_handles_labels()\n",
65+
" unique_labels = dict(zip(labels, handles))\n",
66+
" ax.legend(unique_labels.values(), unique_labels.keys())\n",
67+
" plt.show()"
68+
],
69+
"outputs": [],
70+
"execution_count": null
71+
},
72+
{
73+
"metadata": {},
74+
"cell_type": "markdown",
75+
"source": "The basic distributions are the uniform, normal, Laplace, log-normal, and log-laplace distributions:\n",
76+
"id": "db36a4a93622ccb8"
77+
},
78+
{
79+
"metadata": {},
80+
"cell_type": "code",
81+
"source": [
82+
"plot(Uniform(0, 1))\n",
83+
"plot(Normal(0, 1))\n",
84+
"plot(Laplace(0, 1))\n",
85+
"plot(LogNormal(0, 1))\n",
86+
"plot(LogLaplace(1, 0.5))"
87+
],
88+
"id": "4f09e50a3db06d9f",
89+
"outputs": [],
90+
"execution_count": null
91+
},
92+
{
93+
"metadata": {},
94+
"cell_type": "markdown",
95+
"source": "For those distributions, the distribution parameters provided in the parameter table are used as is, independently of the parameter scale. The parameter scale is only used to transform the sampled parameters:",
96+
"id": "dab4b2d1e0f312d8"
97+
},
98+
{
99+
"metadata": {},
100+
"cell_type": "code",
101+
"source": [
102+
"plot(Normal(1, 2, transformation=LIN))\n",
103+
"plot(Normal(1, 2, transformation=LOG))\n",
104+
"# Note that the log-normal is different from the log-transformed normal distribution:\n",
105+
"plot(LogNormal(1, 2, transformation=LIN))"
106+
],
107+
"id": "f6192c226f179ef9",
108+
"outputs": [],
109+
"execution_count": null
110+
},
111+
{
112+
"metadata": {},
113+
"cell_type": "markdown",
114+
"source": "Prior distributions can also be defined on the parameter scale by using the types `parameterScaleUniform`, `parameterScaleNormal` or `parameterScaleLaplace`. In these cases, 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`):",
115+
"id": "263c9fd31156a4d5"
116+
},
117+
{
118+
"metadata": {},
119+
"cell_type": "code",
120+
"source": [
121+
"plot(Uniform(0, 1, transformation=LOG10))\n",
122+
"plot(ParameterScaleUniform(0, 1, transformation=LOG10))"
123+
],
124+
"id": "5ca940bc24312fc6",
125+
"outputs": [],
126+
"execution_count": null
127+
},
128+
{
129+
"metadata": {},
130+
"cell_type": "markdown",
131+
"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, which may introduce unwanted bias:",
132+
"id": "b1a8b17d765db826"
133+
},
134+
{
135+
"metadata": {},
136+
"cell_type": "code",
137+
"source": "plot(Uniform(0, 1, bounds=(0.1, 0.9)))",
138+
"id": "4ac42b1eed759bdd",
139+
"outputs": [],
140+
"execution_count": null
141+
}
142+
],
143+
"metadata": {
144+
"kernelspec": {
145+
"display_name": "Python 3",
146+
"language": "python",
147+
"name": "python3"
148+
},
149+
"language_info": {
150+
"codemirror_mode": {
151+
"name": "ipython",
152+
"version": 2
153+
},
154+
"file_extension": ".py",
155+
"mimetype": "text/x-python",
156+
"name": "python",
157+
"nbconvert_exporter": "python",
158+
"pygments_lexer": "ipython2",
159+
"version": "2.7.6"
160+
}
161+
},
162+
"nbformat": 4,
163+
"nbformat_minor": 5
164+
}

‎doc/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
example
2424

2525

26+
2627
Indices and tables
2728
==================
2829

0 commit comments

Comments
 (0)