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):已恢复并获得免疫力的群体
模型包含两个基本反应过程:
- 传播过程:易感者与感染者接触后转变为感染者
- 康复过程:感染者康复并获得免疫力
数学建模
在数学上,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提供了三种跳跃类型,按性能从高到低排序:
- MassActionJump:性能最佳,适用于质量作用定律的反应
- ConstantRateJump:中等性能,适用于速率恒定的跳跃
- 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)"])
性能优化建议
- 对于符合质量作用定律的反应,优先使用MassActionJump
- 当跳跃速率在跳跃间恒定不变时,使用ConstantRateJump
- 只有当速率显式依赖时间时才考虑使用VariableRateJump
- 对于VariableRateJump,如果可以提供速率上界,使用有界版本提高性能
扩展应用
JumpProcesses.jl还支持:
- 混合跳跃过程和连续动态(ODE/SDE)的混合系统
- 近似模拟的τ-leaping方法(通过RegularJump)
- 复杂的时间依赖速率函数
- 各种保存控制以减少内存使用
结论
本文详细介绍了使用JumpProcesses.jl构建和模拟离散随机过程的两种主要方法:通过Catalyst.jl的高级接口和直接使用JumpProcesses的底层接口。通过SIR传播模型的示例,我们展示了如何定义跳跃过程、设置模拟参数以及可视化结果。理解不同跳跃类型的性能特点和适用场景,可以帮助用户在模型复杂性和计算效率之间做出合理的选择。
JumpProcesses.jl作为SciML生态系统的重要组成部分,为复杂随机系统的建模和仿真提供了强大而灵活的工具,特别适合生物化学、流行病学、系统生物学等领域的研究应用。
登录后查看全文
热门项目推荐
- QQwen3-Omni-30B-A3B-InstructQwen3-Omni是多语言全模态模型,原生支持文本、图像、音视频输入,并实时生成语音。00
- HHunyuan-MT-7B腾讯混元翻译模型主要支持33种语言间的互译,包括中国五种少数民族语言。00
GitCode-文心大模型-智源研究院AI应用开发大赛
GitCode&文心大模型&智源研究院强强联合,发起的AI应用开发大赛;总奖池8W,单人最高可得价值3W奖励。快来参加吧~0269get_jobs
💼【AI找工作助手】全平台自动投简历脚本:(boss、前程无忧、猎聘、拉勾、智联招聘)Java00AudioFly
AudioFly是一款基于LDM架构的文本转音频生成模型。它能生成采样率为44.1 kHz的高保真音频,且与文本提示高度一致,适用于音效、音乐及多事件音频合成等任务。Python00GOT-OCR-2.0-hf
阶跃星辰StepFun推出的GOT-OCR-2.0-hf是一款强大的多语言OCR开源模型,支持从普通文档到复杂场景的文字识别。它能精准处理表格、图表、数学公式、几何图形甚至乐谱等特殊内容,输出结果可通过第三方工具渲染成多种格式。模型支持1024×1024高分辨率输入,具备多页批量处理、动态分块识别和交互式区域选择等创新功能,用户可通过坐标或颜色指定识别区域。基于Apache 2.0协议开源,提供Hugging Face演示和完整代码,适用于学术研究到工业应用的广泛场景,为OCR领域带来突破性解决方案。00- HHowToCook程序员在家做饭方法指南。Programmer's guide about how to cook at home (Chinese only).Dockerfile08
- PpathwayPathway is an open framework for high-throughput and low-latency real-time data processing.Python00
热门内容推荐
最新内容推荐
PCDViewer-4.9.0-Ubuntu20.04:专业点云可视化与编辑工具全面解析 JDK 8u381 Windows x64 安装包:企业级Java开发环境的完美选择 QT连接阿里云MySQL数据库完整指南:从环境配置到问题解决 SAP S4HANA物料管理资源全面解析:从入门到精通的完整指南 VSdebugChkMatch.exe:专业PDB签名匹配工具全面解析与使用指南 OpenSSL 3.3.0资源下载指南:新一代加密库的全面解析与部署教程 基恩士LJ-X8000A开发版SDK样本程序全面指南 - 工业激光轮廓仪开发利器 ZLIB 1.3 静态库 Windows x64 版本:高效数据压缩解决方案完全指南 SteamVR 1.2.3 Unity插件:兼容Unity 2019及更低版本的VR开发终极解决方案 全球36个生物多样性热点地区KML矢量图资源详解与应用指南
项目优选
收起

OpenHarmony documentation | OpenHarmony开发者文档
Dockerfile
149
1.95 K

deepin linux kernel
C
22
6

旨在打造算法先进、性能卓越、高效敏捷、安全可靠的密码套件,通过轻量级、可剪裁的软件技术架构满足各行业不同场景的多样化要求,让密码技术应用更简单,同时探索后量子等先进算法创新实践,构建密码前沿技术底座!
C
980
395

React Native鸿蒙化仓库
C++
192
274

🎉 (RuoYi)官方仓库 基于SpringBoot,Spring Security,JWT,Vue3 & Vite、Element Plus 的前后端分离权限管理系统
Vue
931
555

openGauss kernel ~ openGauss is an open source relational database management system
C++
145
190

Nop Platform 2.0是基于可逆计算理论实现的采用面向语言编程范式的新一代低代码开发平台,包含基于全新原理从零开始研发的GraphQL引擎、ORM引擎、工作流引擎、报表引擎、规则引擎、批处理引引擎等完整设计。nop-entropy是它的后端部分,采用java语言实现,可选择集成Spring框架或者Quarkus框架。中小企业可以免费商用
Java
8
0

为非计算机科班出身 (例如财经类高校金融学院) 同学量身定制,新手友好,让学生以亲身实践开源开发的方式,学会使用计算机自动化自己的科研/创新工作。案例以量化投资为主线,涉及 Bash、Python、SQL、BI、AI 等全技术栈,培养面向未来的数智化人才 (如数据工程师、数据分析师、数据科学家、数据决策者、量化投资人)。
Jupyter Notebook
75
66

本仓将为广大高校开发者提供开源实践和创新开发平台,收集和展示openHiTLS示例代码及创新应用,欢迎大家投稿,让全世界看到您的精巧密码实现设计,也让更多人通过您的优秀成果,理解、喜爱上密码技术。
C
65
518

为仓颉编程语言开发者打造活跃、开放、高质量的社区环境
Markdown
1.11 K
0