文件交换

图片缩略图

广义SEIR流行病模型(拟合与计算)

version 4.8.8 (13.7 MB) by 大肠Cheynet
具有时间依赖性死亡率和恢复率的扩展SEIR模型的数值实现

48下载

更新2021年2月15日

来自GitHub.

视图版本历史

在github上查看许可证

数值实现了一个具有七状态[2]的广义SEIR模型。除了依赖于函数“lsqcurve”的拟合之外,这个实现是从头开始的。因此,目前的实现可能不同于参考[2]中使用的实现。

这个Matlab实现还包括一些与ref.[2]的主要区别。其中,死亡率和恢复率的表达是时间的分析和经验函数。这种时间依赖性背后的想法是,随着时间的增加,死亡率和恢复率应该趋近于一个常量。如果死亡率保持不变,死亡人数可能会被高估。出生和自然死亡并不是这里的模型。这意味着总人口,包括死亡病例数,保持不变。请注意,ref.[2]是未经同行评审的预印本,我没有足够的资格来判断论文的质量。

本文件包括:

- 用于模拟传染病,恢复和死亡案件的时间历史的函数SEIQRDP.M(其他)
—函数fit_SEIQRDP。m来估计SEIQRDP中使用的8个参数。M在最小二乘意义上。
- 一个示例文件文档.MLX,它呈现数字实现。
—示例文件Example_province_region。mlx使用了约翰·霍普金斯大学为中国湖北省COVID-19流行病收集的数据。
- 一个示例文件example_country.mlx,它使用约翰霍普金斯大学Covid-19 Epidemy [3]的数据收集的数据。
- 由Matteo Secli编写的一个文件“ItalianRegions.mlx”(https://github.com/matteosecli.),我已修改为一个稍微更健壮的拟合。
- 一个示例文件中文进入inovinces.mlx,它说明了Fit_seiqrdp.m如何在用于循环中使用的函数如何安装在不同的中国省份的数据[3]中。
- 一个示例“uncertaintiesissis.mlx”,其说明了拟合有限数据集的危险。
- 一个示例“example_us_cities.mlx”,其说明当“恢复”数据不可用时拟合。
-一个例子模拟emultiplewaves,mlx,说明了适合多个流行波。
- 一个函数getdatacovid,它从[3]读取了约翰霍普金斯大学收集的数据。
- 由Matteo Secli编写的一个函数getdatacovid_ita(https://github.com/matteosecli.),从意大利政府[4]收集意大利COVID-19大流行的最新数据
- getDataCOVID_US函数从[3]收集美国的更新数据
- 一个功能核查时间,绘制拟合和计算的死亡和恢复率(质量检查)
-一个函数getMultipleWaves。m,模拟SEIRQDP模型并将其适合于检测到多个流行波的情况。

如果您需要导入当前提交中未包含的数据,您可以在Matlab工作区中交互导入它们(https://se.mathworks.com/videos/importing-data-from-text-files-interactiony-71076.html.

欢迎任何问题,评论或建议。

参考

[1]https://en.wikipedia.org/wiki/compartmental_models_in_epidemiology#bio-mathemation_deterministic_treatment_of_the_sir_model.

L。[2]Peng,杨,W。,,D,诸葛,C,和香港,L(2020)。基于动态模型的COVID-19在中国的流行分析arXiv预印本arXiv: 2002.06563。

[3]https://github.com/CSSEGISandData/COVID-19

[4]https://github.com/pcm-dpc/covid-19

引用

Cheynet,E.广义SEIR流行病模型(拟合和计算)。Zenodo,2020,Doi:10.5281 / Zenodo.3911854。

查看更多风格

评论和评级(117

Felin Wilta

嗨Cheynet,

理解。谢谢!

大肠Cheynet

嗨felin,

如果最初报告的隔离案例为零,则e,i和q的初始值为为零,因此模拟将在任何时间步骤给出“0”。2020年1月,许多国家的活跃案例数量为零。除了中国和其他一些国家之外,我建议不要在2020年3月之前使用数据。

Felin Wilta

嗨Cheynet,
When I run example_county to simulate the epidemy outbreak based on the fitted parameters, when I use longer Initial and final time (i.e from 22 January 2020 to1st February 2021) I get the values for E,I,Q,R,D to be all zeros. but when I adjust to smaller time period (i.e from January 2021 to February 2021) I get nonzero value. This leads to Comparison of the fitted and real data section in example_country graph to not have the recovered (fitted), active (fitted), deceased (fitted) graph and only have the recovered (reported), active (reported), deceased (reported) graph when values are all zeros for E,I,Q,R,D (longer time period cases). May I know what caused this and how to make the values for E,I,Q,R,D nonzero in longer time period (i.e from 22 January 2020 to1st February 2021)?
期待您的澄清,非常感谢!

Allyson Chong.

嗨Cheynet,
非常感谢你的答案!!在我对这一模型的理解时,你当然很有帮助!再次感谢你!!

大肠Cheynet

嗨,阿廖沙,

抱歉回复晚了。我试着逐一回答以下问题:

1文件example_country.mlx包含来自各个国家的数据,这些数据经常更新。这些数据在这里可以在GitHub上获得:https://github.com/CSSEGISandData/COVID-19。用作“猜测”的参数是执行最小二乘拟合所必需的初始猜测。这些由用户给出作为初始条件,并且可以影响参数的拟合值。单位取决于您正在查看的参数。他们的单位需要确保方程的尺寸均匀性。单位时间是“日”。在文件文档中,在更多详细信息中解释了Lambda和Kappa。

2这些是初始条件。所以它们是"猜测"的价值

3是的,我用的是1小时的时间步长。因为单位时间是“天”,所以我用1/24天作为时间步长。它主要是出于实用目的。我使用了龙格-库塔算法的一个稍微修改的版本,它依赖于参数kappa和lambda在两个时间步之间没有显著变化的假设。如果你使用的时间步长太大,模拟的准确性就会下降,甚至会出现收敛问题。然而,如果你将时间步长增加到12小时,也就是dt = 1/2或甚至dt = 1,也就是一天,代码仍然可以工作得相当好,因为在目前的代码中kappa和lambda随时间平滑地变化

你可以自由地做任何你想做的改变。有几种方法可以使模型适合于数据。我使用的拟合过程不是最简单的方法,因为我先拟合lambda和kappa以限制它们的可能值的范围,然后对参数应用严格的上下限值的一般拟合。一种更直接(但不一定更好)的方法是首先适应lambda和kappa,并且不再允许它们的值更改。然后可以使用第二次fit来检索剩下的参数。

Allyson Chong.

你好再次,Cheynet,
我在这之前就解决了这个问题!我还有几个问题,如果您能帮我澄清我的疑问,我将不胜感激。

1.在文件example_country.mlx和eventized seir模型的部分对实际数据的情况下,我必须输入alpha_guess,beta_guess等的一部分是从国家的数据获得的这些数字我想分析或只是一个受过教育的猜测吗?此外,这些参数使用的单位是什么?例如每天/每周/每1000人口等。
2.在同一节中,有一个'lambda_guess = [0.01,0.001,10];'和'kappa_guess = [0.001,0.001,10];%死亡率”。我能知道这些值是什么意思吗?你是如何得到它们的?
3.在“基于拟合参数模拟流行病爆发”一节中,您使用了1/24的时间步长。
dt = 1/24;我可以知道你是如何选择这个值的,它是如何影响输出的吗?
4.我必须在函数fit_SEIQRDP, SEIQRDP和checkRates的值做任何改变吗?我正在使用你分享的github链接中的确认、死亡和恢复病例的全球时间序列数据。我下载了csv文件并导入了它们,而不是直接从源代码获取。这是我唯一的改变。我需要对我分析的每个国家的参数做任何更改吗?

我仍然是一个初学者,在理解代码方面有些困难,非常感谢您的帮助!!

Allyson Chong.

好的,非常感谢你的帮助,Cheynet!

大肠Cheynet

嗨,阿廖沙,

您看到的消息来自Indr为空的事实。不同地说,这意味着您修改的函数getDataCovid无法正确读取Excel文件。我不知道你如何修改文件,所以我无法更多地帮助。替代方案是您手动导入工作区中的Excel文件,让Matlab创建一个函数来读取Excel文件并重新使用此功能(在次要修改后)。

Allyson Chong.

嗨Cheynet,
我试图改变函数'getDataCOVID.m '中的代码。我没有使用Github存储库中的数据,而是下载了数据并将它们分离到excel文件中。我导入数据并更改代码如下。然后在Example_Country文件中运行这个函数。然后错误'index超过数组元素的数量'出现在下面的行

indR = indR (1);

我该如何解决这个问题?

%%导入数据
Ndays =地板(现在)-datenum(2020、01 22)1;%减去一天,因为数据更新的延迟是24小时
opts = delimitedtextimportoptions(“numvariables”,ndays + 5);
选择。变量名= ["ProvinceState", "CountryRegion", "Lat", "Long", repmat("data",1,Ndays+1)];
选择。变量类型= ["string", "string", repmat("double",1,Ndays+3)];
%指定文件级属性
选择。ExtraColumnsRule =“忽略”;
选择。EmptyLineRule =“读”;
状态={“证实”,“死亡”,“恢复”};
对于II = 1:NUMEL(状态)


如果strcmpi(status {ii},'确认')
tableConfirmed = readtable(“data_confirmed.xlsx”,选择);

{2} elseif strcmpi(地位,“死亡”)
tableDeaths = readtable(“data_deaths.xlsx”,选择);

elsefif strcmpi(status {ii},'恢复')
tableRecovered = readtable(“data_recovered.xlsx”,选择);
其他的
错误(“不明身份”)
结尾
结尾

Time = DateTime(2020,01,22):天(1):DATETIME(DATSTR(楼)),'locale','en_US') - 数据项(1);

Md Nasseef

嗨Cheynet,
谢谢你。我已经安装了最新版本的MATLAB R2020B,现在代码正常工作。

大肠Cheynet

嗨,Md Nasseef,

我认为错误来自您需要至少Matlab R2018来使用可选参数“opts = delimitedtextimportions”的事实。您可以在评论中看到用户的可能解决方案来克服此问题(请参阅4月18日左右的评论)

Md Nasseef

你好Cheynet,
我收到了运行getDataCovid的错误(GetDataCovid(第14行)的错误(第14行)opts = delimitedTextImportoptions(“Numvariables”,NDANES + 5));我正在使用matlab r2017b。你能指导我解决这个问题吗?谢谢你。

大肠Cheynet

代码已更新,但在数据库更改后,还有几个要纠正的小错误。通常,示例现在应该正常工作。

大肠Cheynet

嗨Fayeza,
感谢您的反馈!我打算更新代码以使用户更容易选择拟合的起始和结束时间。这应该在一周内完成

Fayeza Fayeza

嗨Cheynetl,
我明白报告的数据从1月22日开始,但代码从jun开始工作。我们可以对代码进行一些更改来从Say Say开始吗?

徐汇阳

大肠Cheynet

嗨布鲁诺,

向量(a, b, c,……将每个参数的不同猜测分组。这是一种将不同的猜测值分组的紧凑方法。下界和上界是在函数fit_SEIQRDP中设置的,因此用户不需要指定它们。

总共有10个参数,所以不是通过试图同时找到10个最佳值来进行拟合的(至少不是直接)。对病死率(3个参数)和回收率(3个参数)进行了首次拟合。这些参数是用单个测量的速率估计的。这种拟合为参数kappa和lambda提供了狭窄的上下边界。所以在第一步之后,10个参数中的6个几乎是固定的。

布鲁诺•席尔瓦

嗨Cheynet,

谢谢你的回答。你澄清了我的部分疑虑。抱歉,如果我不清楚和的问题,当我定义向量[a, b, c],这意味着什么?例如:'a'是下限,'c'是上限,'b'是速率?或者时间序列是由[下限、中点、上限]来定义的?

大肠Cheynet

嗨布鲁诺,

非常感谢您的反馈!

我不太明白你的问题。所以我试着长话短说:
Lambda_guess和kappa_guess是用来促进拟合的值。假设我们想求函数的最小值。如果用户提供的第一个猜测值接近绝对最小值,算法将适当收敛。如果猜测值不充分,则可能以局部最小值结束。这是最小二乘拟合算法的一个常见问题。

如果您的问题是关于kappa和lambda的“物理”值,则文件文档描述了它们,但我认为文件文档的在线显示中存在一个可视错误.MMLX,因为这些方程式没有出现。

checkRates中的离散值来自于报告的数据(“测量”值)。所以分散性取决于数据源。例如,如果你考虑意大利的数据,离散度就非常小。如果你以法国或比利时为例,分散性会更大。

我恰当地回答了你的问题吗?

布鲁诺•席尔瓦

你好,Cheynet。
首先,祝贺你有趣的工作。我有几个小问题要问:为Lambda_Guess和Kappa_Guess选择vectors是什么意思?
在checkRates测试中,测量值的离散值必须始终呈现类似于example_Country?这些测量值是什么?

大肠Cheynet

嗨拉吉,
Q为活动性病例数。在目前的Matlab提交中,它被称为“隔离的”,因为它假设每个测试呈阳性的病例都是隔离的。这个新的“状态”,表示为Q,可能有两种结果:康复或死亡。请注意,在美国的州和城市中,活跃病例数是/不是可用的,所以我使用的是确诊病例总数减去死亡病例数。虽然不理想,但总比没有好。

Raj基肖尔

你好,Cheynet,我有一个小小的疑问。Q(隔离病例)是否与所有活动性病例相同?还是等于确诊病例总数??

大肠Cheynet

嗨Vikas,
哦,我明白了!谢谢,我没有注意到point(2)的打印错误,我会在下一个版本中更正。

Vikas Sharma.

你好,Cheynet,很高兴看到你考虑我的想法。就(2)而言,正如我之前指出的,在
文档.MLX.illustration where you explain the Numerical solution part.
在任何m文件或实现部分都不存在这种拼写错误。所以结果不会受到这些打字错误的影响。然而,一个
当他看到文档段中的矩阵A.MLX中,当他看到矩阵A.mlx时,不应该混淆。
问候

大肠Cheynet

嗨Vikas,

非常感谢您的反馈!
1)你是对的。我还指出,Kappa_0和Lambda_0应该具有时间的反向的维度。我已经使用适当的定义更新了这些常量的文件文档。
奇怪的是,我仍然看到这个矩阵是7x7。是SEIQRDP函数还是fit_SEIQRDP函数?
你提出了很重要的一点!我花了几个小时研究E0和I0。的确,它们会影响身体的健康。理想情况下,它们也应该是自由参数,但这使得许多未知……我已经更新了示例,稍微修改了您的建议,这似乎与来自中国的数据工作得很好。

Vikas Sharma.

Hello E. Cheynet。许多祝贺精彩的工作。有些错字/建议需要注意:
1. K1和LAMDA1均应具有T ^( - 1)的尺寸,或根据不同的公式进行恢复和死亡率。它是对结果没有任何影响,但数学上,指数术语应该是整体零维度。
2.矩阵A应该是7by7,但它目前显示7by8。
3.当前E0和I0在最终估计方面对印度特别适用于印度的扰动数据。因此,设置初始猜测为5天的孵化时间(在大多数文献中报告)应该对应于以下假设
I0 = 0.1 * (1);%最初感染病例数占报告病例的10%
E0 = 5 *钱数;%最初暴露病例数是最初感染的5倍。
问候

大肠Cheynet

嗨,哈维尔,
这是个好问题。E和我都不知道。所以E0和I0也不知道。如果确诊病例数不是零,那么E0或I0不太可能为零。将它们设为零是错误的,因为这会阻碍模拟。E0和I0的选择是任意的。我选择使用E0和I0等于Q0的值,但也有其他选择。

哈维尔Montoliu

你好大肠Cheynet !很有趣的项目!我是一个初学者,有一个问题,你在选择初始暴露和感染人群的值(E0和I0)时遵循的标准是什么?我看到Q0是初始活动性病例(Q0=确诊(1)-康复(1)-死亡(1)),但为什么将E0和I0设置为初始确诊病例(确诊(1),初始报告病例总数)?
谢谢你!

大肠Cheynet

为了避免误解:这个模型没有机器学习,所以“模型训练”在这里是没有意义的。

大肠Cheynet

嗨Tsardoz,

如果您查看fit_seiqrdp的功能,您将看到拟合算法不同时符合所有参数。例如,对于死亡率,获得第一估计并限制以限制实现局部最小值的风险。对于恢复速率,如果可用,则使用类似的方法。这大大减少了过度装备的问题。可以改善吗?毫无疑问。尽管如此,用户仍然有责任评估配件是否有意义。这部分是“核实条文”功能的目的。

此外,拟合算法依赖于约束参数范围。这就是为什么使用lsqcurvfit而不是nlinfit的原因。最后,评估什么是“现实参数值”几乎没有意义,因为有关哪个参数的信息有限,他们真正的意思,以及数据来源。示例:在此线程中有一个讨论,指出彭等人使用的参数增量。被错误地定义为“隔离时间的倒数”。参数的物理解释是问题,而不是它的价值。

Tsardoz

“如果死亡或恢复率随时间恒定,则会过度装饰。”
恐怕恰恰相反。你用有限的数据尝试和估计的参数越多,就会有越多的模型过度拟合,产生的结果就越差。
目标应该是产生真实的参数值,而不是模型训练集预测的LSE。

大肠Cheynet

嗨Tsardoz,
我忽略了如何计算R0值,这意味着我不知道你是否采用了合适的方法。注意,我(到目前为止)没有包括基本复制数的任何计算。这本身就是一项研究。如果病死率或恢复率随时间不变,就会出现过拟合。明确指出,该模型处理的是随时间变化的恢复和死亡率。合适地使用拟合算法是用户的责任。总之:无用输入,无用输出。

大肠Cheynet

你好,Md Humayun Kabir,
检查函数'getDataCOVID'是否在您的工作区中。Matlab版本不应该影响这个问题。

Md Humayun Kabir

你好e . Cheynet
我的MATLAB是2018b版本。当我运行。mlx文件时,我面临以下问题:
未定义函数或变量'getDataCOVID'。
你能告诉我我在做什么问题吗?
提前谢谢你。

大肠Cheynet

嗨陆克尔,
λ0,λ1,κ0,κ1是通过最小二乘法拟合得到的经验系数。如果可能的话,我使用计算的死亡和康复病例的比率来计算初步估计。最后的估计是利用整个方程组得到的。
随着变量“猜测”表示,这是猜测。它是任意选择的开始。
我没有试图计算基本复制数。在这里我帮不上什么忙。
“tau”变量是一个时间延迟,它改进了随时间变化的采收率的建模。λ0,λ1是可调整的无因次系数。λ0是随着时间增大的回收率。
κ0和κ1也是经验系数。κ0是初始死亡率。κ1隐含地包括关于死亡率的时间延迟的信息。我可能会在以后的版本中用更明确的公式来改变他。

路克鲁尔

你好!
我是初学者,我对你的项目非常感兴趣,这对我有很多帮助。为了更好地了解你的工作,我有几个问题要问。
怎么得到λ0 λ1 κ0 κ1?
你是如何估计初始猜测参数的?
基本再现号码的公式是什么?是否配方β*δ*(1-α)^ t正确?
Lamda_Guess三个变量是什么意思?Kappa_Guess中的两个变量是什么?
如果你能给我一个答复,我将不胜感激

Fryese39.

大肠Cheynet

嗨比尔,
这是个好问题。如果恢复的数量未知,我认为将初始条件设置为恢复E0的数量不为0是一个好主意。我没有试图去看不同的E0值会如何影响被隔离和死亡的计算曲线。你提出的数字可能是现实的。我对这个话题所知不多,无法作出判断。我猜E0的选择会影响参数“lambda”。然而,ode系统是耦合的,因此参数“lambda”涉及多个方程。这意味着E0在拟合过程中只起部分(但可能是重要的?)作用。

比尔·巴德

嗨Etienne,非常感谢您在处理恢复次数未知的情况下的英勇努力。您的例子突出了美国,但我主要有兴趣获得英国的前向预测,我居住的地方(英国也没有公布恢复的数字)。

我在想,对于已恢复的未知情况,您是否尝试过将已恢复的数值设置为死亡数的简单倍数?

我用回收的=死亡* 1.5来试图为英国试过这个;没有其他对输入阵列的调整,并且代码返回非常相似的结果与意大利的延时版本,这具有良好的报告数据和非常合理的曲线适合。

我惊讶地看到,我的英国拟合曲线的恢复/死亡比例为1.5的标量,是时变的,与意大利的一个相似。所以你的解决者似乎有几乎魔法的财产。

无需回复,但对您可能拥有的任何评论感兴趣。

比尔陈

大肠Cheynet

嗨zafar,没有.csv文件下载。函数从GitHub存储库中读取表的内容,并将其存储为CSV。然后将CSV导入本地工作空间,读取和删除。还有替代方法可以使用。

yerlan amanbek.

ZAINULABADIN征服者

需要下载.csv文件

戴尔梁

嗨Cheynet,很抱歉发布脚本。我只是以为你可以为'美国'和'韩国'运行脚本,并看到它有效。您可能对“积极”案件的曲线如何接近高峰时,您可能有兴趣......左右4月30日。有趣的是,许多科学家正在使用所确认的日常变化的数据和预测来预测是一个对称形状的曲线的峰值,而“活跃”追踪的预测曲线与长尾的不对称,我相信更现实。感谢您的意见!

大肠Cheynet

你好戴尔,
我不建议在评论部分发布完整的Matlab脚本。众所周知,美国城市没有“恢复数据”。我对“日本”这个国家没有错。数据正确加载,并通过我上传的最新版本成功地完成了装配。您可能使用了更容易出现bug的旧版本。

戴尔梁

E.Cheynet,运行我发布的代码,我发布了'美国','韩国'和“日本”,你会看到这个问题。'美国'和“韩国的工作,但”日本“没有。当您查看CSV数据文件时,它会出现确认,死亡和恢复的数据,但当您使用位置“日本”时,您的代码不会拿起TableConfired,TableAteAls或TablerEcovered。

大肠Cheynet

你好戴尔,
我不太明白你的意思。如果你谈论Github页面,你不能修改主分支,这是我使用的。

戴尔梁

E.Cheynet,请注意,如果你用我添加的代码运行代码,这些数字需要设置文件夹。也就是说,如果您运行US数据,在Figs文件夹…Figs/US…中必须有一个标记为US的文件子文件夹。或者如果你的位置是韩国,那么必须有一个文件设置为…fig /Korea。

大肠Cheynet

萨米亚,
我无法重现您(使用R2019)得到的错误。检查你是否使用了提交的最新版本。如果您仍然有问题,可以使用Matlab的手动调试工具

Samia Sarothi.

你好e . Cheynet
我的matlab是2019年版本。我面临以下问题:

使用lssqcurvefit时出错(第262行)
函数值和ydata大小不等于。

fit_SEIQRDP错误(第90行)
[多项式系数,~,残余,~,~,~,雅可比矩阵)= lsqcurvefit(@(对位,t) modelFun1 (para t),…

大肠Cheynet

另一种替代方案是使用函数:Websave('dummy.csv',yourfilename);这将读取要在Internet上访问的文件并将其存储为CSV文件。然后,您可以使用传统方法来读取Matlab中的CSV文件。

大肠Cheynet

嗨米Farooq,
检查您的Matlab版本是否足够。我认为至少需要Matlab R2018B来使用“delimitedtextimportoptions”函数。
或者,您可以手动导入您想要的表,并要求Matlab自动为您生成一个函数,您将使用的数据。

m farooq.

m farooq.

你好e . Cheynet

我想跑
[tableConfirmed, tableDeaths tableRecovered,时间]= getDataCOVID ()
得到跟随错误

未定义函数或变量'delimitedTextImportOptions'。

getDataCovid中的错误(第14行)
opts = delimitedtextimportoptions(“numvariables”,ndays + 5);

你能告诉我我在做什么问题吗?
问候

大肠Cheynet

据我所知,速率“delta”的维数是时间的倒数,所以它不应该用百分比表示,而应该用天^(-1)表示。

大肠Cheynet

谢谢迪翁。我在更新的代码中重新定义了变量delta,即传染病病例进入隔离区的速率。

迪翁的话

谢谢E. Cheynet的快速响应和所有努力创建这个奇妙的代码实现。治疗三角洲的速度更有意义。刚刚在论文中和您的代码中的ΔVdelta ^( - 1)明确:在代码中,Delta = 1/8意味着人们以12.5%的速度进入隔离区,而Delta = 1/10则意味着进入10%的速度。因此,增加分母会在逻辑上降低实际检疫数并增加感染。

戴尔梁

E.Cheynet……谢谢你把华盛顿州也包括进来。

大肠Cheynet

嗨dionne,
我同意你的观点,这个定义很奇怪。如果你看看引用[2]的论文,可能会有更好的定义。在他们的论文中,他们没有把它定义为隔离时间的倒数,而是人们进入隔离的速度。这更有意义....

[2] Munz,P.,Hudea,I.,Imad,J.,&Smith,R. J.(2009)。当僵尸攻击!:僵尸感染爆发的数学建模。传染病建模研究进展,4,133-150。

迪翁的话

Delta被定义为“平均隔离时间的倒数”,我用它来表示被感染者隔离时间的倒数,以天为单位。但是,增加德尔塔分母(即增加隔离天数)会导致更多的感染。将delta作为隔离时间,可以增加隔离时间,减少感染,但总体曲线是不现实的。我无法从代码中看出我是否误解了delta的含义(也许它是隔离率,尽管引用的Peng等人的论文将其定义为隔离时间),或者是否存在bug。还有人用过参数吗?

戴尔梁

E.CDHEynet ...在这里是我尝试的......我用examply_us_cities.mlx的新位置语句修改了代码,如下所示:

fprintf(['最近的更新:',datestr(时间(结束)),'\n'])
位置= ",华盛顿";
试一试
indc = find(包含(tableconfirmed.combined_key,location)== 1);
indd =查找(包含(tabledeaths.combined_key,location)== 1);
捕获异常
searchloc = strfind(tableconfirmed.combined_key,location);
(searchLoc indC =找到({}):= = 1);

searchloc = strfind(tableconfirmed.combined_key,location);
(searchLoc indD =找到({}):= = 1);
结尾

这是我的输出:
最新更新日期:04-04-2020
comploy_key.
______________________________

“亚当斯,华盛顿,我们”
“Asotin,华盛顿,美国
“宾顿,华盛顿,美国”



“惠特曼,华盛顿,我们”
“雅吉瓦人、华盛顿、我们”
“离开西澳大利亚,华盛顿,美国”
“未赋值的,华盛顿,美国

错误消息:

矩阵索引超出删除范围。

datetime/parenAssign错误(第58行)
thisdata(rowindices)= [];

Daleexample_Washington(第38行)错误
确认时间(< = minNum) = [];

大肠Cheynet

“example_us_cities.mlx”包含整个状态的一个示例。这里的状态被定义为由与相同状态名称相关联的数据库收集的每个城市的总和。显然,如果我们看看总人口,似乎足够好。

大肠Cheynet

你好戴尔,
这对单个国家来说应该是可能的。然而,你需要数据库的存在和更新,例如在Github。我并没有试图去找。

大肠Cheynet

你好戴尔,
感谢您的信息!
我不确定我理解你的请求,但目前的Matlab提交已经包含了一个美国城市的例子:"Example_US_cities.mlx"

大肠Cheynet

嗨雅各波,

主要原因是脚本删除了数据很少的已确认案例,这对拟合没有意义。根据确诊病例的数量,“确诊”变量有可能变为空。如果(手动)运行两倍于同一行,也可能触发错误,因为向量“Confirmed”的大小发生了变化。在这两种情况下,用户都有责任检查矢量“Confirmed”的大小。

Jacobo白菜

你好,

我不能用example3文件运行canada。它在第40行说:“矩阵索引超出了删除的范围”。

戴尔梁

嗨e.cheynet,我一直在玩你的代码......这太棒了!我现在试图分析个人国家,就像华盛顿州一样。我没有太好的运气。您是否有任何代码,如您发布的多个城市的示例?

大肠Cheynet

对美国一些城市来说,不适合的原因似乎是我修改了回收率,使其保持不变。这不是个好主意。所以我决定回归到时间依赖型的恢复率。

罗纳德曼吕斯兹

尤金·加拉格尔

带有美国示例的新代码在我的Matlab 20a和优化工具箱中运行得很好。新的美国城市现场脚本非常好(有些城市不太适合),但遗憾的是,美国没有关于那些被感染但现在已经恢复的城市的准确数据。中国、意大利和法国正在做些什么来获取美国没有做的恢复数据,这将是一件有趣的事情。

尼克N

你好,谢谢你漂亮的代码,
我在getDatacovid中有错误。m in line: opts = delimitedTextImportOptions

再次感谢。

Rejikumar K

你好,
我是Reji Kumar。
我的matlab比2012年年长。我如何使用它中的代码?任何人都请帮忙。

Ranajoy公司古哈

嗨艾蒂安,
非常感谢你的解释。现在我的疑虑很清楚了。因此,在mlx文件'Confirmed(Reported)'=Confirmed(JHU data) -Recovered(JHU data) -Death(JHU data);从JHU数据,然后匹配quartoned Q(t)。我把您代码中的“已确认(报告)”与JHU数据中的“已确认”数字混淆了。再次感谢你。

Ranajoy公司古哈

Ranajoy公司古哈

嗨艾蒂安,
非常感谢你的解释。我在理解上还有差距。请忽略我的上一篇帖子。这是我的疑问。在mlx文件中生成的图中,确认(报告)、恢复(报告)、死亡(报告)如何与实际JHU数据相关?特别是图中已确认(报告)的数据与JHU中已确认的数据的关系。因为,确认的(报告的)数据点会在图表中某个时间之后下降,但JHU的“确认的”数据从未下降。

大肠Cheynet

嗨Ranajoy,
你提出了一个公平的一点。在John Hopkins University的GitHub存放处,只有三类:康复,确认和死亡。他们对“确认”意味着“已经过测试了”的定义。这就是为什么确认的人数不能减少这个数据库。我猜这个术语“确认(传染性)”不是最好的选择,因为传染病患者不会永远感染。确实,在我上传的示例中,我也使用确认的术语。事实上,确认案件是被隔离的案件,因为我假设在真实情况下,一旦有人“确认”,他们被置于隔离区,这可能是最接近的定义“主动案例=确认-Recovered-Deaths“。

Ranajoy公司古哈

你好,
如果我理解正确,数据是从John Hopkins站点提取的,在John Hopkins github中有四个类别:已确认、死亡、恢复、活动(其中已确认=死亡+恢复+活动;在数据中的任何时间点)。因此,根据约翰·霍普金斯的定义,“确认”永远不会下降,因为它是死亡、恢复和活动的总和。但是,在这个广义SEIR模型中,“确认”图在一段时间后会随着模型而下降,所以当我们匹配模型时,我们不应该使用“活动”数据而不是“确认”数据吗?

大肠Cheynet

嗨,尤金,
约翰霍普金斯大学似乎已经开始在他们的Github页面上上传新的数据集,重点是美国。我会等他们上传完所有的东西,然后再看数据。

所有* MLX文件立即运行。实际数据的绘图仍然显示一些负值,但注意到这些是忽略的。最惊人的图形是ItalianRegions.MLX中的Lombardi Region图表。感染曲线的束缚是显而易见的,但如果被忽略的恢复和死亡状态被忽略了。
由于缺乏恢复数据,您已无法使用SEIR模型将您的数据符合您的遗憾。刚才(4/1下午1点下午4点1日)约翰霍普金斯报告190740美国已确认感染,4127人死亡,7141年恢复。

大肠Cheynet

更新:在更改数据库格式后,已更新意大利地区的数据。如果您尝试运行代码的旧版本,它将无法正常工作。

大肠Cheynet

黑尤金,
我也注意到了这个警告。显然,报告的回收+死亡人数高于一些病例中的“确认”的数量。这很奇怪,但可能是由于报告错误。由于它来自数据库的问题,因此我已经实现了一个快速且肮脏的解决方案。在更新的函数fit_seiqrd中,我只是禁止“实际情况”的数量变为否定并将其设置为零

比尔·巴德

嗨etienne,我修复了我昨天提到的问题,通过去采样在采取“差异”之前进行精细采样的拟合曲线。我现在得到一个非常好的装饰与意大利报告的每日价值观。感谢您的耐心,并希望您获得您的WiFi工作。

大肠Cheynet

你好,
抱歉延迟回复。我暂时丢失了对WiFi(仅移动网)的访问。我试着在不久的将来回复

比尔·巴德

嗨,Etienne,我不想偷懒,但是我昨天才开始使用Matlab(我有IDL的经验)。我希望您能给我一些建议,以解决与您的SEIR建模相关的问题。

我拿了你的榜3并为包括英国在内的几个国家,我居住的国家。然后,我将“Diff”操作员应用于累积阵列(安装和报告)以模拟日情况。

我在代码中的累积报告的阵列和其他来源报告的每日案例之间获得了良好的一致性,所以我的“差异”正常工作。然而,在将“差异”施加到拟合曲线时,它们通过长途方式低估了报告的日常情况。

我追溯这一事实,即拟合曲线比报告值的采样更精确(例如471个元素,而意大利的情况是37个)。

请您给我一些非常一般的建议,如何使Fitted数组与reports数组一致,以便'diff'返回可比较的结果,例如使D数组与Example3中的Deaths一致。

欣赏您可以提供的任何建议,并再次感谢一些伟大的工作

Jacobo白菜

你好,

在用example3文件为加拿大建模时,我得到了一个问题。它在第40行说:“矩阵索引超出了删除的范围”。

尤金·加拉格尔

嗨,我只是在下载最新的代码和数据后ran example4。在今天之前,SEIR模型对所有省份运行很好,但今天中国省份的9个是建模负证实的案件。程序打印警告。

比尔·巴德

嗨,fyi我刚下载了echeynet-seir-aefecf6并在matlab r2020a中运行所有示例,它们都执行精细。

非常感谢和尊重您的一些优秀代码!

大肠Cheynet

嗨,乔丹,
你可以查看time(1)的值。If time(1) >=datetime(2020,04,01),则N= 0。我还注意到,约翰霍普金斯大学的数据库昨天已经更改了,所以示例不会很好地工作。我需要更新示例1到4。

乔丹

在这里不对,因为它返回n = 0,这是正确的吗?谢谢你的帮助!!!
dt = 0.1;%时间步
time1 = datetime(时间(1)):dt: datetime(2020、04 01);
N =元素个数(time1);
t = [0:n-1]。* dt;

尤金·加拉格尔

非常感谢。johns Hopkins网站不会在CSV文件中恢复恢复数据是一种耻辱。让我们希望他们开始掌握恢复的数据,以便您的代码可以使用未来日期运行。johns Hopkins网站首次列出了JH网站上的美国恢复。在他们的互动网页上,今天是他们为美国提供回收的数据(827死,354恢复意大利的第1天:20,499死亡,112,982次恢复):
https://coronavirus.jhu.edu/map.html.

约书亚阿桑摩亚

非常感谢,Cheynet。

大肠Cheynet

嗨约书亚和尤金,
您是对的,URL地址中有一个错字。它现在应该纠正。来自John Hopkins University的GitHub帐户使用更新的数据集,而无需“恢复案例”的数量。这是怜悯,因为这是拟合程序的有价值的信息。因此,我决定仅使用最新的数据到23-03-2020。主要原因是提交的目的只是表明了通用的SEIR模型的配合如何在Matlab中完成。

约书亚阿桑摩亚

请下载新的更新,但我有此错误消息。拜托,我该怎么解决这个问题

使用URLREADWRITE时出错(第92行)
服务器没有找到符合此请求的资源。

urlwrite错误(第52行)
[f,status] = UrlReadWrite(mfilename,catcherrorn,URL,filename,varargin {:});

getDataCOVID错误(第36行)
URLWRITE(文件名,'dummy.csv');

尤金·加拉格尔

etienne。我下载了新文件。示例1运行,但示例2,示例3和示例4所有停止使用以下错误。非常感谢您对此代码的工作。

使用URLREADWRITE时出错(第92行)
服务器没有找到符合此请求的资源。
urlwrite错误(第52行)
[f,status] = UrlReadWrite(mfilename,catcherrorn,URL,filename,varargin {:});
getDataCOVID错误(第36行)
URLWRITE(文件名,'dummy.csv');

大肠Cheynet

示例文件现在应该正常工作。

大肠Cheynet

嗨,尤金,
谢谢你的信息!我看到约翰霍普金斯大学建立了一个新的数据库。然而,它们不再包括“恢复”的数量作为时间的函数。这是奇怪的……今晚我将尝试更新这个示例。

尤金·加拉格尔

etienne,所有4个例子和意大利分析昨晚工作,但是当约翰霍普金斯更新他们的数据时,examply2,examply3和例子4都产生了3-24-2020的时间变量无效的错误(时间向量中的NAT)。意大利示例仍然有效,但它可以从意大利下载其数据。这是错误消息,时间向量的最后一个元素是NAT,而不是有效日期。
使用dateformverify时出错(第18行)
DATESTR将日期号转换为日期向量失败。日期超出范围。
datestr错误(第200行)
S = datef。ormverify (dtnumber dateformstr islocal);
datetime/datestr错误(第801行)
s = datestr (datenum(这),变长度输入宗量{:});

天使苏亚雷斯

你好。

如何解决以下问题:

>> fit_SEIQRDP(Q, R, D, Npop, E0, I0, time, guess, varargin)
尝试以函数的形式执行SCRIPT varargin:

谢谢

大肠Cheynet

嗨,尤金,
谢谢你的消息!我没有注意到这个功能“今天”函数所必需的金融工具箱。在我编写的代码的更新版本中:“楼层(现在)”而不是“Datenum(今天)”

约书亚阿桑摩亚

这些代码非常伟大。谢谢E. Cheynet博士,帮助我们解决这些问题。
请问,你能帮我们怎么画残差吗?

尤金·加拉格尔

这个项目很出色。我将在研究生的概率与统计课程中使用它来演示如何用数据拟合模型的参数。我需要编辑一次才能让程序运行,“today”检索今天的日期只在Mathworks金融工具箱中可用,所以我只是以适当的格式输入了今天的日期。

大肠Cheynet

嗨何塞,
谢谢你的建议。不幸的是,我既没有资格也没有时间来完成这项工作或为之做出贡献。

大肠Cheynet

他雅各波,

由于不同的ePIDEMY位置之间的“延迟”,我认为这不会很好地为世界工作。同样重要的是要记住,人们不应该尝试使用SEIQRDP.M功能进行长期预测。这是由于存在的庞大不确定性,特别是当数据集有限时。

Jacobo白菜

你好,
嗨雅各波,
我已经更新了使用更多Pedagogic描述的示例文件,该变量“猜测”以及所用的数字是什么。
回应:
看来你已经把它修好了,变得更直观了!神奇的谢谢。很好的代码。
___
进一步的问题:
我试图为世界做一个SEIR模型。我使用了您的代码,但可能创建了一个名为“世界”的区域,添加所有病例、死亡等……会有帮助吗?
我通过添加4个。mat文件来解决这个问题,这些文件是从。csv文件中导入的,我在其中添加了区域“world”。
就像一个建议一样,拥有世界地区是一个好主意。

到目前为止,我无法对世界进行适当的预测,曲线并不陡峭。我仍在尝试有几个价值观。

穆里斯

你好!我同意......如果你不介意,虽然速度较慢,我们可以使用这个平台进行协作.....我已经有一些好人工作,所以如果你想有什么希望实施的并希望与其他人并行化,让我知道。

这就是我现在的想法
>创建一个“包装器”围绕“钳子”,允许运行所有行数据(国家和地区)
->输出到.csv或xls的输入和输出从Fit
>将其他指标添加到模拟结果中,例如
->天直到峰值
--->感染的总数
——>总死亡人数
>将数据(可能使用API)与其他人口统计信息(可能来自世界经济论坛或世界卫生组织)合并
>在输出数据中运行一些数据分析,试图找到模式和相关性
->在输入和输出之间运行一个相关矩阵,以评估拟合参数对拟合结果和“后处理”变量的敏感性

你觉得呢?

大肠Cheynet

嗨何塞,
我实际上上传了一个更新的版本,昨天“strfind”替换了“包含”功能,适用于使用较旧的Matlab版本的人替换。我不建议你透露你的电子邮件。你否则会受到重大垃圾邮件。

穆里斯

你好Cheynet,

这就是我解决兼容性问题所做的工作

%%%%%%在2015年以来,自“包含”不会工作以来
% id = strfind (tableRecovered.ProvinceState,“我们”);
%index =查找(不是(cellfun('isempty',ID)))));

对于另一个问题,我只是使用绘图命令而不是原始命令,如下所示:

绘图(时间1,q,'r',time1,r,'c',time1,d,'g',time1,i,'b','linewidth',2);
抓住
情节(时间、Confirmed-Recovered-Deaths‘罗’,时间恢复,“bo”,时间,死亡,“去”,“线宽”,1);

我在巴西领导了一群朋友,在不同的车型上工作,一个是其中之一。如果你能与我们联系,那将是伟大的。我不知道我是否可以在这里发布我的电子邮件,但你可以在Linkdin找到我作为JoséRobertoClarkReis。

谢谢,

乔西

法比安

哦,同时发布。非常感谢Cheynet先生,这太棒了!我是一位科学家在电化学归档,但现在我正在在家里使用我的时间来看看现在的其他重要主题:)

大肠Cheynet

嗨,费边,
我上传了一个新版本,其中不再需要数据。使用“getDataCovid”功能,从上面提到的GitHub存储库直接读取数据。

法比安

嘿,社区,

我可能有个很基本的问题。我从GitHub (.csv文件)获得了最近的疾病统计数据。我如何将数据转移到数据。m文件?
我试图将粘贴附加列复制到Matlab中直接在单个表(恢复,确认,死亡,时间)中,并保存数据内容的更改。
我得到了错误:
使用strfind时出错
单元格必须是字符向量的单元阵列。

我不太明白,因为据我所知,这个函数是用来找到我想要评估的理想国家的。(我没有在示例3的代码中更改任何内容)
当我重新加载原始数据。M文件,代码运行良好。

总的来说:谢谢你的出色工作!

此致,

法比安

大肠Cheynet

嗨雅各波,
我已经更新了使用更多Pedagogic描述的示例文件,该变量“猜测”以及所用的数字是什么。

Jacobo白菜

*我解决了问题* =>我不知道是哪个问题,我只是运行在一个。mlx而不是。m
____
你猜的数字是什么?顺序是什么,我想不出来。
我发现他们是
阿尔法:保护率
贝塔:感染率
伽玛:逆平均潜在时间
三角洲:逆平均隔离时间
Lambda0和lambda1:用于时间依赖性治愈率的系数
Kappa0和Kappa1:在时间依赖性死亡率中使用的系数

Jacobo白菜

我在这里发布了一个请求帮助的帖子:https://la.mathworks.com/matlabcentral/answers/511876-issue-with-seir-model-for-mathlab.
我很期待在我的工作中使用你的软件,我真的希望我能使用它。

我得到了这个错误:

使用seiqrdp错误
输入参数太多。

例子中的错误(第34行)
(年代,E, I, Q, R, D, P] = SEIQRDP(α1,beta1γdelta1,λ₁,Kappa1, Npop, E0,钱数,Q0, R0, D0, t);

Matlab专家

大肠Cheynet

黑,
谢谢你的反馈!
在Matlab 2016b中介绍了“包含”功能,如果您有旧版本,则无法使用它。我认为它是“strfind”的替代方案

对于第二个错误,这也可能是一个兼容性问题,尽管对我来说有点模糊。“DateTime”功能在Matlab 2014B中引入,但可能有一些只能使用更新版本的更新。我确实用Matlab 2019b写了这个函数。或者,您可以使用“DateNum”和功能“DateTick”。

穆里斯

你好,

我正在努力运行非常好的实现。示例1工作正常。示例2,我遇到了至少两个问题。第一个是Funcion“包含”,初步代码。我绕过它并开始使用与湖北省对应的156个价值。它跑了好的。

然后,当我走进去的时候,它停止了:

>>
半机(时间1,q,'r',time1,r,'b',time1,d,'k');
错误使用semilogy
需要一个数字或双可转换参数

我正在使用2015A。这可能是一个兼容性问题吗?

谢谢

MATLAB版本兼容性
用R2020B创建
兼容R2018B及更高版本
平台的兼容性
窗户 macOS Linux

社区宝藏狩猎

在MATLAB中心找到宝藏,并发现社区如何可以帮助你!

开始狩猎!