{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Stochastic SIR model (discrete state, continuous time) in Julia" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using DataFrames\n", "using Distributions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "sir (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function sir(beta,gamma,N,S0,I0,R0,tf)\n", " t = 0\n", " S = S0\n", " I = I0\n", " R = R0\n", " ta=DataArray(Float64,0)\n", " Sa=DataArray(Float64,0)\n", " Ia=DataArray(Float64,0)\n", " Ra=DataArray(Float64,0)\n", " while t < tf\n", " push!(ta,t)\n", " push!(Sa,S)\n", " push!(Ia,I)\n", " push!(Ra,R)\n", " pf1 = beta*S*I\n", " pf2 = gamma*I\n", " pf = pf1+pf2\n", " dt = rand(Exponential(1/pf))\n", " t = t+dt\n", " if t>tf\n", " break\n", " end\n", " ru = rand()\n", " if ru<(pf1/pf)\n", " S=S-1\n", " I=I+1\n", " else\n", " I=I-1\n", " R=R+1\n", " end\n", " end\n", " results = DataFrame()\n", " results[:time] = ta\n", " results[:S] = Sa\n", " results[:I] = Ia\n", " results[:R] = Ra\n", " return(results)\n", "end" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "MersenneTwister(UInt32[0x0000002a], Base.dSFMT.DSFMT_state(Int32[964434469, 1073036706, 1860149520, 1073503458, 1687169063, 1073083486, -399267803, 1072983952, -909620556, 1072836235 … -293054293, 1073002412, -1300127419, 1073642642, 1917177374, -666058738, -337596527, 1830741494, 382, 0]), [1.64879, 1.78639, 1.07348, 1.36027, 1.42523, 1.97645, 1.45162, 1.16015, 1.778, 1.6261 … 1.7593, 1.84751, 1.43425, 1.79251, 1.24761, 1.59121, 1.57693, 1.60592, 1.77807, 1.54728], 382)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "srand(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": [ { "data": { "text/html": [ "
time | S | I | R | |
---|---|---|---|---|
1 | 0.0 | 999.0 | 1.0 | 0.0 |
2 | 5.379752069687035 | 998.0 | 2.0 | 0.0 |
3 | 5.893184678545169 | 997.0 | 3.0 | 0.0 |
4 | 8.471821754633556 | 997.0 | 2.0 | 1.0 |
5 | 8.611176501322 | 996.0 | 3.0 | 1.0 |
6 | 16.767837900838497 | 995.0 | 4.0 | 1.0 |