（应用数学专业论文）两种群食饵—捕食者模型的参数估计的研究.pdf

摘要 摘要 经过一个世纪的发展，生物数学模型的研究得到了广泛的应用，同时也产 生了微分方程的参数估计问题。这是当前统计领域中的热点研究课题。目前利 用观测数据来估计常微分方程参数的方法计算较复杂，与实际数据吻合效果较 差。因此，研究高精度的参数估计方法显得越发重要。 近来，r a m s a y 等人用样条函数作为基函数，给出了一类食饵一捕食者模型 的参数估计，但他们没有充分利用该类模型已有的大量定性性质。本文引进一 种两阶段最小二乘估计法。与以前方法不同的是由于来源于食饵一捕食者模型 的观测数据具有明显的周期性，因此，我们在第一阶段取三角函数作为基函数， 构造了观测数据表示的函数的潜周期模型( h p m ) 估计式，并给出了估计的误差表 达式，结果表明h p m 估计具有良好的性质和数值稳定性。然后，在第二阶段， 通过变换，将食饵一捕食者模型的参数估计问题转化为线性最小二乘问题。这 样得到的估计具有良好的统计性质，并具有计算速度快，计算精度高等优点。 最后，本文针对两种群的食饵一捕食者模型的参数估计给出了数值模拟，通过 统计检验、对比估计值与真实值，并与有关文献中的结果比较，本文提出的两 阶段最小二乘估计法优于现有的结果。 本文由五章构成。第一章简述了问题产生的历史背景、本文的主要工作以 及本文中需要用到的一些定义和引理。在第二章中，我们详细地讨论了周期型 随机数据的数值微分及其性质。在第三章中，讨论了两种群食饵一捕食者模型 的参数估计。在第四章中，我们根据真实参数产生了一组观测数据，然后利用 这组数据，利用本文的方法估计出模型参数，结果优于相关文献中给出的结果， 说明本文的方法是可行的。在第五章中，对全文作了总结，并给出了研究展望。 关键词：食饵一捕食者模型：参数估计；潜周期模型 a b s t r a c t 一一 a b s t r a c t t h er e s e a r c ho fm a t h e m a t i c a lb i o l o g ym o d e lh a sg o tt h ee x t e n s i v ea p p l i c a t i o n t h r o u g ht h ed e v e l o p m e n to fac e n t u r y a tt h es a m et i m e ，e s t i m a t i o no ft h eo r d i n a r y d i f f e r e n t i a le q u a t i o n s ( o d e ) h a s a p p e a r e d p r e s e n t l y , m a n ym e t h o d sf o re s t i m a t i n g p a r a m e t e r sf r o mn o i s yd a t aa r ec o m p u t a t i o n a l l yi n t e n s i v ea n do f t e nf i tp o o r l yw i t ht h e r e a ld a t a c o n s e q u e n t l y , s t u d yf o ram e t h o dw i t hh i g hp r e c i s i o n b e c o m e sm o r ea n d m o r ei m p o r t a n t l a t e l y , r a m s a ye r e h a v eg i v e na l le s t i m a t i o nm e t h o do fp a r a m e t e r si nt h e 2 - s p e c i e sp r e d a t o r - p r e ym o d e lb yt r e a t i n gt h es p i n ef u n c t i o n sa st h eb a s i ss y s t e m h o w e v e r , t h e yd i d n tm a k em o s tu s eo fq u a l i t a t i v ep r o p e r t yo ft h em o d e ls u c h 勰t h e p e r i o d i c i t y i nt h i sp a p e r , w ei n t r o d u c ean e wt h ee s t i m a t i o nm e t h o d 铆o s t a g e e s t i m a t i o nf f s e ) t os t u d yt h ee s t i m a t i o no fp a r a m e t e r si nt h e2 - s p e c i e s p r e d a t o r - p r e v m o d e l f i r s t ，石o ) i ss h o w e da st h eh i d d e np e r i o d i cm o d e l ( h p m ) b yt r e a t i n gt h e t r i g o n o m e t r i cf u n c t i o na st h eb a s i ss y s t e m u s i n gt h ee s t i m a t i o nm e t h o do ft h e h i d d e np e r i o d i cm o d e l ( h p m ) ，w ee s t i m a t e k ib yo b s e r v i n gt h ed i a g r a mo f k ( a ja n de s t i m a t et h ea i l g l ef r e q u e n c y a n dt h eo r i g i n a lp h a s e n e x t ，t h e p r o b l e mi st r a n s f o r m e dt h el e a s ts q u a r ep r o b l e ma n do t h e rp r o b l e m sa r ee s t i m a t e d f i n a l l y , c o m b i n e dw i t hr e a ld a t a ，n u m e r i cs i m u l a t i o nh a sb e e nd o n ea n ds h o w st h a t t h em e t h o di st r u l yf e a s i b l e t h i sp a p e ri sc o m p o s e do ff v e p a r t s i nt h ef i r s tc h a p t e r , w ei n t r o d u c et h e h i s t o r i c a lb a c k g r o u n do ft h ep r o b l e m sw h i c hw i l lb ei n v e s t i g a t e d ，t h em a i nr c s u i t s a n ds o m ed e f i n i t i o n sa n dl e m m a sw h i c ha r er e q u i r e di n t h i s p a p e r i nt h es e c o n d c h a p t e r , w ed i s c u s sn u m e r i c a ld i f f e r e n t i a t i o na n di t sp r o p e r t yo fs t o c h a s t i cd a t u m a n di nt h et h i r dc h a p t e r , w ed i s c u s st h ee s t i m a t i o nm e t h o do fe a c hp a r a m e t e ri l lt h e 2 - s p e c i e sp r e d a t o r 。p r e ym o d e li nd e t a i l a ti nt h ef o r t hc h a p t e r , w ep r o d u c ea g r o u po f o b s e r v a t i o nd a t a b yu s i n gt h er e a lv a l u eo ft h ep a r a m e t e r a c c o r d i n gt ot h e o b s e r v a t i o nd a t aw eu s et h em e t h o di n t r o d u c e di n t h i s p a p e rt o e s t i m a t et h e p a r a m e t e r s ，w h i c hs h o w st h er e s u l t sa r eb e t t e rt h a nt h o s eg i v e ni np e r t i n e n tr e f e r e n c e s a b s t r a c t s ot h em e t h o di n t h i sp a p e ri sf e a s i b l e a tl a s t ，w ec o n c l u d et h ep a p e ra n dg i v et h e r e s e a r c hf o r e c a s t k e yw o r d s ：p r e d a t o r - p r e ym o d e l ；p a r a m e t e re s t i m a t i o n ；h i d d e np e r i o d i cm o d e l i i l 学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工 作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地 方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含 为获得直昌太堂或其他教育机构的学位或证书而使用过的材料。与 我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确 的说明并表示谢意。 学位论文储签名c 瑚：荆 学位论文版权使用授权书 瑚彳日 本学位论文作者完全了解直昌太堂有关保留、使用学位论文 的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁 盘，允许论文被查阅和借阅。本人授权直昌太堂可以将学位论文的全 部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描 等复制手段保存、汇编本学位论文。同时授权中国科学技术信息研究 所将本学位论文收录到中国学位论文全文数据库，并通过网络向 社会公众提供信息服务。 ( 保密的学位论文在解密后适用本授权书) 导师签名：8 靳生 签字醐：川年堋1 日 篼l 章引言 1 1 选题背景 第1 章引言 本课题来源于导师的省自然科学基金项目( 2 0 0 7 g z s 2 3 9 8 ) 。 微分方程模型参数估计是目前国际上概率统计领域的研究热点i 卜引。微分方 程模型是描述实际对象的某些特性随时问或空间演变的过程，进而分析它的变 化规律、预测它的未来性态、研究它的控制手段的重要方法。但是，诸如生物 科学、物理学等领域中建立的微分方程通常没有解析解，实际应用中一般使用 它们的数值解。数值解的精度直接依赖于方程中的各参数以及初始值( 或边界值) 的精度。h i m m e l b l a u 6 j 在1 9 6 7 年首次利用统计方法进行了微分方程参数估计的 研究，b o c k l7 j 于1 9 8 3 年综述了自2 0 世纪6 0 年代以后2 0 年间常微分方程参数辨 识的研究进展，指出现有的利用观测数据估计非线性微分方程模型参数的方法 效果较差，如计算速度慢、拟合精度差、得不到有较好统计解释的区间估计等， 无法满足实际工作的需要。j e r o m e l 4 j 在2 0 0 5 年直接考虑求解参数的非线性最小 二乘估计，利用信赖域方法求解非线性优化问题，其结果虽然较为精确，但计 算过程十分复杂。2 0 0 6 年j i g u oc a o 2 1 和2 0 0 7 年r a m s a y l 4 j 采用函数型数据分析 法( f d a ) 研究了一类微分方程的参数估计方法。这些方法主要存在三点不足，一 是，在计算过程中需要数值求解微分方程；二是，没有充分利用微分方程的大 量定性研究结果，如其相轨线的周期性【8 , 1 1 1 等；三是，不能得出具有统计意义的 置信区间，这在参数估计理论中具有十分重要的作用【9 1 。 常微分方程( o r d i n a r yd i f f e r e n t i a le q u a t i o n ，o d e ) 广泛应用于工程力学、生物 科学、物理学、化学工程学、生理学等领圳1 0 , 1 1 】，所以对微分方程的参数估计问 题的研究具有十分重要的实际意义。假设所研究的系统是封闭的，建立的o d e 模型一般可表示为1 1 2 j ： x 0 ) = 厂b ，t1 0 )( 1 1 ) 一 ，、 其中x o ) = g 。o ) ，j d g ) ) ，立g ) ；a x _ q ) ，口为待定参数向量，f 0 ，丁 。 口f 若外界向系统有输入量“i f ) ，o d e 模型可表示为 第l 章t j i ，j 二= f ( x ，“，t10 )( 1 2 ) 上面两式同样可以处理高阶导数d ”x 的情况。我们知道诸如生物科学、物 理学等领域中建立的o d e 直接依赖j 二方年罕l ，的各参数以及初始值( 或边界值) 的 精度。但是微分方程通常没有解析解，实勋0 迎用中一般使用它们的数值解。求 数值解的方法一般采用r u n g e k u t t a 迭代法，参数估计方法一般采用非线性优化 方法，如最小二乘法( n o n l i n e a rl e a s ts q u a r e s ，n l s ) 等。但这些方法都存在很多缺 陷：第一，迭代过程中误差不断累加、计算遣大、稳定性差、过分依赖初始值； 第二，如果初始值未知，需要将初始值视为参数的形式进行估计，这对估计的 算法要求很高；第三，很多非线性优化算法足局部收敛或局部稳定的，初始值 的选取要充分靠近真实解，否则可能得不到理想的结剁2 , 1 3 】。所以，根据我们观 测到的离散的噪声数据，探讨新的参数估计的方法势在必行。b o c k l 7 l 、 h i m m e l b l a u 6 1 、d eb o o r 与s w a r t z l l 引、g e l m a n 1 5 j 提出了许多新的方法。2 0 0 6 年， r a m s a y ，h o o k e r ，g 和j i g u oc a o 3 l 【1 6 1 提出了一种更有效的方法一广义轮廓法 ( g e n e r a l i z e dp r o f i l i n gm e t h o d ，简称g p m ) ，用g p m 估计出的微分方程模型的解 与实际数据吻合的精度相当高，误差在1 一5 之间，可以基本满足高精度参数 估计的要求。在g p m 中，r a m s a y 等应用了参数级联( p a r a m e t e rc a s c a d e s ) 的处理 方法。例如，考虑一维的情形，设y = ( y 1 一，y 。) 表示x o ) 在时刻l ，一，t 。时的n 个 观测数据。g p m 中以样条函数为基函数，将x ( f ) 线性表示成 x x ( f ) = y c k 仇0 ) = c g )( 1 3 ) 上式中的k 并不是固定的，它随着t 的增大而增大。k 值的确定称为g p m 的定阶问题。系数向量c 称为冗余参数( n u i s a n c ep a r a m e t e r s ) 或局部参数( 1 0 c a l p a r a m e t e r s ) 。设观测误差向量e = ( q ， 乌) i = ( y 。一z 。) ，y 。一x o 。) ) ，且p 服从某 基于仃的联合分布密度函数g ( ei 仃) ( 如实际中，通常假定e 服从期望为o 、标准 差为口的正态分布) 。如果没有时间l o ，t l 端点处的观测数据，或者提供的观测数 据有误差，则初始值、边界值也作为待估计参数，从而得到增广参数向量 0 = b ( o y ，盯，0 j ，0 称为结构参数( s t r u c t u r a lp a r a m e t e r s ) 或全局参数( g l o b a l p a r a m e t e r s ) 。 g p m 法的优点是，光滑度较好，计算快速。不足的是，没有充分 利用方程( 1 1 ) 的大量定性研究成果，如其相轨线的周期性【8 , 1 1 】。 2 第l 章引蠢 另一方面，许多统计学家应用潜周划模犁( h p m ) 参数估计方法对具有周期性 问题进行了研究【1 7 2 5 1 。例如，一维简单的潜嗣期模型一余弦波信号可表示为 x 。一套印。s b p 驴小鼬+ ( 1 4 ) w h i t t l e l l 7 , 1 8 1 提出了基于线性过程= 二。口t m 的周期图和白噪声f ，的周期图 之间的渐进关系式的一种潜频率检验方法。但是w h i t t l e 的方法在潜频率的个数 q 未知i f j ，只能检测出y ，的谱密度厂，( a ) 的估计五以) ，而夕，n ) 在潜频率附近会 出现虚假峰，使检验的功效大大降低，这样就导致潜频率不易被检测出。为了 克服上述困难人们提出了不少新的方法，其中最值得注意的是h a n n a n 提出的方 法。h 狮n a n f 2 0 】通过y 的加窗周期图，给出了以- - 1 中新的检验统计量。然而，h a n n a n 的方法仍然有局限的，它只适用于潜频率个数已知的情形。对于潜频率个数未 知的情形，何书元【2 1 】通过周期图，利用a 一区间分析法给出了潜频率的个数和 潜频率的强相合估计。 在自然界中，许多种群竞争同一有限食物，从而抑制了相互问种群的增长， 这是一个有趣而普遍的现象。于是研究与实际相符的多种群竞争模型的参数估 计方法是很重要的【1 4 ，一7 1 。l o t k a v o l t e r r a 竞争系统是著名的种群模型之一【2 6 1 ，以 两种群生物为例，记a 生物是捕食者，b 生物是被捕食者，假设t 时刻捕食者a 的密度为x ( f ) ，被捕食者b 密度为j ，( f ) ，用下面的方程描述一般的食饵一捕食者 系统 x g ) = x ( f 壮，+ a 2 y ( t ) 1 y ( f ) = y ( t l a ，+ 口。x ( f ) 】 ( 1 5 ) 初始条件为 i x ：竺 。口s ( 1 6 ) 1y ( o ) ；口。 u 朋， 其中a 。，( 1 s 口s6 ) 为模型的待定参数。 与( 1 1 ) 式对应，有 x 。o ) = z o ) ，x ：o ) = y o ) 石( f ) = g ( f ) ，y ( f ) ) ；g ，( f ) ，x ：( f ) ) ，工o p ) = ( 口，口。) 卢一k 。，c 1 6 ) ，为结构参量， ( x ，ri ) = 仁。( f 壮。+ 口：x ：( t ) l x 2 ( f ) k ，+ c ? 1 4 x 。( f d 3 第1 章l j l0 毒f 1 外界向两种群食饵一捕食者系统自输入专 ：“( f ) ( 已知) ，则考虑厂b ，“，tl 卢) 。 从i 刷也食饵一捕食者系统( 1 5 ) 和( 1 6 ) 转化成阳嚣的形式 二o ) = ，g ，tl ) 或x ；r ( x ，蹦，ti ) 1 2 主要工作 ( 1 7 ) ( 1 8 ) 本文旨在研究常微分方程( 1 7 ) 的参数伊i 汁，得到高精度的结构参数。 首先，讨论周期型随机数据的数值微分及其性质。文中给出了周期型数据 的三角最佳逼近多项式，并讨论了其数值微分的估计误差。这是下面进行两阶 段最小二乘估计的理论基础。 其次，用两阶段最小二乘估计法来估计参数罗，与以前方法不同的是充分 利用了食饵一捕食者模型的x it ) 的观测数据y ，= 眵n ，y mj 的周期性，采用最 小二乘估计法时，取三角函数作为基函数，从而可把 z ( f ) = b 。( f ) ，x ：( fj ) 表示成三角函数形式，其中 墨o ) 一耋c 雎c o s ( t + ) ，f = 1 ，2 ( 1 9 ) 采用h p m 法时，通过观察i s n l 图像估计阶k ，并同时估计( 1 9 ) 中的潜频 率和初始相位妒许，这是得到高精度参数估计的关键。然后，把问题转化为最 d , - 乘问题，估计其余参数值。 最后，本文针对两种群的食饵一捕食者模型的参数给出了数值模拟，通过 对比估计值夕与真实值多，可以看出吻合效果较好。 1 3 预备知识 1 3 1 相关定义及引理 定义1 1 口8 】如果时间序列 置) = 置：t g n ) 满足 ( 1 ) 对任何f ，脚? 0 0 ： 4 第1 草驯，i 。 ( 2 ) 对任何t e n ，e x , 一l u ； ( 3 ) 对任何f ，s e n ，e ( 置一肛) ( x ，一) 1 = 以一 就称 彳， 是平稳时问序列，简称平稳序列。称实数列y 。为 x ，) 的自协方差函数。 成= 以r o ，k e z 称为平稳序列 x ， 的自相关系数。 定义1 2 2 8 l 设平稳序列 置) 有自协方差函数 n ) ， ( 1 ) 如果有卜石，万】上的单调不减右连续的函数f ( z ) ，使得 “一f = i f e a a d f ( x ) ，f ( 一兀) ；o ，k z ， ( 1 1 0 ) 就称f ( a ) 是 x ，) 或 n ) 的谱分布函数，简称涪函数。 ( 2 ) 如果有【石，巧】上的非负函数厂( 九) ，使得 “一j = ! 厂( a ) a 眠k e z ， ( 1 1 1 ) 就称f ( a ) 是 x ，) 或 圪) 的谱密度函数或功率谱密度，简称谱密度或功率谱。 定义1 3 1 2 8 l 设 乞) 是个平稳序列，如果对任何f ，s ， 即胂v ( 鸺) = 孟嚣 ( 1 1 2 ) 就称 乞) 是一个白噪声，记作删( 肛，仃2 ) 同样地，如果二维零均值实值噪声 x ( 以，川) ) = x ( 甩，棚) ：( 刀，m ) e n 2 ) 满足 。i z ? r ( ，气，z ，) ) r ( ，t ：，月懂：) = = i c 。r z ，, 【( ，n 。，a 用, 憎m 。a j ) ，。= ( ( ，n ：2 ，内, 曙m 2 2 ) ) ， ( 1 1 3 ) 就称 x ( 以，历) 为二维白噪声。 定义1 4 1 2 8 j 已知观测数据毛，x ：，h ，定义 如( 小丽1 隔n 矿枞， ( 1 1 4 ) 称，( a ) 为观测数据_ ，另：，h 的周期图。 定义1 5 1 2 8 1 设统计量否是口的估计，有如下的定义： ( 1 ) 如果e 蚕一口，就称务是口的无偏估计； ( 2 ) 如果当n 一时，子一0 ，就称否是0 的渐进无偏估计： 5 弟1 蕈引育 ( 3 ) 如粜6 依概率收敛到日，就称0 是0 的相合f 占汁； ( 4 ) 如果蚕a s 收敛到臼，就称蚕是p 的强相合估计 定义1 6 1 8 1 线性平稳序列是有白噪声的线性组合构成的平稳序列，若 皇) 是 一个零均值的线性平稳序列，则称为有色噪声。 定义1 7 0 2 9 1 设x = 伍，x 。) ，y 一化，匕) 是两个随机向量，且满足 e i x 。1 2 九，有砖) 苫砘 ) 。 推论1 2 1 3 1 若存在x 使得p b ) = o ，则p g 。) 一0 ，当a 呻。 定理1 3 1 2 8 1 设k ) 是独立同分布的删( o ，口2 ) ，线性滤波器 c j ) 满足条件 o o 平稳噪声 毒) 由皇2 荟c ，一，定义。则有如下结果 “m s u p 击i ns u p f 萎矿卜网一“p 面丽s u p i 荟矿。is 挑：吖娜， 其中 m ，? 魏叩l 是 皇 的谱密度。 定理1 4 弘卅设模型( 1 4 ) 中的平稳噪声 皇j 满足定理1 3 中的条件，则有如下 的结果： l i m s u p 丽n 3 陪- | = 0 魅， l i m s u p i n 。6 ，- 口肛o ，觚 在参数估计问题中，估计量的口j 收敛速度一般只达到 d ( ( l n l n ) 叫2 ) ， 这罩角频率 的周期图最大估计爻的口j 收敛度可以达到。( ( ，l n l n ) 叫2 ) 1 3 2 潜周期模型 实际问题中有很多的数据表现出明显的周期特性。潜周期模型( h p m ) 是一种 处理实际问题中周期性数据的常用方法【2 5 , 2 8 , 2 9 。例如，如果考虑某地区的气温变 7 第1 苹引言 化，以每个小时的甲均左i 温为一个记录数据，则观测数据j r 。，x ：，j 以每年的 3 6 5 x2 4 d , 时为一个火周期，以每天的2 4 小时为一个小周期。对于这种具有明显 周期性的数据町以用h p m 法水研列2 8 1 。 考虑简单的潜剧期模型一余弦波信号 x 。= 月c o s h 驴小舅，t e n + ( 1 1 6 ) 其中0 c o 。 甜： o 。取三角函数作为基 函数，有 z p ) = c o + 【c 2 s i n ( k wf ) + c 2 c o s ( k c of ) 】= 靠铱) = c 。驴 ( 2 6 ) 篇局 其中，c 。，k = 0 工，k 为待估参数，国= 2 x t ，c = ( c 。) ，矽= 协o ) ) 。 考虑参数的最小二乘估计，即最1 1 , 化 脚啡) 2 再【y 一荟c 舭，) 】2 ( 2 7 ) = ( y fc ) ( y - fc ) = l i y fc l l 2 ( 2 8 ) ( 2 9 ) 其中，f = ( 识( f f ) ) | 】i fx ( x + 1 ) 。 对( 2 8 ) 式关于c 求导，得到正规方程 2 f fc 一2 f y 一0 解得参数c 的最小二乘估计为： ；( f f ) 。1 f y( 2 1 0 ) 直接验证即知，当t ；h ， y ，歹= 1 ，2 ，m 】为等距观测，且m = 2 n - 1 时， ( 2 1 0 ) 式中参数c 。即( 2 4 1 ) 中4 。 下面考虑，o ) 的数值微分。 由f d a 理论，可取z ( f ) 作为离散数据的微分近似值，故对( 2 4 ) 式两边求导， 并令它近似f o ) 的导数，得： n 7 。 ) = 罗4 。g - o ns i n ( o j , x + 吼) ( 2 1 1 ) 1 2 第2 章周期一! 随机数引 | ，j 数值微分 2 2 误差估计 下面讨论以上数值微分的误差。 为方便起见，不妨假设f ( x ) 是定义住r ：周期为1 的函数， n e n ，a = 0 ；x o 0 是参数。 在周期为1 ，且属于日2 ( o ，1 ) 的周期函数中，找- - i y ，使得 m ( y ) 一! n 冀( 元一y ( x - ) ) 2 , ol l y l l l 2 r 。力 达到最小，然后将y 。作为的近似。这等价于在周期为1 ，且属于日2 ( o ，1 ) ，及 昙毫( 元一y ( _ ) ) 2 s 6 2 的光滑周期函数中，求函数弦，使l l y 1 1 2 r ( 吣) si l y 8 2 r ( 。，1 ) ，然后 1 3 第2 兰望型型堕垫塾塑塑塾堕垡坌一 一一一 9 4 y 。作为，的近似。 - 3 1 理2 1 1 3 5 3 6 1 设西是由( 2 1 2 ) 定义的泛函。对于任意满足) ( o ) 2y ( 1 ) ， _ ) ，( o ) ；y ( 1 ) ，y 。( o ) = y “( 1 ) 的函数y ，成立 m ( 弦) s 巾( y ) ， 即弘是泛函的极小点。 - 3 1 理2 2 1 3 s ，3 6 被插值函数， ) 是定义在( o ，1 ) i r - i b 七的光滑剧期函数， 7 ) = 4 ，c o s ( j x + 妒， 是( o ，1 ) 区间上的周期型插值函数，且，( 鼍) = 歹瓴 ，i = 0 ，1 ，则 l l 于”一厂”l i 二。一i i ，”l l ：伊m s l i ，? ；。山 ( 2 1 3 ) 定理2 2 设被插值函数厂 ) 是( o ，1 ) 区间上的光滑周期函数。，o ) 的一组观 测数据为厂( ) ，记y o ) 是基于厂( t ) ，i = o ，1 ，n 的刷期犁插值函数， l 厂“) 一y “) l 6 ， 6 为观测误差。夕 ) 是西( y ) 泛函极小元，则帆( o ，1 ) 有误差估计式 i l e l l 玖蚺彤吼，s 乒瓦i 砥+ y 弘 证明记误差 p g ) 一l y ( 工) 一，( x ) l ( 2 1 4 ) 利用连续型周期函数在一个周期长度内的两个端点上的函数值必定相等， 。 易知 p ( 0 ) ；e ( 1 ) ，e ( o ) = p ( 1 ) ，e ”( 0 ) = e u ( 1 ) 要得到被插函数厂o ) 的最佳逼近函数厂 ) ，厂o ) 相对于y o ) 必须满足条 件： 俐艮， m 艮， ( 2 1 5 ) 1 4 禾i j 用弓i 理2 2 1 3 6 1 有： 第2 章周期型随机数据的数 f “数分 l k ”i 臣：。舯= 0 y 。一歹”0 ：。m s l l y ” i 2 ( 0 , d4 - l i 夕”i l i ：。m s2 8 y ”l i i ：。，。， 俐乙u ，= z 已2 出一n 荟- i e 2 ( 喜) j i l = 其中岛x i ) x i + i ) ，“= 0 ，1 ，一1 ) ，办一t + 。一鼍= 万l 。 我们假设观测值厂 。) 与真实值y ) 之间存在误差6 利用重要不等式：( 口+ 6 ) 2s2 ( a 2 + 6 2 ) ，有 则 ( 2 1 6 ) ( 2 1 7 ) 耋嘶弘乏n ( y 一，) 2 墨2 耋( ( ，一地) ) 2 + ( ，( 一y ) 2 ) 叫+ 1 ) 6 2 这是因为( 夕 ) 一厂“) ) 2s6 2 和( ，“) 一y 瓴) ) 2s 62 。将之代入( 2 1 7 ) 可得 整理可得 即 ，s4 ( + 1 ) 6 2 + 薹( 已2 ( 岛) _ 6 2 n 1 ( 玉) 弦 一 s 4 ( n + 1 ) 6 2 4 ( n + 1 ) 6 2 = 4 + 1 6 2 + n s8 ( 5 24 - 2 h + 蓦( e 2 ( 岛) 一p 2 ( + 劾n 荟- i r 栅矗出 劾n 荟- i ，韶缸： 玉) 弦 8 9 ：蚴等r 向，而 i e e 缸 怫蜊：秘出) ；y o l e 。1 2 出) ： = 8 , 5 2 + 2 hi e i l z ( o a ) 叫k 。) ( 1 | p l i l ：( o a ) - - 8 e i l 。j ， ) 2 s8 62 + 2 h i i e l i 产。 e 印，1 ) 2 b + 劢er 以o ，1 ) 1 5 ( 2 1 8 ) 、- ， 誓 ，- 2 p x 向 丘 、，、， t ，l 2 p一 、， 莓 ，l 2 e ，i m v 向 2 6、，l 2 、，、， x ，l y 一 、- 、 x l ，- 、 x 向 ，k v 向 第2 章周期型随机数据的数值微分 现在讨论矧b m i i e i i ：帆，= z p 坨出e e i _ o 叫饥；0 叫 d xs oe 怍i 缸州l e i 以， 将( 2 1 6 ) 、( 2 1 8 ) 两式代入i ：式得 m 舯se 岐叩，l i e 峙吣) s 压m b ，1 ) ( 2 届+ 2 , 怛k ，) ( 2 1 9 ) 整理得 m m 一2 ( 磊。i l y 弘。 1 ) ) 咿k 1 ) n4 6 i l y 弘”)( 2 2 0 ) 配方得 ( 1 p i i r ( 。，1 ) 一2 矗。l i y ”0 r ( 。1 ) ) 2 4 6 i l y ”8 f ( 。1 ) + ( 2 厅i l y ”k ：( 。1 ) ) 2 于是得到误差估计式 i l e l l 姗，= 杪一，眠s 4 6 杪舯+ 2 肛m ，+ 扬l yi 啪， 证毕。 由误差估计的结果知， 于0 ) = 4c o s ( x + ) 作为插值函数，其误差估 万：o 计存在一个上界，不仅函数值比较接近，其导数值也比较接近。 推论1 若6 = 0 ，即观测数据无误差，则慨( o ，1 ) ，有舰7 o ) = ，o ) 。 推论2 用7 。 ) 作，0 ) 的近似值，其计算是数值稳定的。 1 6 第3 单炳种群食饵捕食者模型参数的两阶段最小二乘讨 第3 章两种群食饵捕食者模型参数的两阶段最, j , - - 乘估计 不失一股性，考虑两种群的食饵一捕食者模型 二= r ( x ，ti 卢)( 3 1 ) 其中 石( f ) = b 。p ) ，工：( f ) ) l ，f = ( 厂1 ，厶y ， = ( a 。，口。，j 。( o ) ，x ：( o ) y 一( 口。，口。，口，口6 ) ) ，f 【o ，l 】，乙尺 ，满足定理1 1 中的条件，使得( 3 1 ) 式存在惟一解，记为x ( t ，卢) 。 设乙= k d ，z 帆) 为一( f ) 在时刻，t 加，的观测数据，有 z 口= x i ( t ，) + 8 j , f = 1 ，2 ，歹t 1 ，2 ，咒 ( 3 2 ) 其中，i i d ( o ，o r 2 ) ，0 o r 0 ，处于临界状态，不能用判 断非线性方程平衡点稳定性的准则讨论( 4 1 ) 的平衡点昂的情况。用分析相轨线的 方法来讨论平衡点的稳定性。方程的相轨线为 c i = x - l o l n x - 0 妙+ 2 i n y ， 即 c = x l o y 之e 叭p 。， ( 4 3 ) c 为任意常数。 为了研究( 4 1 ) 式的确定的相轨线的图形，记 妒( x ) = 一0 。e 。 妒( y ) = y - 2 e 仉抄 若妒和妒的极大值分别记作和妒。，则不难确定石。，y o 满足 驴( 工。) = ，x 。= 1 0 ( 4 4 ) 妒( y 。) = 妒。，y 。= 2 0 ( 4 5 ) 显然仅当( 4 3 ) 式右端常数cs 妒。时，相轨线才有定义。 当c = 妒。时，z = x 。，y = y 。，n ( 4 。4 ) ，( 4 5 ) - e jp o ( 1 0 , 2 0 ) 比较可知( z 。，y 。) 正 是平衡点p o ，所以p o 是相轨线的退化点，昂是稳定平衡点。 下面给出方程( 4 1 ) ，( 4 2 ) 的两阶8 n d - 乘估计法模拟。 函数l 如 ) | 的图形如下 第4 章数值杖拟 圉4 3 ( a ) 序列z o ) 的f s 。( ) 图43 m ) 序列y q ) 的i 如( ) f 计算得到峰群的个数为4 ，从而周期个数的估计f = 4 。 分别计算出x ( t ) 和y q ) 的潜周期模型逼近，结果分别为 和 。口) = 98 9 3 1 + 善j ，c 。s ( 西f 十庐，) ，+ 第4 章数值模拟 y ( t ) = 2 0 8 2 0 7 + 荟- ，c 。s ( 击r + 妒) + ( 4 6 ) $ t 1 ( 4 7 ) - - 符系数值见表4 1 和表4 2 。 表4 1 式( 4 6 ) 中各系数的值 k 123 4 咄 3 8 6 7 77 7 4 4 81 1 6 1 9 01 5 6 0 6 0 吼 0 7 9 2 50 6 3 3 10 3 2 7 50 6 9 8 3 4 6 9 3 0 52 5 3 8 60 9 2 0 60 3 7 7 6 表4 2 式( 4 7 ) 中各系数的值 七1234 蛾 3 8 6 4 87 7 4 7 91 1 6 1 8 81 5 4 8 9 6 饥 0 6 6 7 81 3 4 2 11 8 3 4 12 2 6 2 2 4 2 9 2 1 3 51 6 4 1 4 18 4 0 0 64 0 9 0 3 由表4 1 和表4 2 ，可求得x ( f ) 和y ( t ) 的周期为 _ 2 _ x 。1 6 2 4 5 q ( 4 7 ) 与前面的估计值大致相等。误差的自协方差序列如图4 4 所示，从图中看出，自 协方差函数接近于零，表明模型通过统计检验，方程( 4 6 ) 乖- 1 1 ( 4 7 ) 可作为x ( f ) 、y ( t ) 的逼近。作出x ( f ) 、y ( t ) 的观测值和利用潜周期模型得出的估计值，见图4 5 ， 近似效果很好。 2 4 第4 章数值模拟 b l 舢1 00 2 c 删d eo f x 图4 4 潜周期模型自协方差函数统计检验酗 q 卿d a t a 镕p f 咖e d 训o f x 51,t1厂_万西，r了j 一, 铲, ! y i l i i i t 矿i ! 矿i ! i , j 0 石百南而 2 r谕1 i l m 8 胛d n 由h 镕p m 瞳c i 酣l u e o f y l wj 一，吲一一 图4 5 观测序列和潜周期模型逼近序j 7 u ( o 为原始数据为预测数据1 然后，利用线性最小二乘法，求出各参数估计值。综合z ( 0 ) 、y ( 0 ) 的估计 值，结果列于表43 。 一 一 一 。 一5 ， 。 5 5 j _口e8 第4 章数值模拟 表4 3 参数估计值及其误差 待估参数x a 1一a 2a 3一口4口5a 6 真实值 20 11 011 37 2 估计值 2 1 2 9 50 1 0 1 89 9 8 9 00 9 8 2 31 3 0 2 57 2 4 9 4 相对误差( ) 6 4 7 5 1 8 0 0 0 0 1 1 0 0 1 7 7 0 00 1 9 2 30 6 8 6 1 置信区间 ( 1 9 6 9 6 ，( 0 0 9 8 3 ，( 7 4 9 2 0 ，( 0 7 4 9 2 ， ( 口= 0 9 5 )2 3 2 3 0 )0 1 0 5 4 ) 1 2 4 6 58 ) 1 2 1 5 3 ) 检验统计量 r 2 = 0 9 4 3 4 6 ，r 2 - - 0 8 2 0 1 ， f = 2 4 8 6 3 6 6 6 ，p = 0 0 0 0 0f = 6 7 9 0 4 3 3 ，p = 0 0 0 0 0 本文计算所用计算机配置环境为：w i n d o w s x p , i n t e lp e n t i u m2 4 0 g h z ，2 5 6 m b 内存，计算平台为m a t l a b 7 2 。整个程序耗费的计算时间为6 6 9 0 3 3 5 秒，计算 速度很快，明显优于文献 3 ， 4 中非线性优化算法。 由表4 3 最后一行各统计量的值可知，模型的各参数的估计通过检验，拟合 效果很好，又置信区间( 置信度o 9 5 ) 不含0 点，且均包含真值，所以估计结果是 可靠的，结果具有很高的置信度。从表4 3 还可知，各