捕食模型参数的确定1

捕食模型参数的确定假设有一个生态系统，其中含有两种生物，即：A生物和B生物，其中A生物是捕食者，B生物是被捕食者。请建立捕食数学模型，利用有关数据，解决以下问题：1)在观测数据（DATA1）无误差的情况下，确定模型中的参数，并分析误差。2)在观测资料有误差（时间变量不含有误差）的情况下，请分别利用观测数据DATA2和DATA3，确定参数在某种意义下的最优解，并与仿真结果比较，进而改进你们的数学模型。3)假设连观测资料的时间变量也含有误差，试利用数据DATA4，建立数学模型，确定参数在某种意义下的最优解。通过对此生态系统的观测，可以得到相关的观测数据。观测数据的格式依次为：观测时刻jt、A生物数目)(jtx、B生物数目)(jty对于生态系统中的两种生物A和B，A生物为捕食者，B生物为被捕食者。在某一段时期内，A生物的数量与B生物的数量之间存在一定的关系。根据已知条件，可将（15）式改写为如下形式：12()dxxydt（1）34()dyyxdt（2）0506()()xtyt其中16kk为模型的待定参数。进行变换可得：3412()()yxdydxxy（3）3412()()dxxdyyyx即（4）积分得：10203040lnln)()(lnln)()0yyyyxxxx（可将上述表达式改写成n元齐次线性方程组的形式，如下所示：mnA0（5）上述n元齐次线性方程组有非零解的充分必要条件是系数矩阵的秩R(A)n。我们首先用DATA1中的3组数据确定，，，，4321aaaa程序clearA=zeros(3,4);A(1,1)=log(0.1374480266382216/60);A(1,2)=60-0.1374480266382216;A(1,3)=-log(11.750840650304518/10);A(1,4)=10-11.750840650304518;A(2,1)=log(7.108705996120129/60);A(2,2)=7.108705996120129-60;A(2,3)=-log(3.4133367257849176/10);A(2,4)=10-3.4133367257849176;A(3,1)=log(0.4251595082899424/60);A(3,2)=0.4251595082899424-60;A(3,3)=-log(20.80921881438798/10);A(3,4)=10-20.80921881438798;r=rank(A);%rank(A)=rn时,该方程有无穷多个解,求它的一个基本解y=null(A,r)表1)41(a'kk的值'1a'2a'3a'4a-0.0478-0.0042-0.99250.11253314140000222222(lnln)lnlnyyyxxyxx（28）如设：31400000222(lnln)yyxx，112，422，332，1x=lny，2x=x，3lnxx，则（28）式可以写为如下形式；0112233yxxx（29）对于（29）式中因变量y是自变量123xxxx的线性函数。可以建立起因变量y的多元线性回归模型，①利用数据文件data2.txt,首先绘制出x(t)、y(t)与时间t的关系图象，和y(t)对x(t)的散点图，如图3所示（程序见gg1）。利用MATLAB工具箱中的regress(y,x)命令，计算回归系数β的最小二乘估计ˆβ及其置信区间，计算结果如表1所示，计算程序名为gg5。由它得到的模型为：ˆy-139.4322355563643+19.884217746922571x-9.9882903314159892x+99.804970720159603x（51）结果分析：表1显示2R=0.9971088966600568是指因变量y的99.71%可由模型（51）确定，F值远远超过F检验的临界值，p远小于，因而模型（30）从整体上看是可用的。表1的回归系数给出了模型（30）中0，1，2，3的估计值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。图3y(t)对x(t)的散点图和x(t)、y(t)与时间t的关系图得出β的估计值即可确定i（i=1,2,3,4）之间的相互关系：112，332，422，再利用残差ˆˆ(1,,150)iiiieyyi达到最小值进行最优化计算来确定2的取值，从而确定其它参数值，模型（30）的残差与2x的散点图如图4所示。表1data2.txt中β的计算结果回归参数参数估计值参数置信区间0-139.4322355563643[-142.1292050800128-136.7352660327159]119.88421774692257[19.7090526267255320.05938286711962]2-9.988290331415989[-10.20688059603675-9.769700066795226]399.80497072015960[97.65938535619841101.9505560841208]2R=0.9971088966600568F=16784.58380705823p0.0001从图中可以看出随着2x的增加，残差ie也逐渐增大。通过优化计算得出在捕食者A和被捕食者B共同存在，相互竞争达到动态平衡的状况下参数2的最优值为2=0.1，从而可以确定参数1，3，4的最优值。1=-1.988421774692257，2=0.1，3=9.980497072015961，4=-0.998829033141599，5=12.962285633035274，6=72.12301583344777利用上述数据对模型（30）的拟合效果如图5所示，相关程序为lorenzep.m和gg8.m。从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30）式所示的模型时，假设回归变量123xxxx对因变量y的影响是相互独立的。然而根据直觉和经验可以猜想，1x、2x、3x之间的交互作用会对y产生影响，不妨简单地利用1x、2x、3x的乘积代表它们的交互作用，于是将模型（30）增加一项得到：011223341232(0,)yxxxxxxN（52）图4残差ie与2x的散点图图5模型（30）的数据拟合②利用数据文件DATA3.txt,首先绘制出x(t)、y(t)与时间t的关系图象，和y(t)对x(t)的散点图，如图6所示（程序见gg4）。利用MATLAB工具箱中的regress(y,x)命令，计算回归系数β的最小二乘估计ˆβ及其置信区间，计算结果如表2所示，计算程序名为gg6。图6y(t)对x(t)的散点图和x(t)、y(t)与时间t的关系图由它得到的模型为：y-120.3100735257612+17.9086484171411x-8.8031135947385142x+87.89290670810463x（53）结果分析：表2显示2R=0.9228752184692365是指因变量y的92.29%可由模型（53）确定，F值远远超过F检验的临界值，p远小于，因而模型（30）从整体上看是可用的。表2的回归系数给出了模型（30）中0，1，2，3的估计表2DATA3.txt中β的计算结果回归参数参数估计值参数置信区间0-120.3100735257612[-124.4170908611838-116.2030561903386]117.9086484171410[17.6460181852192018.17127864906279]2-8.803113594738514[-9.138961160791356-8.467266028685673]387.89290670810462[84.6065452045026591.17926821170659]2R=0.9228752184692365F=5967.045867867118p0.0001值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。得出β的估计值即可确定i（i=1,2,3,4）之间的相互关系：112，332，422，再利用残差ˆˆ(1,,1500)iiiieyyi达到最小值进行最优化计算来确定2的取值，从而确定其它参数值，模型（30）的残差与2x的散点图如图7所示。从图中可以看出随着2x的增加，残差ie也逐渐增大。通过优化计算得出在捕食者A和被捕食者B共同存在，相互竞争达到动态平衡的状况下参数2的最优值为2=0.115，从而可以确定参数1，3，4的最优值。1=-2.059494567971215，2=0.115，3=10.10768427143203，4=-1.012358063394929，5=12.82305785349266，6=73.40728069086651利用上述数据对模型（30）的拟合效果如图8所示，相关程序为lorenzep1.m和gg10.m。图7残差ie与2x的散点图从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30）式所示的模型时，假设回归变量123xxxx对因变量y的影响是相互独立的。然而根据直觉和经验可以猜想，1x、2x、3x之间的交互作用会对y产生影响，不妨简单地利用1x、2x、3x的乘积代表它们的交互作用，于是将模型（30）增加三项得到：011223341251362341232(0,)yxxxxxxxxxxxxN（54）