预测多元时间序列
这个例子展示了如何在猎物拥挤的情况下，对捕食者和猎物种群测量的数据进行多元时间序列预测。捕食者-被捕食者种群变化动态采用线性和非线性时间序列模型。比较了这些模型的预测性能。
数据描述
该数据是一个双变量时间序列，由1-捕食者1-猎物种群(以千为单位)组成，每年收集10次，持续20年。有关数据的详细信息，请参见三种生态种群系统:MATLAB和C mexo - file时间序列建模．
加载时间序列数据。
负载PredPreyCrowdingDataZ = iddata(y，[]，0.1，“TimeUnit”，“年”，“Tstart”, 0);
z是一个iddata对象，其中包含两个输出信号，日元而且y2，分别表示捕食者种群和被捕食者种群。的OutputData的属性z以201 × 2矩阵的形式包含总体数据，例如z.OutputData (: 1)捕食者的数量和z.OutputData (: 2)是猎物的数量。
绘制数据图。
情节(z)标题(“捕食者-猎物种群数据”) ylabel (“人口(千)”）
数据显示，由于拥挤，食肉动物的数量有所下降。
使用前半部分作为识别时间序列模型的估计数据。
Ze = z(1:120);
利用剩余数据搜索模型顺序，并验证预测结果。
Zv = z(121:结束);
线性模型的估计
将时间序列建模为线性自回归过程。可以使用以下命令以多项式形式或状态空间形式创建线性模型基于“增大化现实”技术(仅适用于标量时间序列)，arx，armax，n4sid而且党卫军．由于线性模型不捕获数据偏移量(非零条件均值)，首先对数据进行趋势化。
[zed, Tze] =趋势(ze, 0);[zvd, Tzv] =趋势(zv, 0);
识别需要型号订单的说明。对于多项式模型，可以使用arxstruc命令。自arxstruc只适用于单输出模型，对每个输出分别执行模型顺序搜索。
Na_list = (1:10)';V1 = arxstruc(zed(:，1，:)，zvd(:，1，:)，na_list);na1 = selstruc(V1,0);V2 = arxstruc(zed(:，2，:)，zvd(:，2，:)，na_list);na2 = selstruc(V2,0);
的arxstrucCommand分别表示阶7和阶8的自回归模型。
使用这些模型顺序来估计任意选择交叉项的多方差ARMA模型。
Na = [na1 na1-1;na2-1钠);Nc = [na1;钠);sysARMA = armax(zed，[na nc])
sysARMA =离散ARMA模型:模型输出“日元”:一个(z) y_1 (t) = - ai (z) y_i (t) + C (z) e_1 (t) (z) = 1 - 0.885 z ^ 1 - 0.1493 z ^ 2 + 0.8089 z ^ 3 - 0.2661 z ^ 4 - 0.9487 z ^ 5 + 0.8719 z ^ 6 - 0.2896 z ^ 7 A₂(z) = 0.3433 z ^ 1 - 0.2802 z ^ 2 - 0.04949 z ^ 3 + 0.1018 z ^ 4 - 0.02683 z ^ 5 - 0.2416 z ^ 6 C (z) = 1 - 0.4534 z ^ 1 - 0.4127 z ^ 2 + 0.7874 z ^ 3 + 0.298 z ^ 4 - 0.8684 z z ^ 6 ^ 5 + 0.6106 + 0.3616 z ^ 7模型输出“y2”:(z) y_2 (t) = - ai (z) y_i (t) + C (z) e_2 (t) (z) = 1 - 0.5826 z ^ 1 - 0.4688 z ^ 2 - 0.5949 z ^ 3 - 0.0547 z ^ 4 + 0.5062 z ^ 5 + 0.4024 z ^ 6 - 0.01544 z ^ 7 - 0.1766 z ^ 8 A_1 (z) = 0.2386 z ^ 1 + 0.1564 z ^ 2 - 0.2249 z ^ 3 - 0.2638 z ^ 4 - 0.1019 z ^ 5 - 0.07821 z ^ 6 + 0.2982 z ^ 7 C (z) = 1 - 0.1717 z ^ 1 - 0.09877 z ^ 2 - 0.5289 z ^ 3 - 0.24 z ^ 4 + 0.06555 z ^ 5 + 0.2217 z ^ 6 - 0.05765 z ^ 7 - 0.1824 z ^ 8样时间:0.1年参数化:多项式订单:na=[7 6;7 8] nc=[7;8]自由系数数:43使用"polydata"， "getpvec"， "getcov"表示参数及其不确定度。状态:在时域数据“zed”上使用ARMAX估计。拟合估计数据:[89.85;90.97]%(预测焦点)FPE: 3.814e-05, MSE: 0.007533
计算提前10步(1年)的预测输出，以在估计数据的时间跨度内验证模型。由于用于估计的数据是去趋势的，因此需要为有意义的预测指定那些偏移量。
predOpt = predictOptions(“OutputOffset”Tze.OutputOffset ');yhat1 = predict(sysARMA,ze,10, predOpt);
的预测Command预测测量数据在时间跨度内的响应，是验证估计模型质量的工具。时间响应t用t = 0时的测量值计算，…, t-10。
绘制预测响应和实测数据。
情节(泽,yhat1)标题(“与测量数据相比，10步预测反应”）
注意，预测响应的生成和与测量数据的绘图，可以自动使用比较命令。
compareOpt = compareOptions(“OutputOffset”Tze.OutputOffset ');比较(泽sysARMA 10 compareOpt)
使用比较还以百分比形式显示了拟合优度的归一化均方根(NRMSE)度量。
在验证数据后，对模型的输出进行预测sysARMA100步(10年)以外的估计数据，并计算输出标准差。
forecastOpt = forecastOptions(“OutputOffset”Tze.OutputOffset ');[yf1,x01,sysf1,ysd1] = forecast(sysARMA, ze, 100, forecastOpt);
yf1预测的响应，作为一个返回iddata对象。yf1。OutputData包含预测值。
sysf1系统是否类似于sysARMA而是状态空间形式。模拟sysf1使用sim卡在初始条件下，x01，重现预测的响应，yf1．
ysd1是标准差矩阵。它测量了由于数据中添加性扰动的影响而导致的预测不确定性sysARMA。NoiseVariance)，参数不确定性(如getcov (sysARMA))和与将过去数据映射到预测所需的初始条件的过程相关的不确定性(见data2state)．
绘制模型的测量、预测和预测输出sysARMA．
t = yf1.SamplingInstants;te = ze.SamplingInstants;t0 = z.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat1.y (: 1),.．.t, yf1.y (: 1),“米”，.．.t, yf1.y (: 1) + ysd1 (: 1),“k——”，.．.t, yf1.y (: 1) -ysd1 (: 1),“k——”)包含(的时间(年)）;ylabel (“食肉动物数量，以千计”）;次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat1.y (:, 2),.．.t, yf1.y (:, 2),“米”，.．.t, yf1.y (:, 2) + ysd1 (:, 2),“k——”，.．.t, yf1.y (:, 2) -ysd1 (:, 2),“k——”）使数字变大。FIG = gcf;p = fig.位置;fig.Position = [p (1), (2) - p (4) * 0.2, p (3) * 1.4, p (4) * 1.2);包含(的时间(年)）;ylabel (“猎物数量，以千计”）;传奇({“测量”，“预测”，“预测”，“预测不确定性(1 sd)”}，.．.“位置”，“最佳”）
该图显示，使用线性ARMA模型(增加了对偏移量的处理)进行预测在一定程度上是有效的，结果显示，与12-20年时间跨度内的实际人口相比，具有很高的不确定性。这表明种群变化动态可能是非线性的。
估计非线性黑盒模型
在本节中，我们通过在模型中添加非线性元素来扩展前一节的黑盒识别方法。这是使用非线性ARX建模框架进行的。非线性ARX模型在概念上类似于线性模型(例如sysARMA上图)，他们使用两步过程计算预测输出:
将测量信号映射到有限维回归空间。
使用静态函数将回归函数映射到预测输出，静态函数可以是线性的或非线性的。
非线性可以在两个地方添加——在回归集(步骤1)和/或在静态映射函数(步骤2)中。在本节中，我们将探索这两个选项。首先，我们创建了非线性ARX模型的线性回归形式，其中所有非线性都仅限于回归定义;将回归函数转换为预测输出的静态映射函数是线性的(准确地说，是仿射的)。然后，我们创建了一个模型，使用高斯过程函数作为静态映射函数，但保持其回归函数线性。
多项式AR模型
创建一个2输出无输入的非线性ARX模型。
L =线性回归(ze.OutputName，{1,1});sysNLARX = idnlarx(ze。OutputName,泽。InputName, L， []);sysNLARX。Ts = 0.1;sysNLARX。TimeUnit =“年”；
sysNLARX是不使用非线性函数的一阶非线性ARX模型;它使用一阶回归量的加权和来预测模型响应。
getreg (sysNLARX)
ans =2 x1细胞{“y1 (t - 1)”}{}“y2 (t - 1)”
为了引入非线性函数，在模型中加入多项式回归函数。
创建阶为2的回归量，并包括交叉项(上面列出的标准回归量的乘积)。下载188bet金宝搏将这些回归量添加到模型中。
P =多项式回归(zee . outputname，{1,1}，2,false,true);sysNLARX.Regressors(end+1) = P;getreg (sysNLARX)
ans =5 x1细胞{y1 (t - 1)的}{y2 (t - 1)的}{y1 (t - 1) ^ 2的}{y2 (t - 1) ^ 2的}{“y1 (t - 1) * y2 (t - 1)}
使用估计数据估计模型的系数(回归权重和偏移量)，泽．
sysNLARX = nlarx(ze,sysNLARX)
sysNLARX = 2个输出的非线性多元时间序列模型输出:y1, y2回归:1。变量y1 y2中的线性回归函数。2阶回归变量y1, y2输出函数:输出1:线性带偏移输出2:线性带偏移采样时间:0.1年状态:终止条件:近(局部)最小值，(范数(g) < tol)..迭代次数:0，函数求值次数:1使用NLARX对时域数据“ze”进行估计。拟合估计数据:[88.34;88.91]%(预测焦点)FPE: 3.265e-05, MSE: 0.01048更多信息在模型的“报告”属性。
计算提前10步的预测输出来验证模型。
yhat2 = predict(sysNLARX,ze,10);
预测模型的输出比估计数据高出100步。
yf2 = forecast(sysNLARX,ze,100);
对于非线性ARX模型，未计算预测响应的标准差。这些数据是不可用的，因为在这些模型的估计过程中没有计算参数协方差信息。
绘制测量、预测和预测输出。
t = yf2.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat2.y (: 1),.．.t, yf2.y (: 1),“米”)包含(的时间(年)）;ylabel (“捕食者数量(数千)”）;标题(《ARX模型的多项式回归预测结果》)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat2.y (:, 2),.．.t, yf2.y (:, 2),“米”)传说(“测量”，“预测”，“预测”)包含(的时间(年)）;ylabel (“猎物数量(数千)”）;
图显示，使用非线性ARX模型的预测结果比使用线性模型的预测结果更好。非线性黑盒建模不需要关于控制数据的方程的先验知识。
请注意，为了减少回归量的数量，您可以使用主成分分析选择(转换后的)变量的最佳子集(参见主成分分析)或特征选择(见sequentialfs)的统计和机器学习工具箱™。
估计一个高斯过程模型
模型模型sysNLARX使用线性映射函数;非线性只包含在回归量中。一个更灵活的方法是选择一个非平凡的映射函数，比如高斯过程函数(GP)。创建一个非线性ARX模型，只使用线性回归量，但添加了一个高斯过程函数来将回归量映射到输出。
复制sysNLARX，这样我们就不必再次创建模型结构sysGP = sysNLARX;移除多项式回归集sysGP.Regressors(2) = [];将线性映射函数替换为GP函数。默认情况下，GP使用平方指数核。GP = idGaussianProcess;sysGP。OutputFcn = GP;disp (sysGP)
带属性的2x0 idnlarx数组:Regressors: [1x1 linearRegressor] OutputFcn: [2x1 idGaussianProcess] RegressorUsage: [2x4 table]归一化:[1x1 struct] TimeVariable: 't' NoiseVariance: [2x2 double] InputName: {0x1 cell} InputUnit: {0x1 cell} InputGroup: [1x1 struct] OutputName: {2x1 cell} OutputUnit: {2x1 cell} OutputGroup: [1x1 struct]备注:[0x1 string] UserData:[]名称:" Ts: 0.1000 TimeUnit: 'years'报告:[1x1 idresults.nlarxhw]
您现在有了一个基于线性回归器和GP映射函数的模板模型sysGP。它的参数是它的平方指数核的两个系数。使用nlarx来估计这些参数。
sysGP = nlarx(ze, sysGP)
sysGP = 2个输出的非线性多元时间序列模型输出:y1, y2回归函数:变量y1, y2中的线性回归函数输出函数:输出1:使用SquaredExponential核的高斯过程函数输出2:使用SquaredExponential核的高斯过程函数采样时间:0.1年状态:在时域数据“ze”上使用NLARX估计。拟合估计数据:[88.07;88.84]%(预测焦点)FPE: 3.369e-05, MSE: 0.01083更多信息在模型的“报告”属性。
getpvec (sysGP)
ans =10×1-0.6671 0.0196 -0.0000 -0.0571 -3.6447 -0.0221 -0.5752 0.0000 0.1725 -3.1068
与前面一样，计算预测和预测的响应并绘制结果。
计算提前10步的预测输出来验证模型。
yhat3 = predict(sysGP,ze,10);
预测模型的输出比估计数据高出100步。
yf3 = forecast(sysGP,ze,100);
绘制测量、预测和预测输出。
t = yf2 . samplinginstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat3.y (: 1),.．.t, yf3.y (: 1),“米”)包含(的时间(年)）;ylabel (“捕食者数量(数千)”）;标题(《高斯过程模型的预测结果》)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat3.y (:, 2),.．.t, yf3.y (:, 2),“米”)传说(“测量”，“预测”，“预测”)包含(的时间(年)）;ylabel (“猎物数量(数千)”）;
结果表明，与线性模型相比，该模型的精度有所提高。使用基于GP的模型的一个优点是，与其他形式(如小波网络(idWaveletNetwork)或Sigmoid网络(idSigmoidNetwork)．这使得它们易于训练，参数估计的不确定性较低。
非线性灰盒模型的估计
如果你有系统动力学的先验知识，你可以用非线性灰盒模型拟合估计数据。一般来说，了解动态的性质有助于提高模型质量，从而提高预测精度。对于捕食者-猎物动力学，捕食者的变化(日元)和猎物(y2)人口可表示为:
有关方程的更多信息，请参见三种生态种群系统:MATLAB和C mexo - file时间序列建模．
在此基础上建立非线性灰盒模型。
指定一个描述捕食者-猎物系统模型结构的文件。该文件将状态导数和模型输出指定为时间、状态、输入和模型参数的函数。选择两个输出(捕食者和猎物种群)作为状态，推导出动态的非线性状态空间描述。
文件名=“predprey2_m”；
指定模型顺序(输出、输入和状态的数量)
Order = [2 0 2];
指定参数的初始值
，
，
，
,
，并表示所有参数都要估计。请注意，在估计黑盒模型时，不存在指定参数初始猜测的需求sysARMA而且sysNLARX．
参数= struct(“名字”, {“生存因素，捕食者”“死亡因素，掠食者”.．.“生存因素，猎物”“死亡因素，猎物”.．.“拥挤因素，猎物”}，.．.“单位”, {1 /年的1 /年的1 /年的1 /年的1 /年的}，.．.“价值”，{-1.1 0.9 1.1 0.9 0.2}，.．.“最低”，{-Inf -Inf -Inf -Inf -Inf -Inf}，.．.“最大”，{Inf Inf Inf Inf Inf Inf}，.．.“固定”，{假假假假假});
类似地，指定模型的初始状态，并指出两个初始状态都要估计。
InitialStates = struct(“名字”, {捕食者种群的“猎物人口”}，.．.“单位”, {的大小(千)的大小(千)}，.．.“价值”{1.8 - 1.8},.．.“最低”， {0 0}，.．.“最大”{正正无穷},.．.“固定”,{假假});
将模型指定为连续时间系统。
Ts = 0;
创建具有指定结构、参数和状态的非线性灰盒模型。
sysGrey = idnlgrey(文件名，顺序，参数，初始化状态，Ts，“TimeUnit”，“年”）;
估计模型参数。
sysGrey = nlgreyest(ze,sysGrey);礼物(sysGrey)
sysGreyGrey = ' preprey2_m '定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, x(t)， p1，…， p5) y(t) = H(t, x(t)， p1，…，p5）+e(t) with 0 input(s), 2 state(s), 2 output(s), and 5 free parameter(s) (out of 5). States: Initial value x(1) Predator population(t) [Size (thou..] xinit@exp1 2.01325 (estimated) in [0, Inf] x(2) Prey population(t) [Size (thou..] xinit@exp1 1.99687 (estimated) in [0, Inf] Outputs: y(1) y1(t) y(2) y2(t) Parameters: Value Standard Deviation p1 Survival factor, predators [1/year] -0.995895 0.0125269 (estimated) in [-Inf, Inf] p2 Death factor, predators [1/year] 1.00441 0.0129368 (estimated) in [-Inf, Inf] p3 Survival factor, preys [1/year] 1.01234 0.0135413 (estimated) in [-Inf, Inf] p4 Death factor, preys [1/year] 1.01909 0.0121026 (estimated) in [-Inf, Inf] p5 Crowding factor, preys [1/year] 0.103244 0.0039285 (estimated) in [-Inf, Inf] Status: Termination condition: Change in cost was less than the specified tolerance.. Number of iterations: 6, Number of function evaluations: 7 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "ze". Fit to estimation data: [91.21;92.07]% FPE: 8.613e-06, MSE: 0.005713 More information in model's "Report" property. Model Properties
计算提前10步的预测输出来验证模型。
yhat4 = predict(sysGrey,ze,10);
预测模型的输出超出估计数据100步，并计算输出标准差。
[yf4,x04,sysf4,ysd4] = forecast(sysGrey,ze,100);
绘制测量、预测和预测输出。
t = yf4.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat4.y (: 1),.．.t, yf4.y (: 1),“米”，.．.t, yf4.y (: 1) + ysd4 (: 1),“k——”，.．.t, yf4.y (: 1) -ysd4 (: 1),“k——”)包含(的时间(年)）;ylabel (“捕食者数量(数千)”）;标题(《灰盒模型的预测结果》)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat4.y (:, 2),.．.t, yf4.y (:, 2),“米”，.．.t, yf4.y (:, 2) + ysd4 (:, 2),“k——”，.．.t, yf4.y (:, 2) -ysd4 (:, 2),“k——”)传说(“测量”，“预测”，“预测”，“预测不确定性(1 sd)”)包含(的时间(年)）;ylabel (“猎物数量(数千)”）;
结果表明，采用非线性灰盒模型进行预测，预测结果良好，预测输出不确定性低。
比较预测表现
比较从已确定的模型中获得的预测响应，sysARMA，sysNLARX,sysGrey．前两个是离散时间模型和sysGrey是一个连续时间模型。
CLF plot(z,yf1,yf2,yf3,yf4)图例({“测量”，“ARMA”，“多项式AR”，“全科医生”，“方框”})标题(“预测反应”) xlim([7 23])
非线性ARX模型的预测效果优于线性模型。在非线性灰盒模型中加入动力学知识，进一步提高了模型的可靠性，从而提高了预测精度。
请注意，灰盒建模中使用的方程与非线性ARX模型中使用的多项式回归量密切相关。如果你用一阶差分近似控制方程中的导数，你将得到类似于:
在那里, 是否有些参数与原参数有关 以及用于差分的样本时间。这些方程表明，在构建非线性ARX模型时，第一个输出(y1(t-1)和*y1(t-1)*y2(t-1))只需要2个回归量，第二个输出只需要3个回归量。即使没有这样的先验知识，线性回归模型结构使用多项式回归仍然是一个流行的选择在实践中。
用非线性灰盒模型预测200年以上的数值。
[yf5，~，~，ysd5] = forecast(sysGrey, ze, 2000);
绘制数据的后一部分(显示1sd的不确定度)
t = yf5.SamplingInstants;N = 700:2000;次要情节(1、2、1);情节(t (N) yf5.y (N, 1),“米”，.．.t (N), yf5.y (N - 1) + ysd5 (N, 1),“k——”，.．.t (N), yf5.y (N, 1) -ysd5 (N, 1),“k——”)包含(的时间(年)）;ylabel (“捕食者数量(数千)”）;标题(“200年后的人口预测”) ax = gca;斧子。YLim = [0.8 1];斧子。XLim = [82 212];次要情节(1、2、2);情节(t (N) yf5.y (N, 2),“米”，.．.t (N), yf5.y (N, 2) + ysd5 (N, 2),“k——”，.．.t (N), yf5.y (N, 2) -ysd5 (N, 2),“k——”)传说(人口预测的，“预测不确定性(1 sd)”)包含(的时间(年)）;ylabel (“猎物数量(数千)”）;Ax = gca;斧子。YLim = [0.9 1.1];斧子。XLim = [82 212];
图中显示，预计掠食性种群数量将达到890左右的稳定状态，而猎物种群数量预计将达到990。
相关的话题
您也可以从以下列表中选择一个网站: