计算科学 - 使用 SDEint python 包的随机 SIR - 吾爱随笔录

使用 SDEint python 包的随机 SIR

计算科学 Python 数字麻木的随机颂歌

2021-12-23 12:38:10

我想使用 SDEint 包给出以下随机 SIR 模型的数值解（绘图）。即，SDE 系统。

{\begin{cases} d S = - β S I d t - σ S I d W \\ d I = (β S I - γ I) d t + σ S I d W \\ d R = γ I d t \end{cases}

$\begin{cases} dS = -\beta SIdt - \sigma SIdW \\ dI = (\beta SI - \gamma I)dt + \sigma SIdW \\ dR = \gamma Idt \end{cases}$

这是我到目前为止一直在尝试的代码。下面我将解释我的问题。

import matplotlib.pyplot as plt
import numpy as np
import sdeint

tspan = np.linspace(0, 5, 1001)

y0 = np.array([900/1000, 100/1000, 0.,])
S = y0[0]
I = y0[1]
R = y0[2]

def f(y, t):

    f0 = -.80*S*I
    f1 = .80*S*I - 1.1*I
    f2 = 1.1*I
    return np.array([f0, f1, f2,])


def G(y,t):
    return np.array([[-0.9*S*I, 0, 0],
              [0.9*S*I, 0, 0],
              [0,0,0]])

result = sdeint.itoint(f, G, y0, tspan )

plt.plot(result)
plt.show()

$dS$ 和 $dI$ 受到相同噪音的影响 $dW$ ，以前我使用过 np.diag 但那是针对独立的 Wiener 进程但它们在这里不是独立的，它们应该是相同的。也就是说，如果我只看系统的扩散部分，我可以写成

[\begin{matrix} - σ S I \\ σ S I \\ 0 \end{matrix}] d W

$\left[ \begin{matrix} -\sigma SI \\ \sigma SI \\ 0 \end{matrix} \right]dW$

不确定如何将其输入我在我的代码中的 G 定义中使用以下矩阵/数组

[\begin{matrix} - σ S I & 0 & 0 \\ σ S I & 0 & 0 \\ 0 & 0 & 0 \end{matrix}]

$\left[ \begin{matrix} -\sigma SI & 0 & 0 \\ \sigma SI & 0 & 0 \\ 0&0&0 \end{matrix} \right]$

当乘以三个维纳过程的列向量时，将产生我上面想要的结果。

我在运行脚本时注意到的关键问题是，我得到了绘制的路径，特别是 $I$ , 它下降到下面 $0$ 情况并非如此，这些数值解应该保持非负。

我的第二个问题是绘制的路径 $S$ 和 $I$ 是作为彼此的反射出来的，这不应该完全是这种情况，因为它们在各自的定义中是不对称的。

任何输入，或对另一个包或方法的建议表示赞赏。

1个回答

此模型在本教程中使用Julia 的 DifferentialEquations.jl实现。这是该代码的一个版本：

using DifferentialEquations
using StochasticDiffEq
using DiffEqCallbacks
using Random
using SparseArrays
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

A = zeros(3,2)

# Make `g` write the sparse matrix values
function sir_noise!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    ifrac = β*c*I/N*S
    #rfrac = γ*I
    du[1,1] = -sqrt(ifrac)
    du[2,1] = sqrt(ifrac)
    #du[2,2] = -sqrt(rfrac)
    #du[3,2] = sqrt(rfrac)
end;

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[2]
end;
function affect!(integrator)
  integrator.u[2] = 0.0
end;
cb = ContinuousCallback(condition,affect!);

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0]; # S,I,R
p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob_sde = SDEProblem(sir_ode!,sir_noise!,u0,tspan,p,noise_rate_prototype=A);
sol_sde = solve(prob_sde,LambaEM(),callback=cb);

df_sde = DataFrame(sol_sde(t)')
df_sde[!,:t] = t;

@df df_sde plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

整个存储库描述了如何实现不同的流行病模型，从延迟跳跃方程到 SDE，因此它是一个很好的资源。

请注意，这里有两件事：

在化学朗之万方程 (CLE) 中，噪声项不是反应，而是反应的平方根。您的方程式缺少平方根，这是更容易将其拉为负数的一件事。事实上，由于与没有正解的事实类似的原因，您的方程不能保证有u' = -sqrt(u)正解。
这使用了ContinuousCallback更容易施加积极性的机制。

有关展示如何执行 (GPU) 并行性之类的扩展示例，请参阅 QuantEcon 随机微分方程说明中的扩展讨论。

其它你可能感兴趣的问题

上一篇将距离矩阵转换回原始数据下一篇Backward Euler + Quasi Newton(Broyden) 方法无法求解 Van der Pol 方程 (Stiff ODE)