diff --git a/10_StochasticContinuous_practical.qmd b/10_StochasticContinuous_practical.qmd
index 521e435..cc72daf 100644
--- a/10_StochasticContinuous_practical.qmd
+++ b/10_StochasticContinuous_practical.qmd
@@ -54,16 +54,16 @@ SIR_gillespie <- function(init_state, parms, tf) {
       time <- time + rexp(n = 1, rate = sum(rates))
       ## check if next event is supposed to happen before end time
       if (time <= tf) {
-        ## generate cumulative sum of rates, to determine the type of the next
-        ## event
-        cumulative_rates <- cumsum(rates)
+        ## determine the type of the next
+        next_event <- sample(x = length(rates), size = 1, prob = rates)
+        ## change to name
+        next_event <- names(rates)[next_event]
         ## determine type of next event
-        type <- runif(n = 1, min = 0, max = sum(rates))
-        if (type < cumulative_rates["infection"]) {
+        if (next_event == "infection") {
           ## infection
           S <- S - 1
           I <- I + 1
-        } else if (type < cumulative_rates["recovery"]) {
+        } else if (next_event == "recovery") {
           ## recovery
           I <- I - 1
           R <- R + 1
diff --git a/10_StochasticContinuous_solutions.qmd b/10_StochasticContinuous_solutions.qmd
index c364cb6..0b6b310 100644
--- a/10_StochasticContinuous_solutions.qmd
+++ b/10_StochasticContinuous_solutions.qmd
@@ -46,16 +46,16 @@ SIR_gillespie <- function(init_state, parms, tf) {
       time <- time + rexp(n = 1, rate = sum(rates))
       ## check if next event is supposed to happen before end time
       if (time <= tf) {
-        ## generate cumulative sum of rates, to determine the type of the next
-        ## event
-        cumulative_rates <- cumsum(rates)
+        ## determine the type of the next
+        next_event <- sample(x = length(rates), size = 1, prob = rates)
+        ## change to name
+        next_event <- names(rates)[next_event]
         ## determine type of next event
-        type <- runif(n = 1, min = 0, max = sum(rates))
-        if (type < cumulative_rates["infection"]) {
+        if (next_event == "infection") {
           ## infection
           S <- S - 1
           I <- I + 1
-        } else if (type < cumulative_rates["recovery"]){
+        } else if (next_event == "recovery") {
           ## recovery
           I <- I - 1
           R <- R + 1
diff --git a/docs/10_StochasticContinuous_practical.html b/docs/10_StochasticContinuous_practical.html
index 3c6147c..89cc7e3 100644
--- a/docs/10_StochasticContinuous_practical.html
+++ b/docs/10_StochasticContinuous_practical.html
@@ -2,12 +2,12 @@
 <html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
 
 <meta charset="utf-8">
-<meta name="generator" content="quarto-1.4.555">
+<meta name="generator" content="quarto-1.5.57">
 
 <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
 
 
-<title>Modern Techniques in Modelling - 10. Stochastic continuous models</title>
+<title>10. Stochastic continuous models – Modern Techniques in Modelling</title>
 <style>
 code{white-space: pre-wrap;}
 span.smallcaps{font-variant: small-caps;}
@@ -33,7 +33,7 @@
 }
 @media print {
 pre > code.sourceCode { white-space: pre-wrap; }
-pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
+pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
 }
 pre.numberSource code
   { counter-reset: source-line 0; }
@@ -101,7 +101,7 @@
   }
 }</script>
 
-  <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
+  <script src="https://cdnjs.cloudflare.com/polyfill/v3/polyfill.min.js?features=es6"></script>
   <script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml-full.js" type="text/javascript"></script>
 
 <script type="text/javascript">
@@ -146,7 +146,7 @@
     </a>
   </div>
             <div id="quarto-search" class="" title="Search"></div>
-          <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
+          <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" role="menu" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
   <span class="navbar-toggler-icon"></span>
 </button>
           <div class="collapse navbar-collapse" id="navbarCollapse">
@@ -164,7 +164,7 @@
 <span class="menu-text">Timetable</span></a>
   </li>  
   <li class="nav-item dropdown ">
-    <a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="button" data-bs-toggle="dropdown" aria-expanded="false">
+    <a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="link" data-bs-toggle="dropdown" aria-expanded="false">
  <span class="menu-text">Practicals</span>
     </a>
     <ul class="dropdown-menu" aria-labelledby="nav-menu-practicals">    
@@ -204,7 +204,7 @@
   </li>  
 </ul>
           </div> <!-- /navcollapse -->
-          <div class="quarto-navbar-tools">
+            <div class="quarto-navbar-tools">
 </div>
       </div> <!-- /container-fluid -->
     </nav>
@@ -299,16 +299,16 @@ <h3 class="anchored" data-anchor-id="setting-up-the-model-algorithm">Setting up
 <span id="cb1-39"><a href="#cb1-39" aria-hidden="true" tabindex="-1"></a>      time <span class="ot">&lt;-</span> time <span class="sc">+</span> <span class="fu">rexp</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">rate =</span> <span class="fu">sum</span>(rates))</span>
 <span id="cb1-40"><a href="#cb1-40" aria-hidden="true" tabindex="-1"></a>      <span class="do">## check if next event is supposed to happen before end time</span></span>
 <span id="cb1-41"><a href="#cb1-41" aria-hidden="true" tabindex="-1"></a>      <span class="cf">if</span> (time <span class="sc">&lt;=</span> tf) {</span>
-<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a>        <span class="do">## generate cumulative sum of rates, to determine the type of the next</span></span>
-<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a>        <span class="do">## event</span></span>
-<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a>        cumulative_rates <span class="ot">&lt;-</span> <span class="fu">cumsum</span>(rates)</span>
-<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine type of next event</span></span>
-<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a>        type <span class="ot">&lt;-</span> <span class="fu">runif</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">min =</span> <span class="dv">0</span>, <span class="at">max =</span> <span class="fu">sum</span>(rates))</span>
-<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a>        <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"infection"</span>]) {</span>
+<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine the type of the next</span></span>
+<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a>        next_event <span class="ot">&lt;-</span> <span class="fu">sample</span>(<span class="at">x =</span> <span class="fu">length</span>(rates), <span class="at">size =</span> <span class="dv">1</span>, <span class="at">prob =</span> rates)</span>
+<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a>        <span class="do">## change to name</span></span>
+<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a>        next_event <span class="ot">&lt;-</span> <span class="fu">names</span>(rates)[next_event]</span>
+<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine type of next event</span></span>
+<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a>        <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"infection"</span>) {</span>
 <span id="cb1-48"><a href="#cb1-48" aria-hidden="true" tabindex="-1"></a>          <span class="do">## infection</span></span>
 <span id="cb1-49"><a href="#cb1-49" aria-hidden="true" tabindex="-1"></a>          S <span class="ot">&lt;-</span> S <span class="sc">-</span> <span class="dv">1</span></span>
 <span id="cb1-50"><a href="#cb1-50" aria-hidden="true" tabindex="-1"></a>          I <span class="ot">&lt;-</span> I <span class="sc">+</span> <span class="dv">1</span></span>
-<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a>        } <span class="cf">else</span> <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"recovery"</span>]) {</span>
+<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a>        } <span class="cf">else</span> <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"recovery"</span>) {</span>
 <span id="cb1-52"><a href="#cb1-52" aria-hidden="true" tabindex="-1"></a>          <span class="do">## recovery</span></span>
 <span id="cb1-53"><a href="#cb1-53" aria-hidden="true" tabindex="-1"></a>          I <span class="ot">&lt;-</span> I <span class="sc">-</span> <span class="dv">1</span></span>
 <span id="cb1-54"><a href="#cb1-54" aria-hidden="true" tabindex="-1"></a>          R <span class="ot">&lt;-</span> R <span class="sc">+</span> <span class="dv">1</span></span>
@@ -682,18 +682,7 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
     }
     return false;
   }
-  const clipboard = new window.ClipboardJS('.code-copy-button', {
-    text: function(trigger) {
-      const codeEl = trigger.previousElementSibling.cloneNode(true);
-      for (const childEl of codeEl.children) {
-        if (isCodeAnnotation(childEl)) {
-          childEl.remove();
-        }
-      }
-      return codeEl.innerText;
-    }
-  });
-  clipboard.on('success', function(e) {
+  const onCopySuccess = function(e) {
     // button target
     const button = e.trigger;
     // don't keep focus
@@ -725,7 +714,29 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
     }, 1000);
     // clear code selection
     e.clearSelection();
+  }
+  const getTextToCopy = function(trigger) {
+      const codeEl = trigger.previousElementSibling.cloneNode(true);
+      for (const childEl of codeEl.children) {
+        if (isCodeAnnotation(childEl)) {
+          childEl.remove();
+        }
+      }
+      return codeEl.innerText;
+  }
+  const clipboard = new window.ClipboardJS('.code-copy-button:not([data-in-quarto-modal])', {
+    text: getTextToCopy
   });
+  clipboard.on('success', onCopySuccess);
+  if (window.document.getElementById('quarto-embedded-source-code-modal')) {
+    // For code content inside modals, clipBoardJS needs to be initialized with a container option
+    // TODO: Check when it could be a function (https://github.com/zenorocha/clipboard.js/issues/860)
+    const clipboardModal = new window.ClipboardJS('.code-copy-button[data-in-quarto-modal]', {
+      text: getTextToCopy,
+      container: window.document.getElementById('quarto-embedded-source-code-modal')
+    });
+    clipboardModal.on('success', onCopySuccess);
+  }
     var localhostRegex = new RegExp(/^(?:http|https):\/\/localhost\:?[0-9]*\//);
     var mailtoRegex = new RegExp(/^mailto:/);
       var filterRegex = new RegExp('/' + window.location.host + '/');
@@ -733,7 +744,7 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
         return filterRegex.test(href) || localhostRegex.test(href) || mailtoRegex.test(href);
     }
     // Inspect non-navigation links and adorn them if external
- 	var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool)');
+ 	var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool):not(.about-link)');
     for (var i=0; i<links.length; i++) {
       const link = links[i];
       if (!isInternal(link.href)) {
diff --git a/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-5-1.png b/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-5-1.png
index 5adf522..bd97809 100644
Binary files a/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-5-1.png and b/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-6-1.png b/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-6-1.png
index 53e5eda..78327ca 100644
Binary files a/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-6-1.png and b/docs/10_StochasticContinuous_practical_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/docs/10_StochasticContinuous_solutions.html b/docs/10_StochasticContinuous_solutions.html
index 53aa52d..7e25070 100644
--- a/docs/10_StochasticContinuous_solutions.html
+++ b/docs/10_StochasticContinuous_solutions.html
@@ -2,12 +2,12 @@
 <html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
 
 <meta charset="utf-8">
-<meta name="generator" content="quarto-1.4.555">
+<meta name="generator" content="quarto-1.5.57">
 
 <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
 
 
-<title>Modern Techniques in Modelling - 10. Stochastic continuous time models (solutions)</title>
+<title>10. Stochastic continuous time models (solutions) – Modern Techniques in Modelling</title>
 <style>
 code{white-space: pre-wrap;}
 span.smallcaps{font-variant: small-caps;}
@@ -33,7 +33,7 @@
 }
 @media print {
 pre > code.sourceCode { white-space: pre-wrap; }
-pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
+pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
 }
 pre.numberSource code
   { counter-reset: source-line 0; }
@@ -117,7 +117,7 @@
     </a>
   </div>
             <div id="quarto-search" class="" title="Search"></div>
-          <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
+          <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" role="menu" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
   <span class="navbar-toggler-icon"></span>
 </button>
           <div class="collapse navbar-collapse" id="navbarCollapse">
@@ -135,7 +135,7 @@
 <span class="menu-text">Timetable</span></a>
   </li>  
   <li class="nav-item dropdown ">
-    <a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="button" data-bs-toggle="dropdown" aria-expanded="false">
+    <a class="nav-link dropdown-toggle" href="#" id="nav-menu-practicals" role="link" data-bs-toggle="dropdown" aria-expanded="false">
  <span class="menu-text">Practicals</span>
     </a>
     <ul class="dropdown-menu" aria-labelledby="nav-menu-practicals">    
@@ -175,7 +175,7 @@
   </li>  
 </ul>
           </div> <!-- /navcollapse -->
-          <div class="quarto-navbar-tools">
+            <div class="quarto-navbar-tools">
 </div>
       </div> <!-- /container-fluid -->
     </nav>
@@ -260,16 +260,16 @@ <h2 class="anchored" data-anchor-id="practical-1.-stochastic-simulation-with-the
 <span id="cb1-39"><a href="#cb1-39" aria-hidden="true" tabindex="-1"></a>      time <span class="ot">&lt;-</span> time <span class="sc">+</span> <span class="fu">rexp</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">rate =</span> <span class="fu">sum</span>(rates))</span>
 <span id="cb1-40"><a href="#cb1-40" aria-hidden="true" tabindex="-1"></a>      <span class="do">## check if next event is supposed to happen before end time</span></span>
 <span id="cb1-41"><a href="#cb1-41" aria-hidden="true" tabindex="-1"></a>      <span class="cf">if</span> (time <span class="sc">&lt;=</span> tf) {</span>
-<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a>        <span class="do">## generate cumulative sum of rates, to determine the type of the next</span></span>
-<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a>        <span class="do">## event</span></span>
-<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a>        cumulative_rates <span class="ot">&lt;-</span> <span class="fu">cumsum</span>(rates)</span>
-<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine type of next event</span></span>
-<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a>        type <span class="ot">&lt;-</span> <span class="fu">runif</span>(<span class="at">n =</span> <span class="dv">1</span>, <span class="at">min =</span> <span class="dv">0</span>, <span class="at">max =</span> <span class="fu">sum</span>(rates))</span>
-<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a>        <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"infection"</span>]) {</span>
+<span id="cb1-42"><a href="#cb1-42" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine the type of the next</span></span>
+<span id="cb1-43"><a href="#cb1-43" aria-hidden="true" tabindex="-1"></a>        next_event <span class="ot">&lt;-</span> <span class="fu">sample</span>(<span class="at">x =</span> <span class="fu">length</span>(rates), <span class="at">size =</span> <span class="dv">1</span>, <span class="at">prob =</span> rates)</span>
+<span id="cb1-44"><a href="#cb1-44" aria-hidden="true" tabindex="-1"></a>        <span class="do">## change to name</span></span>
+<span id="cb1-45"><a href="#cb1-45" aria-hidden="true" tabindex="-1"></a>        next_event <span class="ot">&lt;-</span> <span class="fu">names</span>(rates)[next_event]</span>
+<span id="cb1-46"><a href="#cb1-46" aria-hidden="true" tabindex="-1"></a>        <span class="do">## determine type of next event</span></span>
+<span id="cb1-47"><a href="#cb1-47" aria-hidden="true" tabindex="-1"></a>        <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"infection"</span>) {</span>
 <span id="cb1-48"><a href="#cb1-48" aria-hidden="true" tabindex="-1"></a>          <span class="do">## infection</span></span>
 <span id="cb1-49"><a href="#cb1-49" aria-hidden="true" tabindex="-1"></a>          S <span class="ot">&lt;-</span> S <span class="sc">-</span> <span class="dv">1</span></span>
 <span id="cb1-50"><a href="#cb1-50" aria-hidden="true" tabindex="-1"></a>          I <span class="ot">&lt;-</span> I <span class="sc">+</span> <span class="dv">1</span></span>
-<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a>        } <span class="cf">else</span> <span class="cf">if</span> (type <span class="sc">&lt;</span> cumulative_rates[<span class="st">"recovery"</span>]){</span>
+<span id="cb1-51"><a href="#cb1-51" aria-hidden="true" tabindex="-1"></a>        } <span class="cf">else</span> <span class="cf">if</span> (next_event <span class="sc">==</span> <span class="st">"recovery"</span>) {</span>
 <span id="cb1-52"><a href="#cb1-52" aria-hidden="true" tabindex="-1"></a>          <span class="do">## recovery</span></span>
 <span id="cb1-53"><a href="#cb1-53" aria-hidden="true" tabindex="-1"></a>          I <span class="ot">&lt;-</span> I <span class="sc">-</span> <span class="dv">1</span></span>
 <span id="cb1-54"><a href="#cb1-54" aria-hidden="true" tabindex="-1"></a>          R <span class="ot">&lt;-</span> R <span class="sc">+</span> <span class="dv">1</span></span>
@@ -527,18 +527,7 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
     }
     return false;
   }
-  const clipboard = new window.ClipboardJS('.code-copy-button', {
-    text: function(trigger) {
-      const codeEl = trigger.previousElementSibling.cloneNode(true);
-      for (const childEl of codeEl.children) {
-        if (isCodeAnnotation(childEl)) {
-          childEl.remove();
-        }
-      }
-      return codeEl.innerText;
-    }
-  });
-  clipboard.on('success', function(e) {
+  const onCopySuccess = function(e) {
     // button target
     const button = e.trigger;
     // don't keep focus
@@ -570,7 +559,29 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
     }, 1000);
     // clear code selection
     e.clearSelection();
+  }
+  const getTextToCopy = function(trigger) {
+      const codeEl = trigger.previousElementSibling.cloneNode(true);
+      for (const childEl of codeEl.children) {
+        if (isCodeAnnotation(childEl)) {
+          childEl.remove();
+        }
+      }
+      return codeEl.innerText;
+  }
+  const clipboard = new window.ClipboardJS('.code-copy-button:not([data-in-quarto-modal])', {
+    text: getTextToCopy
   });
+  clipboard.on('success', onCopySuccess);
+  if (window.document.getElementById('quarto-embedded-source-code-modal')) {
+    // For code content inside modals, clipBoardJS needs to be initialized with a container option
+    // TODO: Check when it could be a function (https://github.com/zenorocha/clipboard.js/issues/860)
+    const clipboardModal = new window.ClipboardJS('.code-copy-button[data-in-quarto-modal]', {
+      text: getTextToCopy,
+      container: window.document.getElementById('quarto-embedded-source-code-modal')
+    });
+    clipboardModal.on('success', onCopySuccess);
+  }
     var localhostRegex = new RegExp(/^(?:http|https):\/\/localhost\:?[0-9]*\//);
     var mailtoRegex = new RegExp(/^mailto:/);
       var filterRegex = new RegExp('/' + window.location.host + '/');
@@ -578,7 +589,7 @@ <h2 class="anchored" data-anchor-id="practical-2.-the-adaptivetau-package">Pract
         return filterRegex.test(href) || localhostRegex.test(href) || mailtoRegex.test(href);
     }
     // Inspect non-navigation links and adorn them if external
- 	var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool)');
+ 	var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool):not(.about-link)');
     for (var i=0; i<links.length; i++) {
       const link = links[i];
       if (!isInternal(link.href)) {
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-1.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-1.png
index 1fc4b22..13b69d5 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-1.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-1.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-2.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-2.png
index f1bd0a3..afd7815 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-2.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-2.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-3.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-3.png
index 2b233a3..6559297 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-3.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-3.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-4.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-4.png
index 5917ac9..c703b6b 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-4.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-4.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-5.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-5.png
index 8463f0d..bd7a0e7 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-5.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-5.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-6.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-6.png
index 7f76c67..6ce0135 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-6.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-1-6.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-1.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-1.png
index 9f4b753..9c4bb41 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-1.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-1.png differ
diff --git a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-2.png b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-2.png
index 60fa885..0e0def8 100644
Binary files a/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-2.png and b/docs/10_StochasticContinuous_solutions_files/figure-html/unnamed-chunk-2-2.png differ
diff --git a/docs/search.json b/docs/search.json
index 49f7440..18bbf10 100644
--- a/docs/search.json
+++ b/docs/search.json
@@ -88,14 +88,14 @@
     "href": "10_StochasticContinuous_solutions.html",
     "title": "10. Stochastic continuous time models (solutions)",
     "section": "",
-    "text": "library(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## generate cumulative sum of rates, to determine the type of the next\n        ## event\n        cumulative_rates &lt;- cumsum(rates)\n        ## determine type of next event\n        type &lt;- runif(n = 1, min = 0, max = sum(rates))\n        if (type &lt; cumulative_rates[\"infection\"]) {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (type &lt; cumulative_rates[\"recovery\"]){\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\n## Plot the result\nggplot(r) +\n    geom_line(aes(time, I, colour = \"I\"))\n## Re-run this from the line `r &lt;- SIR_gillespie(...)` a few times to \n## convince yourself the output is different every time.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\n## We store the simulations in a data frame, traj, which\n## contains the results from multiple simulation runs and an additional column\n## that represents the simulation index\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\n## Next, plot the multiple simulation runs\nggplot(mlr,\n       aes(x = time, y = value, group = i, color = compartment)) +\n  geom_line() +\n  facet_wrap(~compartment)\n## Plot the distribution of overall outbreak sizes.\n\n## create data frame of outbreak sizes\noutbreak_sizes &lt;- mlr %&gt;%\n  filter(compartment==\"R\") %&gt;%\n  group_by(i) %&gt;%\n  filter(time==max(time)) %&gt;%\n  select(i, size = value)\n\n## plot it as a histogram\nggplot(outbreak_sizes, aes(x = size)) +\n  geom_histogram()\n## determine number of large outbreaks (defined as larger than 50)\noutbreak_sizes %&gt;%\n  filter(size &gt; 50) %&gt;%\n  nrow\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(\n    mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n## We now compare the two averages to the deterministic trajectory\n\n## Model function (from ODE practical)\nlibrary(deSolve)\nSIR_model &lt;- function(times, state, parms){\n  ## Define variables\n  S &lt;- state[\"S\"]\n  I &lt;- state[\"I\"]\n  R &lt;- state[\"R\"]\n  N &lt;- S + I + R\n                                        # Extract parameters\n  beta &lt;- parms[\"beta\"]\n  gamma &lt;- parms[\"gamma\"]\n                                        # Define differential equations\n  dS &lt;- - (beta * S * I) / N\n  dI &lt;- (beta * S * I) / N - gamma * I\n  dR &lt;- gamma * I\n  res &lt;- list(c(dS, dI, dR ))\n  return(res)\n}\n\node_output_raw &lt;-\n  ode(y = init.values, times = seq(0, tmax), func = SIR_model, parms = parms,\n      method = \"rk4\")\n## Convert to data frame for easy extraction of columns\node_output &lt;- as.data.frame(ode_output_raw)\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment==\"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")"
+    "text": "library(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## determine the type of the next\n        next_event &lt;- sample(x = length(rates), size = 1, prob = rates)\n        ## change to name\n        next_event &lt;- names(rates)[next_event]\n        ## determine type of next event\n        if (next_event == \"infection\") {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (next_event == \"recovery\") {\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\n## Plot the result\nggplot(r) +\n    geom_line(aes(time, I, colour = \"I\"))\n## Re-run this from the line `r &lt;- SIR_gillespie(...)` a few times to \n## convince yourself the output is different every time.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\n## We store the simulations in a data frame, traj, which\n## contains the results from multiple simulation runs and an additional column\n## that represents the simulation index\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\n## Next, plot the multiple simulation runs\nggplot(mlr,\n       aes(x = time, y = value, group = i, color = compartment)) +\n  geom_line() +\n  facet_wrap(~compartment)\n## Plot the distribution of overall outbreak sizes.\n\n## create data frame of outbreak sizes\noutbreak_sizes &lt;- mlr %&gt;%\n  filter(compartment==\"R\") %&gt;%\n  group_by(i) %&gt;%\n  filter(time==max(time)) %&gt;%\n  select(i, size = value)\n\n## plot it as a histogram\nggplot(outbreak_sizes, aes(x = size)) +\n  geom_histogram()\n## determine number of large outbreaks (defined as larger than 50)\noutbreak_sizes %&gt;%\n  filter(size &gt; 50) %&gt;%\n  nrow\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(\n    mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n## We now compare the two averages to the deterministic trajectory\n\n## Model function (from ODE practical)\nlibrary(deSolve)\nSIR_model &lt;- function(times, state, parms){\n  ## Define variables\n  S &lt;- state[\"S\"]\n  I &lt;- state[\"I\"]\n  R &lt;- state[\"R\"]\n  N &lt;- S + I + R\n                                        # Extract parameters\n  beta &lt;- parms[\"beta\"]\n  gamma &lt;- parms[\"gamma\"]\n                                        # Define differential equations\n  dS &lt;- - (beta * S * I) / N\n  dI &lt;- (beta * S * I) / N - gamma * I\n  dR &lt;- gamma * I\n  res &lt;- list(c(dS, dI, dR ))\n  return(res)\n}\n\node_output_raw &lt;-\n  ode(y = init.values, times = seq(0, tmax), func = SIR_model, parms = parms,\n      method = \"rk4\")\n## Convert to data frame for easy extraction of columns\node_output &lt;- as.data.frame(ode_output_raw)\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment==\"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")"
   },
   {
     "objectID": "10_StochasticContinuous_solutions.html#practical-1.-stochastic-simulation-with-the-gillespie-algorithm",
     "href": "10_StochasticContinuous_solutions.html#practical-1.-stochastic-simulation-with-the-gillespie-algorithm",
     "title": "10. Stochastic continuous time models (solutions)",
     "section": "",
-    "text": "library(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## generate cumulative sum of rates, to determine the type of the next\n        ## event\n        cumulative_rates &lt;- cumsum(rates)\n        ## determine type of next event\n        type &lt;- runif(n = 1, min = 0, max = sum(rates))\n        if (type &lt; cumulative_rates[\"infection\"]) {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (type &lt; cumulative_rates[\"recovery\"]){\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\n## Plot the result\nggplot(r) +\n    geom_line(aes(time, I, colour = \"I\"))\n## Re-run this from the line `r &lt;- SIR_gillespie(...)` a few times to \n## convince yourself the output is different every time.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\n## We store the simulations in a data frame, traj, which\n## contains the results from multiple simulation runs and an additional column\n## that represents the simulation index\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\n## Next, plot the multiple simulation runs\nggplot(mlr,\n       aes(x = time, y = value, group = i, color = compartment)) +\n  geom_line() +\n  facet_wrap(~compartment)\n## Plot the distribution of overall outbreak sizes.\n\n## create data frame of outbreak sizes\noutbreak_sizes &lt;- mlr %&gt;%\n  filter(compartment==\"R\") %&gt;%\n  group_by(i) %&gt;%\n  filter(time==max(time)) %&gt;%\n  select(i, size = value)\n\n## plot it as a histogram\nggplot(outbreak_sizes, aes(x = size)) +\n  geom_histogram()\n## determine number of large outbreaks (defined as larger than 50)\noutbreak_sizes %&gt;%\n  filter(size &gt; 50) %&gt;%\n  nrow\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(\n    mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n## We now compare the two averages to the deterministic trajectory\n\n## Model function (from ODE practical)\nlibrary(deSolve)\nSIR_model &lt;- function(times, state, parms){\n  ## Define variables\n  S &lt;- state[\"S\"]\n  I &lt;- state[\"I\"]\n  R &lt;- state[\"R\"]\n  N &lt;- S + I + R\n                                        # Extract parameters\n  beta &lt;- parms[\"beta\"]\n  gamma &lt;- parms[\"gamma\"]\n                                        # Define differential equations\n  dS &lt;- - (beta * S * I) / N\n  dI &lt;- (beta * S * I) / N - gamma * I\n  dR &lt;- gamma * I\n  res &lt;- list(c(dS, dI, dR ))\n  return(res)\n}\n\node_output_raw &lt;-\n  ode(y = init.values, times = seq(0, tmax), func = SIR_model, parms = parms,\n      method = \"rk4\")\n## Convert to data frame for easy extraction of columns\node_output &lt;- as.data.frame(ode_output_raw)\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment==\"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")"
+    "text": "library(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## determine the type of the next\n        next_event &lt;- sample(x = length(rates), size = 1, prob = rates)\n        ## change to name\n        next_event &lt;- names(rates)[next_event]\n        ## determine type of next event\n        if (next_event == \"infection\") {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (next_event == \"recovery\") {\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\n## Plot the result\nggplot(r) +\n    geom_line(aes(time, I, colour = \"I\"))\n## Re-run this from the line `r &lt;- SIR_gillespie(...)` a few times to \n## convince yourself the output is different every time.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\n## We store the simulations in a data frame, traj, which\n## contains the results from multiple simulation runs and an additional column\n## that represents the simulation index\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\n## Next, plot the multiple simulation runs\nggplot(mlr,\n       aes(x = time, y = value, group = i, color = compartment)) +\n  geom_line() +\n  facet_wrap(~compartment)\n## Plot the distribution of overall outbreak sizes.\n\n## create data frame of outbreak sizes\noutbreak_sizes &lt;- mlr %&gt;%\n  filter(compartment==\"R\") %&gt;%\n  group_by(i) %&gt;%\n  filter(time==max(time)) %&gt;%\n  select(i, size = value)\n\n## plot it as a histogram\nggplot(outbreak_sizes, aes(x = size)) +\n  geom_histogram()\n## determine number of large outbreaks (defined as larger than 50)\noutbreak_sizes %&gt;%\n  filter(size &gt; 50) %&gt;%\n  nrow\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(\n    mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment==\"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n## We now compare the two averages to the deterministic trajectory\n\n## Model function (from ODE practical)\nlibrary(deSolve)\nSIR_model &lt;- function(times, state, parms){\n  ## Define variables\n  S &lt;- state[\"S\"]\n  I &lt;- state[\"I\"]\n  R &lt;- state[\"R\"]\n  N &lt;- S + I + R\n                                        # Extract parameters\n  beta &lt;- parms[\"beta\"]\n  gamma &lt;- parms[\"gamma\"]\n                                        # Define differential equations\n  dS &lt;- - (beta * S * I) / N\n  dI &lt;- (beta * S * I) / N - gamma * I\n  dR &lt;- gamma * I\n  res &lt;- list(c(dS, dI, dR ))\n  return(res)\n}\n\node_output_raw &lt;-\n  ode(y = init.values, times = seq(0, tmax), func = SIR_model, parms = parms,\n      method = \"rk4\")\n## Convert to data frame for easy extraction of columns\node_output &lt;- as.data.frame(ode_output_raw)\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment==\"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")"
   },
   {
     "objectID": "10_StochasticContinuous_solutions.html#practical-2.-the-adaptivetau-package",
@@ -347,7 +347,7 @@
     "href": "10_StochasticContinuous_practical.html#practical-1.-stochastic-simulation-with-the-gillespie-algorithm",
     "title": "10. Stochastic continuous models",
     "section": "Practical 1. Stochastic simulation with the Gillespie algorithm",
-    "text": "Practical 1. Stochastic simulation with the Gillespie algorithm\nIn this practical session, we will simulate the SIR model using the Gillespie algorithm.\n\nSetting up the model algorithm\nThe first bit of code below is a function that simulates an SIR model using the Gillespie algorithm.\n\nlibrary(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## generate cumulative sum of rates, to determine the type of the next\n        ## event\n        cumulative_rates &lt;- cumsum(rates)\n        ## determine type of next event\n        type &lt;- runif(n = 1, min = 0, max = sum(rates))\n        if (type &lt; cumulative_rates[\"infection\"]) {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (type &lt; cumulative_rates[\"recovery\"]) {\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\nTake 10 minutes to familiarise yourself with the SIR_gillespie() function and make sure you understand its inner workings.\nNow, run the SIR_gillespie() function using the commands below.\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\nThe data frame r now contains the results of the simulation.\n\n\n\n\n\n\nExercise 1\n\n\n\nInspect the contents of the data frame and write some code to plot the results; you can use plot or ggplot for this. Run the simulation and plot the results multiple times to convince yourself the output is different every time.\n\n\n\n\nSummarising multiple stochastic runs\nNext, we’ll run multiple simulations at once.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\nThe code above runs the simulation 100 times and stores the resulting model trajectories in a data frame, traj, which contains the results from multiple simulation runs and an additional column that represents the simulation index.\n\n\n\n\n\n\nExercise 2\n\n\n\nWrite some code to plot the multiple simulation runs on the same plot.\n\n\nYou’ll notice that some outbreaks die out very quickly, while some others grow to affect large parts of the population.\n\n\n\n\n\n\nExercise 3\n\n\n\nPlot the distribution of overall outbreak sizes. What proportion of outbreaks dies out quickly?\n\n\nNext, we calculate the mean and standard deviation of each state across the multiple runs. To do that, we need the value of the different states at pre-defined time steps, whereas SIR_gillespie only returns the times at which certain events happened. In order to, for example, extract the values of the trajectory at integer time points, we can use the following:\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\nNow, we calculate a summary trajectory containing the mean and standard deviation (sd) of the number of infectious people at every time step:\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment == \"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n\n\n\n\n\n\n\n\nYour plot should look something like the above.\nAs a second summary, we can also consider only trajectories that have not gone extinct, that is, where \\(I&gt;0\\):\n\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment == \"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nExercise 4\n\n\n\nWe have plotted the mean \\(\\pm\\) standard deviation; is this a reasonable summary statistic? What else could one plot? Implement your own plot of trajectory plus uncertainty.\n\n\n\n\nComparing the stochastic model to the equivalent deterministic ODE model\nWe will now compare the stochastic model runs and their averages to the equivalent deterministic ODE model trajectory.\n\n\n\n\n\n\nExercise 5\n\n\n\nUsing code from the ODE practical, load the deSolve package, define the model function (e.g. SIR_model from the ODE practical), run the model with the same parameters as the Gillespie model, and store the output in a data frame called ode_output (remember, you need to convert the matrix returned by deSolve’s ode() to a data frame using as.data.frame()).\n\n\nYou can then plot both the stochastic and deterministic model runs using the following code:\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment == \"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")\n\n\n\n\n\n\n\nExercise 6\n\n\n\nCompare the deterministic trajectory to some of the individual stochastic trajectories and to their mean. Does the deterministic trajectory capture “typical” behaviour of the stochastic trajectories?\n\n\n\n\n\n\n\n\nExercise 7\n\n\n\nRepeat the experiment with different values of the parameters. When does deterministic vs stochastic make more or less of a difference?\n\n\n\n\n\n\n\n\nExercise 8\n\n\n\nRewrite the model code to be an SEIR model. How do the results differ?"
+    "text": "Practical 1. Stochastic simulation with the Gillespie algorithm\nIn this practical session, we will simulate the SIR model using the Gillespie algorithm.\n\nSetting up the model algorithm\nThe first bit of code below is a function that simulates an SIR model using the Gillespie algorithm.\n\nlibrary(ggplot2) ## for plotting\nlibrary(dplyr) ## for manipulation of results\nlibrary(tidyr) ## for storing multiple simulation runs in a data frame\n\n## Function SIR_gillespie.\n## This takes three arguments:\n## - init_state: the initial state\n##   (a named vector containing the number in S, I and R)\n## - parms: the parameters\n##   (a named vector containing the rates beta and gamma)\n## - tf: the end time\nSIR_gillespie &lt;- function(init_state, parms, tf) {\n\n  time &lt;- 0 ## initialise time to 0\n\n  ## assign parameters to easy-access variables\n  beta &lt;- parms[[\"beta\"]]\n  gamma &lt;- parms[[\"gamma\"]]\n\n  ## assign states to easy-access variables\n  S &lt;- init_state[[\"S\"]]\n  I &lt;- init_state[[\"I\"]]\n  R &lt;- init_state[[\"R\"]]\n  N &lt;- S + I + R\n\n  ## create results data frame\n  results_df &lt;- data.frame(time = 0, t(init_state))\n\n  ## loop until end time is reached\n  while (time &lt; tf) {\n    ## update current rates\n    rates &lt;- c(\n      infection = beta * S * I / N, \n      recovery = gamma * I\n    )\n\n    if (sum(rates) &gt; 0) { ## check if any event can happen\n      ## time of next event\n      time &lt;- time + rexp(n = 1, rate = sum(rates))\n      ## check if next event is supposed to happen before end time\n      if (time &lt;= tf) {\n        ## determine the type of the next\n        next_event &lt;- sample(x = length(rates), size = 1, prob = rates)\n        ## change to name\n        next_event &lt;- names(rates)[next_event]\n        ## determine type of next event\n        if (next_event == \"infection\") {\n          ## infection\n          S &lt;- S - 1\n          I &lt;- I + 1\n        } else if (next_event == \"recovery\") {\n          ## recovery\n          I &lt;- I - 1\n          R &lt;- R + 1\n        }\n      } else { ## next event happens after end time\n        time &lt;- tf\n      }\n    } else { ## no event can happen - go straight to end time\n      time &lt;- tf\n    }\n    ## add new row to results data frame\n    results_df &lt;- rbind(results_df, c(time = time, S = S, I = I, R = R))\n  }\n  ## return results data frame\n  return(results_df)\n}\n\nTake 10 minutes to familiarise yourself with the SIR_gillespie() function and make sure you understand its inner workings.\nNow, run the SIR_gillespie() function using the commands below.\n\ninit.values &lt;- c(S = 249, I = 1, R = 0) ## initial state\nparms &lt;- c(beta = 1, gamma = 0.5) ## parameter vector\ntmax &lt;- 20 ## end time\n\n## run Gillespie simulation\nr &lt;- SIR_gillespie(init_state = init.values, parms = parms, tf = tmax)\n\nThe data frame r now contains the results of the simulation.\n\n\n\n\n\n\nExercise 1\n\n\n\nInspect the contents of the data frame and write some code to plot the results; you can use plot or ggplot for this. Run the simulation and plot the results multiple times to convince yourself the output is different every time.\n\n\n\n\nSummarising multiple stochastic runs\nNext, we’ll run multiple simulations at once.\n\n## Run multiple simulation runs and plot a few of them\nnsim &lt;- 100 ## number of trial simulations\n\ntraj &lt;- tibble(i = 1:nsim) %&gt;%\n  rowwise() %&gt;%\n  mutate(trajectory = list(as.data.frame(\n           SIR_gillespie(init.values, parms, tmax)))) %&gt;%\n  unnest(trajectory)\n\n## convert to long data frame\nmlr &lt;- traj %&gt;%\n  gather(compartment, value, 3:ncol(.))\n\nThe code above runs the simulation 100 times and stores the resulting model trajectories in a data frame, traj, which contains the results from multiple simulation runs and an additional column that represents the simulation index.\n\n\n\n\n\n\nExercise 2\n\n\n\nWrite some code to plot the multiple simulation runs on the same plot.\n\n\nYou’ll notice that some outbreaks die out very quickly, while some others grow to affect large parts of the population.\n\n\n\n\n\n\nExercise 3\n\n\n\nPlot the distribution of overall outbreak sizes. What proportion of outbreaks dies out quickly?\n\n\nNext, we calculate the mean and standard deviation of each state across the multiple runs. To do that, we need the value of the different states at pre-defined time steps, whereas SIR_gillespie only returns the times at which certain events happened. In order to, for example, extract the values of the trajectory at integer time points, we can use the following:\n\n## Extract the values of the trajectory at integer time points\ntimeTraj &lt;- mlr %&gt;%\n  group_by(i, compartment) %&gt;%\n  summarise(traj = list(data.frame(\n              time = seq(0, tmax, by = 0.1),\n              value = approx(x = time, y = value, xout = seq(0, tmax, by = 0.1),\n                           method=\"constant\")$y))) %&gt;%\n  unnest(traj)\n\nNow, we calculate a summary trajectory containing the mean and standard deviation (sd) of the number of infectious people at every time step:\n\n## Calculate summary trajectory with mean & sd of infectious people over time\nsumTraj &lt;- timeTraj %&gt;%\n  filter(compartment == \"I\") %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value))\n\n## plot\nggplot(sumTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3)\n\n\n\n\n\n\n\n\nYour plot should look something like the above.\nAs a second summary, we can also consider only trajectories that have not gone extinct, that is, where \\(I&gt;0\\):\n\n## Only consider trajectories that have not gone extinct:\nsumTrajAll &lt;- sumTraj %&gt;%\n  mutate(trajectories=\"all\")\n\nsumTrajGr0 &lt;- timeTraj %&gt;%\n  filter(compartment == \"I\", value &gt; 0) %&gt;%\n  group_by(time) %&gt;%\n  summarise(mean = mean(value),\n            sd = sd(value)) %&gt;%\n  mutate(trajectories=\"greater_than_zero\")\n\niTraj &lt;- bind_rows(sumTrajAll, sumTrajGr0)\n\n## plot\nggplot(iTraj, aes(x = time, y = mean, ymin = pmax(0, mean-sd), ymax = mean+sd,\n                  colour = trajectories, fill = trajectories)) +\n  geom_line() +\n  geom_ribbon(alpha = 0.3) +\n  scale_color_brewer(palette=\"Set1\")\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nExercise 4\n\n\n\nWe have plotted the mean \\(\\pm\\) standard deviation; is this a reasonable summary statistic? What else could one plot? Implement your own plot of trajectory plus uncertainty.\n\n\n\n\nComparing the stochastic model to the equivalent deterministic ODE model\nWe will now compare the stochastic model runs and their averages to the equivalent deterministic ODE model trajectory.\n\n\n\n\n\n\nExercise 5\n\n\n\nUsing code from the ODE practical, load the deSolve package, define the model function (e.g. SIR_model from the ODE practical), run the model with the same parameters as the Gillespie model, and store the output in a data frame called ode_output (remember, you need to convert the matrix returned by deSolve’s ode() to a data frame using as.data.frame()).\n\n\nYou can then plot both the stochastic and deterministic model runs using the following code:\n\n## Combine into one big data frame\nallTraj &lt;- ode_output %&gt;%\n  gather(compartment, value, 2:ncol(.)) %&gt;% ## convert to long format\n  filter(compartment == \"I\") %&gt;%\n  rename(mean = value) %&gt;% ## in deterministic, mean = value\n  mutate(trajectories=\"deterministic\", ## label trajectories\n         sd = 0) %&gt;% ## in deterministic, sd = 0\n  bind_rows(iTraj)\n\n## plot\nggplot(allTraj, aes(x = time, y = mean, colour = trajectories)) +\n  geom_line() +\n  scale_color_brewer(palette=\"Set1\")\n\n\n\n\n\n\n\nExercise 6\n\n\n\nCompare the deterministic trajectory to some of the individual stochastic trajectories and to their mean. Does the deterministic trajectory capture “typical” behaviour of the stochastic trajectories?\n\n\n\n\n\n\n\n\nExercise 7\n\n\n\nRepeat the experiment with different values of the parameters. When does deterministic vs stochastic make more or less of a difference?\n\n\n\n\n\n\n\n\nExercise 8\n\n\n\nRewrite the model code to be an SEIR model. How do the results differ?"
   },
   {
     "objectID": "10_StochasticContinuous_practical.html#practical-2.-the-adaptivetau-package",
