时间戳

首页 > 解决方案 > 在 Julia 中创建随机 SIR 模型

julia - 在 Julia 中创建随机 SIR 模型

问题描述

我是 julia 的新手，想通过以下方法创建一个随机 SIR 模型：http: //epirecip.es/epicookbook/chapters/sir-stochastic-discretestate-continuoustime/julia

我写了自己的解释，几乎相同：

# 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)

错误：

不精确错误：Int64（2.508057234147307）

Stacktrace：[1] Int64 at .\float.jl:709 [inlined] [2] convert at .\number.jl:7 [inlined] [3] push！在 .\array.jl:868 [内联] [4] stochasticSIR(::Array{Float64,1}, ::Array{Int64,1}, ::Int64) 在 .\In[9]:33 [5] In[9]:51 的顶级范围

有人可以解释为什么我会收到此错误吗？它没有告诉我是哪一行（我使用的是 Jupyter 笔记本），我不明白。

标签： julia

解决方案

第一个错误

您必须将您的引用限定time为sirData.time

错误消息有点令人困惑，因为time它也是一个函数Base，所以它自动在范围内。

第二个错误

您需要将数据表示为Float64，因此您必须明确键入输入数组：

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

或者，您可以使用浮点文字创建数组：[999.0,1,0]. 如果您创建一个仅包含文字整数的数组，Julia 将创建一个整数数组。

推荐阅读