首页
/ SciML/JumpProcesses.jl 教程:离散随机过程与Gillespie方法详解

SciML/JumpProcesses.jl 教程:离散随机过程与Gillespie方法详解

2025-06-04 04:56:34作者:董宙帆

引言

本文将深入讲解如何使用SciML生态系统中的JumpProcesses.jl包来定义和模拟连续时间跳跃过程(也称为点过程)。在生物学领域,这类模型常被称为随机化学动力学模型或Gillespie模型。我们将通过一个经典的SIR传播模型示例,全面介绍JumpProcesses.jl的核心功能和使用方法。

模型背景:SIR传播动力学

SIR模型是流行病学中描述传播动态的经典模型,它将人群分为三类:

  • 易感者(Susceptible, S):可能被影响的健康人群
  • 感染者(Infected, I):已影响并可传播的群体
  • 康复者(Recovered, R):已恢复并获得免疫力的群体

模型包含两个基本反应过程:

  1. 传播过程:易感者与感染者接触后转变为感染者
  2. 康复过程:感染者康复并获得免疫力

数学建模

在数学上,SIR模型可以用跳跃过程来描述。设S(t)、I(t)、R(t)分别表示t时刻三类人群的数量,模型可表示为:

\begin{aligned}
S(t) &= S(0) - Y_1\left(\int_0^t \beta S(s^-)I(s^-)ds\right) \\
I(t) &= I(0) + Y_1\left(\int_0^t \beta S(s^-)I(s^-)ds\right) - Y_2\left(\int_0^t \nu I(s^-)ds\right) \\
R(t) &= R(0) + Y_2\left(\int_0^t \nu I(s^-)ds\right)
\end{aligned}

其中Y₁和Y₂是独立的单位泊松过程,β是传播率参数,ν是康复率参数。

实现方法一:使用Catalyst构建反应网络

对于化学/生物系统,推荐使用Catalyst.jl来方便地构建模型:

using Catalyst
sir_model = @reaction_network begin
    β, S + I --> 2I  # 传播过程
    ν, I --> R       # 康复过程
end

定义参数和初始条件:

p = (:β => 0.1/1000, :ν => 0.01)
u₀ = [:S => 990, :I => 10, :R => 0]
tspan = (0.0, 250.0)

构建离散问题并转换为跳跃问题:

prob = DiscreteProblem(sir_model, u₀, tspan, p)
jump_prob = JumpProblem(sir_model, prob, Direct())

求解并可视化:

sol = solve(jump_prob, SSAStepper())
plot(sol)

实现方法二:直接使用JumpProcesses底层接口

JumpProcesses提供了三种跳跃类型,按性能从高到低排序:

  1. MassActionJump:性能最佳,适用于质量作用定律的反应
  2. ConstantRateJump:中等性能,适用于速率恒定的跳跃
  3. VariableRateJump:最灵活但性能最低,适用于速率随时间变化的跳跃

使用ConstantRateJump

定义传播和康复两个跳跃过程:

β = 0.1/1000.0
ν = 0.01
p = (β, ν)

# 传播过程
rate1(u, p, t) = p[1] * u[1] * u[2]  # β*S*I
function affect1!(integrator)
    integrator.u[1] -= 1  # S减少1
    integrator.u[2] += 1  # I增加1
end
jump1 = ConstantRateJump(rate1, affect1!)

# 康复过程
rate2(u, p, t) = p[2] * u[2]  # ν*I
function affect2!(integrator)
    integrator.u[2] -= 1  # I减少1
    integrator.u[3] += 1  # R增加1
end
jump2 = ConstantRateJump(rate2, affect2!)

构建并求解问题:

u₀ = [999, 10, 0]
tspan = (0.0, 250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump1, jump2)
sol = solve(jump_prob, SSAStepper())
plot(sol, label=["S(t)" "I(t)" "R(t)"])

性能优化建议

  1. 对于符合质量作用定律的反应,优先使用MassActionJump
  2. 当跳跃速率在跳跃间恒定不变时,使用ConstantRateJump
  3. 只有当速率显式依赖时间时才考虑使用VariableRateJump
  4. 对于VariableRateJump,如果可以提供速率上界,使用有界版本提高性能

扩展应用

JumpProcesses.jl还支持:

  • 混合跳跃过程和连续动态(ODE/SDE)的混合系统
  • 近似模拟的τ-leaping方法(通过RegularJump)
  • 复杂的时间依赖速率函数
  • 各种保存控制以减少内存使用

结论

本文详细介绍了使用JumpProcesses.jl构建和模拟离散随机过程的两种主要方法:通过Catalyst.jl的高级接口和直接使用JumpProcesses的底层接口。通过SIR传播模型的示例,我们展示了如何定义跳跃过程、设置模拟参数以及可视化结果。理解不同跳跃类型的性能特点和适用场景,可以帮助用户在模型复杂性和计算效率之间做出合理的选择。

JumpProcesses.jl作为SciML生态系统的重要组成部分,为复杂随机系统的建模和仿真提供了强大而灵活的工具,特别适合生物化学、流行病学、系统生物学等领域的研究应用。

登录后查看全文
热门项目推荐

热门内容推荐

项目优选

收起
docsdocs
OpenHarmony documentation | OpenHarmony开发者文档
Dockerfile
149
1.95 K
kernelkernel
deepin linux kernel
C
22
6
openHiTLSopenHiTLS
旨在打造算法先进、性能卓越、高效敏捷、安全可靠的密码套件,通过轻量级、可剪裁的软件技术架构满足各行业不同场景的多样化要求,让密码技术应用更简单,同时探索后量子等先进算法创新实践,构建密码前沿技术底座!
C
980
395
ohos_react_nativeohos_react_native
React Native鸿蒙化仓库
C++
192
274
RuoYi-Vue3RuoYi-Vue3
🎉 (RuoYi)官方仓库 基于SpringBoot,Spring Security,JWT,Vue3 & Vite、Element Plus 的前后端分离权限管理系统
Vue
931
555
openGauss-serveropenGauss-server
openGauss kernel ~ openGauss is an open source relational database management system
C++
145
190
nop-entropynop-entropy
Nop Platform 2.0是基于可逆计算理论实现的采用面向语言编程范式的新一代低代码开发平台,包含基于全新原理从零开始研发的GraphQL引擎、ORM引擎、工作流引擎、报表引擎、规则引擎、批处理引引擎等完整设计。nop-entropy是它的后端部分,采用java语言实现,可选择集成Spring框架或者Quarkus框架。中小企业可以免费商用
Java
8
0
金融AI编程实战金融AI编程实战
为非计算机科班出身 (例如财经类高校金融学院) 同学量身定制,新手友好,让学生以亲身实践开源开发的方式,学会使用计算机自动化自己的科研/创新工作。案例以量化投资为主线,涉及 Bash、Python、SQL、BI、AI 等全技术栈,培养面向未来的数智化人才 (如数据工程师、数据分析师、数据科学家、数据决策者、量化投资人)。
Jupyter Notebook
75
66
openHiTLS-examplesopenHiTLS-examples
本仓将为广大高校开发者提供开源实践和创新开发平台,收集和展示openHiTLS示例代码及创新应用,欢迎大家投稿,让全世界看到您的精巧密码实现设计,也让更多人通过您的优秀成果,理解、喜爱上密码技术。
C
65
518
CangjieCommunityCangjieCommunity
为仓颉编程语言开发者打造活跃、开放、高质量的社区环境
Markdown
1.11 K
0