1、文章编号:1000鄄0887(2013)03鄄0252鄄07 訫应用数学和力学编委会,ISSN 1000鄄0887经经典SIR模型辨识和参数估计问题*何艳辉,摇唐三一(陕西师范大学数学与信息科学学院,西安710062)(我刊编委唐三一来稿)摘要:摇可辨识与否是模型的基本特性,也是研究模型参数估计的基础.在传染病动力学中,SIR模型仍是最常用的模型.该文研究了如何用高阶导数法与多观测点法判定SIR模型的可辨识性.研究表明针对SIR模型的可辨识技巧中多观测点法优于高阶导数法.不仅从理论上判定了SIR模型的可辨识性,而且结合流感疫情数据通过参数估计进一步验证了SIR模型的可辨识性.文中发展的技巧和方
2、法有望推广到其他类型的传染病模型辨识和参数确定上.关摇键摇词:摇可辨识性;摇高阶导数法;摇多观测点法;摇辨识函数;摇参数估计中图分类号:摇 O213摇 摇 摇文献标志码:摇 ADOI: 10. 3879/ j. issn. 1000鄄0887. 2013. 03. 005引摇 摇言日益完善的各种疾病观测数据为辨识性分析提供了保障,也为针对传染病动力学模型的辨识性研究提供了数据支撑.模型可辨识性的直观理解即在给定可观测值与初值时,模型中所有参数是不是有且仅有唯一估计值.对于模型可辨识性严格的数学定义,可参见文献1.为了使模型具有实际应用,理论上应该先从模型可辨识性入手.如果是可辨识的,才能论及模
3、型参数估计与相应的统计推断.然而实际操作中更多的是直接估计模型参数并做统计推断,此时一旦有参数不可辨识,则所做的统计推断是不合理的.因为如果模型中有参数不可辨识,那么它的估计值就不唯一,即该参数不收敛.本文运用高阶导数法2鄄3与多观测点法4,通过构建适当的辨识函数,在先从代数的角度给出了SIR(susceptible鄄infectious鄄recovered)模型的可辨识性条件.然后结合英国流感疫情数据5计算得到所需的辨识等式判定了SIR模型是可辨识的.最后,采用非线性最小二乘法(NLS)来估计SIR模型中可辨识的参数,进一步验证了其可辨识性.需要说明,本文均用圆点数或圆括号里的数表示各阶导数
4、.1摇 SIR模型可辨识性假设:1)一类传染病在某一地区中传播,且所研究地区的种群的总数N是常数;2) S类、I类、 R类种群都具有出生能力且新出生的个体都属于易感者类;3)每一类都具有自然出生率252摇应用数学和力学,第34卷第3期摇 2013年3月15日出版摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇Applied Mathematics and Mechanics摇摇 摇 Vol. 34,No. 3,Mar. 15,2013*收稿日期:摇 2013鄄01鄄15基金项目:摇国家自然科学基金资助项目(11171199)作者简介:摇何艳辉(1987 ),女,江西人,硕士(E鄄ma
5、il:yanhuiqun1988163. com);唐三一(1970 ),男,教授,博士(通讯作者. E鄄mail: sytang snnu. edu. cn).与死亡率.模型所考虑地区中的种群数量N分成3类:易感者类(susceptible, S类)、染病者类(infec鄄tive, I类)和移除者类(removed, R类),即有N = S + I + R,则该SIR模型为摇 摇S = bN - 茁SI - dS,I = 茁SI - 酌I - dI,R = 酌I - dR ,(1)其中,自然出生率与死亡率分别用b,d表示,单位时间内一个病人传染易感者的概率用茁表示,用酌表示单位时间内从染病
6、者中的移出率,并用兹 =(b,茁,d,酌)表示该模型参数向量.在实际生活中,我们容易收集到病人数I的数据,而易感者S和移除者R的数据难以获得,所以I是观测状态变量,S和R是未观测(或潜在)状态变量.代数上系统(1)的可辨识性定义是指对于某个时间t* ,某个正整数k,存在一个函数椎:R4伊 R4(k+1) 寅 R4满足以下条件:摇 摇 det 鄣椎鄣兹 屹 0,摇 摇 椎(兹,I,I, ,I(k) = 0, (2)其中, I, ,I(k)是I在0,t* 上关于时间t的各阶导数,且分别称等式(2)为辨识方程,函数椎( )为辨识函数.由上述定义可知,要分析模型(1)的可辨识性需先构建适当的辨识函数.
7、而且,为了消去其未观测(或潜在)状态变量(即S,R),应对其观测状态变量(即I)求高阶导数.由模型(1)前两个等式易得I的二阶导数(记为I)为摇 摇 I = 茁SI + 茁SI - (酌 + d)I =摇 摇 摇 摇 茁SI + 茁(bN - 茁SI - dS)I - (酌 + d)I =摇 摇 摇 摇 1I (I + 酌I + dI)(I - 茁I2 - dI) + 茁bNI - (酌 + d)I =摇 摇 摇 摇 I2I + 茁bNI - 茁II - 茁(酌 + d)I2 - dI - d(酌 + d)I . (3)显然上式已经不含未观测状态变量S与R了,否则还需求I的更高阶导数.为了方便
8、,我们记c= 茁bN,浊 = 酌 + d,易知兹 = (b,茁,d,酌)与兹* = (c,茁,浊,d)之间存在一个一一映射,则等式(3)右端可表示为摇 摇 f(t,兹* ,I,I) = I2I + cI - 茁II - 茁浊I2 - dI - d浊I . (4)因此,我们可以得到摇 摇鄣f鄣c = I,鄣f鄣茁 =- II - I2浊,鄣f鄣d =- I - I浊,鄣f鄣浊 =- 茁I2 - dI .由上述过程可知,要判断兹*或兹中4个参数的可辨识性,我们需要得到4个可辨识的等式,即辨识性的关键是如何构建辨识函数.接下来,论文将介绍两种如何通过等式(3)和(4)来352何摇 摇艳摇 摇辉摇 摇
9、 摇唐摇 摇三摇 摇一构建辨识函数的方法.1.1摇高阶导数法(higher鄄order derivative method, HODM)该方法的基本思想是通过不断求观测状态变量的更高阶导数,直至消去所有的未观测状态变量,从而得到一系列关于模型参数的辨识函数.根据参考文献中对该方法过程的描述2鄄3,我们可以构建出模型(1)对应的辨识函数:摇 摇 椎0 = (I - f, I(3) - f, , I(5) - f (3)忆 = 0. (5)如果det(鄣椎0 / 鄣兹* ) 屹 0,则由隐函数定理易知兹* = (c,茁,浊,d)可辨识,即如果摇 摇 rank 鄣椎0鄣兹 * = rank鄣f鄣c鄣
10、f鄣茁鄣f鄣d鄣f鄣浊鄣f (1)鄣c鄣f (1)鄣茁鄣f (1)鄣d鄣f (1)鄣浊鄣f (2)鄣c鄣f (2)鄣茁鄣f (2)鄣d鄣f (2)鄣浊鄣f (3)鄣c鄣f (3)鄣茁鄣f (3)鄣d鄣f (3)鄣浊= 4, (6)则兹*就是可辨识的.易知辨识函数(5)中含有I的5阶导数,则计算时至少需要6个I的观测数据.而且需要指出,当兹*为高维参数时,该方法的高阶导数很复杂且式(6)中仍含参数,使得判定模型参数的辨识性很困难,比如通过式(4)计算得到式(6)中的项鄣f (3) / 茁可以表示为摇 摇 鄣f(3)茁 =- 2 (I)2 - 2II(3) - 3II(3) - II(4) - 2
11、浊(3II + II(3),其他项由计算公式可以类似得到.因此,下面将提供另一种计算量相对较少的方法.1.2摇多观测点法(multiple time points method, MTPM)假设已有(I,I,I)在时刻t1, ,t4的值,并用(Ii,Ii,Ii)表示(I,I,I)在t = ti,i = 1, ,4的值.令摇 摇 f1 = f(t1,兹* ,I1,I1), ,f4 = f(t4,兹* ,I4,I4),则由式(3)和式(4)可得摇 摇 椎1 = (I1 - f1, ,I4 - f4)忆 = 0. (7)如果det(鄣椎1 / 鄣兹* ) 屹 0,则由隐函数定理易知兹* = (c,茁
12、,浊,d)可辨识,该条件等价于摇 摇 rank 鄣椎1鄣兹 * = rank鄣f1鄣c鄣f1鄣茁鄣f1鄣d鄣f1鄣浊鄣f2鄣c鄣f2鄣茁鄣f2鄣d鄣f2鄣浊鄣f3鄣c鄣f3鄣茁鄣f3鄣d鄣f3鄣浊鄣f4鄣c鄣f4鄣茁鄣f4鄣d鄣f4鄣浊= 4. (8)同理,通过式(4)计算得到式(8)中的项鄣f3 / 茁可以表示为摇 摇 鄣f3茁 =- I3I3 - I23浊,其他项由计算公式可以类似得到.注意到该方法计算I至少需要3个I的观测数据,那么得到模型(1)对应的辨识函数至少需要6个I的观测数据,与HODM结论一致.452经典SIR模型辨识和参数估计问题2摇 SIR模型可辨识性应用本节将用英国传染病
13、监测中心在1978年发布的有关流感病人的统计数据5来验证SIR模型的可辨识性.该数据统计了两周内每天英格兰北部一所男孩寄宿学校流感爆发和流行情况,总共有763个学生,从发生一个染病者开始统计染病学生的人数,总共统计了15天,每一天统计得到的染病人数分别为摇 摇 摇 摇 1,3,7,25,72,222,282,256,233,189,123,70,25,11,4.由于该数据收集时间很短,学生没有因病死亡,即人口总数是常数,所以我们进一步假设b = d = 0,从而得到相应的简化后SIR模型为摇 摇S =- 茁SI,I = 茁SI - 酌I,R = 酌I .(9)由本文第1节内容可知,要分析模型(
14、9)的辨识性,首先需构建相应的辨识函数.由模型(9)中前两个等式易得I的二阶导数为摇 摇 I = 茁SI + 茁SI - 酌I =摇 摇 摇 摇 茁SI + 茁( - 茁SI)I - 酌I =摇 摇 摇 摇 1I (I + 酌I)I - 茁I2 1I (I + 酌I) - 酌I =摇 摇 摇 摇 I2I - 茁II - 茁酌I2. (10)记滋 = 茁酌,易知棕 = (茁,酌)与棕* = (茁,滋)之间存在一个一一映射,则等式右端可表示为摇 摇 f(t,棕* ,I,I) = I2I - 茁II - 滋I2. (11)因此,我们可以得到摇 摇鄣f鄣茁 =- II,鄣f鄣滋 =- I2 .2.1摇
15、 HODM判定模型(9)可辨识由HODM可以构建出模型(9)对应的辨识函数为椎0 = (I - f,I(3) - f (1)忆 = 0.如果det(鄣椎0 / 鄣棕* ) = I2(I2 - II) 屹 0,则由隐函数定理易知棕*可辨识,即等价为摇 摇 rank 鄣椎0鄣棕 * = rank鄣f鄣茁鄣f鄣滋鄣f (1)鄣茁鄣f (1)鄣滋= rank - II - I2- I2 - II - 2II = 2.又根据参考文献6的研究上式可化为摇 摇 rank 鄣椎0鄣棕 * = rank I 0III I- I I = rankI I- I I = 2. (12)从流感爆发后的观测数据可知,病人
16、数量I和其变化速率即I在不断变化,所以I和I都不为0,显然式(12)成立,因此称模型(9)或棕*可辨识.552何摇 摇艳摇 摇辉摇 摇 摇唐摇 摇三摇 摇一2.2摇 MTPM判定模型(9)可辨识同样地,我们可由MTPM构建出模型(9)对应的辨识函数椎1 = (I1 - f1,I2 - f2)忆 = 0.如果det(鄣椎1 / 鄣棕* ) = I1I1I22 - I21I2I2 屹 0,则由隐函数定理易知棕*可辨识,即摇 摇 rank 鄣椎1鄣棕 * = rank鄣f1鄣茁鄣f1鄣滋鄣f2鄣茁鄣f2鄣滋= rank - I1I1 - I21- I2I2 - I22= 2.同理上式可化为摇 摇 r
17、ank 鄣椎1鄣棕 * = rank I1(I1 - I1) I21I2(I2 - I2) I22= rank I1 00 I2I1 - I1 I1I2 - I2 I2=摇 摇 摇 摇 rank I1 - I1 I1I2 - I2 I2= rank I1 I1I2 I2= 2. (13)同理可知I1,I2和I1,I2不都相等,即有式(13)成立,因此称模型(9)或棕*可辨识.2.3摇参数估计在2. 1节与2. 2节我们分别用了两种可辨识性技巧来研究模型(9)的可辨识性,以及在观测数据零误差时判定可辨识性所需的最少样本数据点.然而实际中不可能达到零误差的理想状态,因此,在实际判定模型可辨识性时,
18、需要更多的样本数据点.本小节通过Monte Carlo(MC)随机模拟得到了足够多的样本观测值,然后采用最常用的NLS来估计参数,结果验证了模型(9)是可辨识的.对于仅有观测变量I的模型(9),其观测模型即摇 摇 I(ti) = I(ti;棕) + 着i,摇 摇 i = 1,2, ,n, (14)其中, I(ti;棕)表示含未知参数棕的模型(9)在ti处的解,I(ti)为对应的观测值,随机误差着i通常设为相互独立同服从0均值,滓2方差的Gauss分布,即着i N(0,滓2) .此时,估计参数用NLS与极大似然法(MLE)效果是等价的.易知,NLS最小化的目标函数为摇 摇 ISSR = 移ni
19、= 1I(ti) - I(ti;棕)2, (15)其中, ISSR表示残差平方和.接下来,用MC方法在Matlab中模拟了模型(9)的观测值I(t) .由之前给出的实际数据5可给定初值(S0,I0,R0) = (762,1,0),又根据已有研究结果5选取初值参数值即(茁,酌) =(0郾 002 2,0. 452 9),再通过MC随机模拟得到I在时间点(单位:d)t = (0,1,2, ,14)处的值I(t;棕),然后在等式(14)的基础上对I(t;棕)加入Gauss随机误差着就得到了样本观测值I(t) .最后我们采用NLS方法利用模拟得到的观测值I(t)重新估计模型(9)中的参数.表1摇 滓
20、= 0. 01时,参数估计结果Table 1摇 Estimation result for 滓 = 0. 01parammeters est. value std. error true value茁 0. 002 2 0. 000 8-5E 0. 002 2酌 0. 452 7 0. 351 9-5E 0. 452 9摇 摇选取滓 =0. 01时,每个时间点模拟100次得到的结果见表1与图1.作为比较,取滓 =100,在每个时间点模拟1 000次得到相应结果见表2与图2.652经典SIR模型辨识和参数估计问题表2摇 滓 = 100时,参数估计结果Table 2摇 Estimation res
21、ult for 滓 = 100parammeters est. value std. error true value茁 0. 002 2 0. 000 2 0. 002 2酌 0. 456 6 0. 054 6 0. 452 9图1摇 滓 = 0. 01时,加入Gauss随机误差的观测值、估计的参数对应的解曲线与实际观测数据对比图Fig. 1摇 Initial values S0, I0 and R0 are known and 滓 = 0. 01 (The solid line represents thesolution with noise, the dashed line repre
22、sents the solution of estimated parameterswithout noise, and the stars represent the actual observed data of the infectious)图2摇 滓 = 100时,加入Gauss随机误差的观测值、估计的参数对应的解曲线与实际观测数据对比图Fig. 2摇 Initial values S0, I0 and R0 are known and 滓 = 100 (The solid line represents thesolution with noise, the dashed line
23、represents the solution of estimated parameterswithout noise, and the stars represent the actual observed data of the infectious)摇 摇表1与表2中信息显示,参数估计值与数据真实值非常接近,即使观测数据数据噪音很强,即滓值很大时,参数估计值同样非常接近参数真实值.从图1与图2可以看出,估计曲线与实际观测数据也拟合得非常好.因此,我们在数值上验证了模型(9)是可辨识的.3摇结果与讨论论文重点探讨了SIR模型的辨识性,通过包括HODM与MTPM在内的方法来构建辨识函752
24、何摇 摇艳摇 摇辉摇 摇 摇唐摇 摇三摇 摇一数,从数学上分析了其辨识性.结果表明两种方法得出的结论一致,即该模型是可辨识的,且它们在分析时所需的最少数据点相同.研究表明针对SIR模型的可辨识技巧中MTPM优于HODM .结合文献5中给出的实际数据,在Matlab中用我们采用MC方法进行数值研究,并通过NLS估计模型参数,所得结论进一步验证了SIR模型的可辨识性.需要指出:当模型的参数向量维数很高时,通过MTPM或HODM来构建辨识函数变得很复杂,判定也相当困难,建议寻求更优的方法,我们将在以后的工作中做进一步的研究与讨论.参考文献(References):1摇 Conte G, Moog C
25、 H, Perdon A M. Nonlinear Control Systems: an Algebraic SettingM.London: Springer, 1999.2摇 Xia X, Moog C H. Identifiability of nonlinear systems with application to HIV/AIDS modelsJ. IEEE Transactions on Automatic Control, 2003, 48(2): 330鄄336.3摇 Jeffrey A M, Xia X. Identifiability of HIV/AIDS Model
26、s(Chaptern). In: Deterministic andStochastic Models of AIDS Epidemics and HIV Infections With InterventionM. Singapore:World Scientific Publishing, 2005.4摇 Wu H, Zhu H, Miao H, Perelson A S. Parameter identifiability and estimation of HIV/AIDSdynamic modelsJ. Bulletin of Mathematical Biology, 2008,
27、70(3): 785鄄799.5摇肖燕妮,周义仓,唐三一.生物数学原理M.西安:西安交通大学出版社, 2012. (XIAO Yan鄄ni, ZHOU Yi鄄cang, TANG San鄄yi. The Principle of BiomathematicsM. Xi爷an: Xi爷an Jiao鄄tong University Press, 2012(in Chinese)6摇 Roger A H, Charles R J.矩阵分析M.杨奇译.机械工业出版社, 1985. (Roger A H,Charles R J. Matrix AnalysisM. YANG Qi, Transl. M
28、echanical Industry Press, 1985. (inChinese)Identification and Parameter Estimationfor Classical SIR ModelHE Yan鄄hui,摇 TANG San鄄yi(School of Mathematics and Information Science, Shaanxi Normal University,Xi爷an 710062, P. R. China)Abstract: Whether a model can be identified is a basic characteristic o
29、f the model before stud鄄ying parameter estimation. Until recently, the classical susceptible鄄infectious鄄recovered (SIR)model is still one of the most commonly used models. In present work the algebraic identifi鄄ability of the SIR model by using high鄄order derivative method (HODM) and multiple timepo
30、ints method (MTPM) was studied. The results indicate that the SIR model can be identifiedif only the infectious was reported, and MTPM is much better than HODM. Using the data ofthe flu, the least square method was adopted to estimate the parameters of the SIR model. Theresult further confirmed that the SIR model was identifiable. The methods developed here couldbe applied to investigate other type models and left those for future studies.Key words: identifiability; HODM; MTPM; identification function; parameter estimation852经典SIR模型辨识和参数估计问题
