From c8d9bfe9cb29afc09cce6b0a46ca8ee22c7256a3 Mon Sep 17 00:00:00 2001
From: Romain Ragonnet <romain.ragonnet@monash.edu>
Date: Wed, 2 Aug 2023 14:27:52 +1000
Subject: [PATCH] Create notebook to plot sojourn time distributions

---
 .../utils/sojourn_distribution.ipynb          | 81 +++++++++++++++++++
 1 file changed, 81 insertions(+)
 create mode 100644 notebooks/user/rragonnet/utils/sojourn_distribution.ipynb

diff --git a/notebooks/user/rragonnet/utils/sojourn_distribution.ipynb b/notebooks/user/rragonnet/utils/sojourn_distribution.ipynb
new file mode 100644
index 000000000..2d5d92be6
--- /dev/null
+++ b/notebooks/user/rragonnet/utils/sojourn_distribution.ipynb
@@ -0,0 +1,81 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy import stats\n",
+    "from numpy import linspace\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Average sojourn time\n",
+    "mean_period = 5.\n",
+    "\n",
+    "# number of serial compartments to compare (list can be extended)\n",
+    "n_comps_to_compare = [1, 4]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, len(n_comps_to_compare), figsize=(5.*len(n_comps_to_compare), 4))\n",
+    "x_min, x_max = 0., 14\n",
+    "x = linspace(x_min, x_max, 1000)\n",
+    "\n",
+    "for i, n_comps in enumerate(n_comps_to_compare):\n",
+    "    ax = axes[i]\n",
+    "    distri = stats.gamma(a=n_comps, scale=mean_period / n_comps)\n",
+    "    ax.fill_between(x, distri.pdf(x), alpha=.5, color='blue')\n",
+    "    ax.set_xlabel(\"days\")\n",
+    "    ax.set_title(f\"{n_comps} compartment(s)\")\n",
+    "\n",
+    "    # Print some percentages\n",
+    "    print(f\"Using {n_comps} compartment(s):\")\n",
+    "    for quantile in [7., 10., 14.]:\n",
+    "        perc_greater = round(100* (1 - distri.cdf(quantile)))\n",
+    "        print(f\"\\t{perc_greater}% have sojourn time > {quantile} days\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "summer2",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.11"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}