{ "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": [ "
timeSIR
10.0999.01.00.0
25.379752069687035998.02.00.0
35.893184678545169997.03.00.0
48.471821754633556997.02.01.0
58.611176501322996.03.01.0
616.767837900838497995.04.01.0
" ], "text/plain": [ "6×4 DataFrames.DataFrame\n", "│ Row │ time │ S │ I │ R │\n", "├─────┼─────────┼───────┼─────┼─────┤\n", "│ 1 │ 0.0 │ 999.0 │ 1.0 │ 0.0 │\n", "│ 2 │ 5.37975 │ 998.0 │ 2.0 │ 0.0 │\n", "│ 3 │ 5.89318 │ 997.0 │ 3.0 │ 0.0 │\n", "│ 4 │ 8.47182 │ 997.0 │ 2.0 │ 1.0 │\n", "│ 5 │ 8.61118 │ 996.0 │ 3.0 │ 1.0 │\n", "│ 6 │ 16.7678 │ 995.0 │ 4.0 │ 1.0 │" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "head(sir_out)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Visualisation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "using StatPlots" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "50\n", "\n", "\n", "100\n", "\n", "\n", "150\n", "\n", "\n", "200\n", "\n", "\n", "0\n", "\n", "\n", "250\n", "\n", "\n", "500\n", "\n", "\n", "750\n", "\n", "\n", "1000\n", "\n", "\n", "Time\n", "\n", "\n", "Number\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "S\n", "\n", "\n", "\n", "I\n", "\n", "\n", "\n", "R\n", "\n", "\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@df sir_out plot(:time, [:S :I :R], xlabel=\"Time\",ylabel=\"Number\")" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.3", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }