{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Stochastic SIR model (discrete state, continuous time) in R" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "library(reshape2)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "sir <- function(beta, gamma, N, S0, I0, R0, tf) {\n", " time <- 0\n", " S <- S0\n", " I <- I0\n", " R <- R0\n", " ta <- numeric(0)\n", " Sa <- numeric(0)\n", " Ia <- numeric(0)\n", " Ra <- numeric(0)\n", " while (time < tf) {\n", " ta <- c(ta, time)\n", " Sa <- c(Sa, S)\n", " Ia <- c(Ia, I)\n", " Ra <- c(Ra, R)\n", " pf1 <- beta * S * I\n", " pf2 <- gamma * I\n", " pf <- pf1 + pf2\n", " dt <- rexp(1, rate = pf)\n", " time <- time + dt\n", " if (time > tf) {\n", " break\n", " }\n", " ru <- runif(1)\n", " if (ru < (pf1/pf)) {\n", " S <- S - 1\n", " I <- I + 1\n", " } else {\n", " I <- I - 1\n", " R <- R + 1\n", " }\n", " if (I == 0) {\n", " break\n", " }\n", " }\n", " results <- data.frame(time = ta, S = Sa, I = Ia, R = Ra)\n", " return(results)\n", "}" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "set.seed(42)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "sir_out <- sir(0.1/1000,0.05,1000,999,1,0,200)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "if(dim(sir_out)[1]==1){\n", " sir_out <- sir(0.1/1000,0.05,1000,999,1,0,200)\n", "}" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
time | S | I | R |
---|---|---|---|
0.000000 | 999 | 1 | 0 |
1.891201 | 998 | 2 | 0 |
3.470562 | 997 | 3 | 0 |
4.169704 | 997 | 2 | 1 |
8.149657 | 996 | 3 | 1 |
11.145904 | 995 | 4 | 1 |