Creating a stochastic SIR model in Julia

Question

I am new to julia and want to create a Stochastic SIR model by following: http://epirecip.es/epicookbook/chapters/sir-stochastic-discretestate-continuoustime/julia

I have written my own interpretation which is nearly the same:

# Following the Gillespie algorthim:
# 1. Initialization of states & parameters
# 2. Monte-carlo step. Random process/step selection.
# 3. Update all states. e.g., I = I + 1 (increase of infected by 1 person). Note: only +/- by 1.
# 4. Repeat until stopTime.

# p - Parameter array: β, ɣ for infected rate and recovered rate, resp. 
# initialState - initial states of S, I, R information.
# stopTime - Total run time. 
using Plots, Distributions

function stochasticSIR(p, initialState, stopTime)

    # Hold the states of S,I,R separately w/ a NamedTuple. See '? NamedTuple' in the REML for details
    # Populate the data storage arrays with the initial data and initialize the run time

    sirData = (dataₛ = [initialState[1]], dataᵢ = [initialState[2]], dataᵣ = [initialState[3]], time = [0]);


    while sirData.time[end] < stopTime

        if sirData.dataᵢ[end] == 0 # If somehow # of infected = 0, break the loop.
            break
        end   

        # Probabilities of each process (infection, recovery). p[1] = β and p[2] = ɣ
        probᵢ = p[1] * sirData.dataₛ[end] * sirData.dataᵢ[end];
        probᵣ = p[2] * sirData.dataᵣ[end];
        probₜ  = probᵢ + probᵣ; # Total reaction rate

        # When the next process happens
        k    = rand(Exponential(1/probₜ));
        push!(sirData.time, sirData.time[end] + k) # time step by k

        # Probability that the reaction is:
        #    probᵢ, probᵣ resp. is: probᵢ / probₜ, probᵣ / probₜ
        randNum = rand();

        # Update the states by randomly picking process (Gillespie algo.)
        if randNum < (probᵢ / probₜ)
            push!(sirData.dataₛ, sirData.dataₛ[end] - 1);
            push!(sirData.dataᵢ, sirData.dataᵢ[end] + 1);
        else
            push!(sirData.dataᵢ, sirData.dataᵢ[end] - 1);
            push!(sirData.dataᵣ, sirData.dataᵣ[end] +1)
        end

    end

end



sirOutput = stochasticSIR([0.0001, 0.05], [999,1,0], 200)
#plot(hcat(sirData.dataₛ, sirData.dataᵢ, sirData.dataᵣ), sirData.time)

Error:

InexactError: Int64(2.508057234147307)

Stacktrace: [1] Int64 at .\float.jl:709 [inlined] [2] convert at .\number.jl:7 [inlined] [3] push! at .\array.jl:868 [inlined] [4] stochasticSIR(::Array{Float64,1}, ::Array{Int64,1}, ::Int64) at .\In[9]:33 [5] top-level scope at In[9]:51

Could someone please explain why I receive this error? It does not tell me what line (I am using Jupyter notebook) and I do not understand it.

David Varela David Varela · Accepted Answer · 2020-02-12T17:51:03

First error

You have to qualify your references to time as sirData.time

The error message is a bit confusing because time is a function in Base as well, so it is automatically in scope.

Second error

You need your data to be represented as Float64, so you have to explictly type your input array:

sirOutput = stochasticSIR([0.0001, 0.05], Float64[999,1,0], 200)

Alternatively, you can create the array with float literals: [999.0,1,0]. If you create an array with only literal integers, Julia will create an integer array.

Creating a stochastic SIR model in Julia

2 Answers

First error

Second error