Этот пример показывает, как выполнять многомерное прогнозирование временных рядов данных, измеренных из популяций хищников и добычи в сценарии скопления добычи. Динамика изменения популяции хищника-добычи моделируется с использованием линейных и нелинейных моделей временных рядов. Сравниваются прогнозные характеристики этих моделей.
Данные представляют собой двумерный временной ряд, состоящий из 1-хищных 1-хищных популяций (в тысячах), собираемых 10 раз в год в течение 20 лет. Дополнительные сведения об этих данных см. в разделе Три системы экологического населения: MATLAB и C MEX-File Modeling of Time-Series.
Загрузите данные временных рядов.
load PredPreyCrowdingData z = iddata(y,[],0.1,'TimeUnit','years','Tstart',0);
z является iddata объект, содержащий два выходных сигнала, y1 и y2, которые относятся к популяциям хищников и добычи соответственно. OutputData имущество z содержит данные о населении в виде матрицы 201 на 2, так что z.OutputData(:,1) является популяцией хищника и z.OutputData(:,2) является популяцией добычи.
Постройте график данных.
plot(z) title('Predator-Prey Population Data') ylabel('Population (thousands)')
Данные свидетельствуют о сокращении популяции хищников из-за скопления людей.
Использовать первую половину в качестве оценочных данных для идентификации моделей временных рядов.
ze = z(1:120);
Оставшиеся данные используются для поиска заказов моделей и проверки результатов прогнозирования.
zv = z(121:end);
Моделирование временного ряда как линейного авторегрессионного процесса. Линейные модели можно создавать в полиномиальной форме или в форме «состояние-пространство» с помощью таких команд, как ar (только для скалярных временных рядов), arx, armax, n4sid и ssest. Так как линейные модели не фиксируют смещения данных (ненулевое условное среднее), сначала уменьшите данные.
[zed, Tze] = detrend(ze, 0); [zvd, Tzv] = detrend(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);
arxstruc команда предлагает авторегрессивные модели порядков 7 и 8 соответственно.
Используйте эти порядки моделей для оценки модели ARMA с несколькими дисперсиями, где перекрестные термины были выбраны произвольно.
na = [na1 na1-1; na2-1 na2]; nc = [na1; na2]; sysARMA = armax(zed,[na nc])
sysARMA =
Discrete-time ARMA model:
Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + C(z)e_1(t)
A(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_2(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^-5 + 0.6106 z^-6 + 0.3616 z^-7
Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + C(z)e_2(t)
A(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
Sample time: 0.1 years
Parameterization:
Polynomial orders: na=[7 6;7 8] nc=[7;8]
Number of free coefficients: 43
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARMAX on time domain data "zed".
Fit to estimation data: [89.85;90.97]% (prediction focus)
FPE: 3.814e-05, MSE: 0.007533
Вычислите прогнозируемый результат на 10 шагов вперед (на 1 год) для проверки достоверности модели в течение периода времени данных оценки. Поскольку данные были отклонены для оценки, необходимо указать эти смещения для значимых прогнозов.
predOpt = predictOptions('OutputOffset',Tze.OutputOffset');
yhat1 = predict(sysARMA,ze,10, predOpt); predict команда предсказывает отклик на протяжении периода времени измеряемых данных и является инструментом для проверки качества оцененной модели. Ответ в момент времени t вычисляется с использованием измеренных значений в моменты времени t = 0,..., t-10.
Постройте график прогнозируемой реакции и измеренных данных.
plot(ze,yhat1)
title('10-step predicted response compared to measured data')
Следует отметить, что генерация прогнозируемого ответа и его построение с использованием измеренных данных могут быть автоматизированы с помощью compare команда.
compareOpt = compareOptions('OutputOffset',Tze.OutputOffset');
compare(ze,sysARMA,10,compareOpt)
График, созданный с помощью compare также показывает нормализованное среднеквадратичное (NRMSE) измерение доброты посадки в процентах.
После проверки данных спрогнозируйте выходные данные модели. sysARMA 100 шагов (10 лет) сверх оценочных данных и вычисление выходных стандартных отклонений.
forecastOpt = forecastOptions('OutputOffset',Tze.OutputOffset');
[yf1,x01,sysf1,ysd1] = forecast(sysARMA, ze, 100, forecastOpt);yf1 является прогнозируемым ответом, возвращенным как iddata объект. yf1.OutputData содержит прогнозируемые значения.
sysf1 является системой, аналогичной sysARMA но находится в форме state-space. Моделирование sysf1 с использованием sim команда, с начальными условиями, x01воспроизводит прогнозируемый ответ, yf1.
ysd1 - матрица стандартных отклонений. Он измеряет неопределенность прогнозирования из-за влияния аддитивных возмущений в данных (измеренных sysARMA.NoiseVariance), неопределенность параметров (как сообщается getcov(sysARMA)) и неопределенности, связанные с процессом сопоставления прошлых данных с исходными условиями, необходимыми для прогнозирования (см. data2state).
Постройте график измеренных, прогнозируемых и прогнозируемых выходных данных для модели sysARMA.
t = yf1.SamplingInstants; te = ze.SamplingInstants; t0 = z.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat1.y(:,1),... t,yf1.y(:,1),'m',... t,yf1.y(:,1)+ysd1(:,1),'k--', ... t,yf1.y(:,1)-ysd1(:,1), 'k--') xlabel('Time (year)'); ylabel('Predator population, in thousands'); subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat1.y(:,2),... t,yf1.y(:,2),'m',... t,yf1.y(:,2)+ysd1(:,2),'k--', ... t,yf1.y(:,2)-ysd1(:,2),'k--') % Make the figure larger. fig = gcf; p = fig.Position; fig.Position = [p(1),p(2)-p(4)*0.2,p(3)*1.4,p(4)*1.2]; xlabel('Time (year)'); ylabel('Prey population, in thousands'); legend({'Measured','Predicted','Forecasted','Forecast Uncertainty (1 sd)'},... 'Location','best')
Графики показывают, что прогнозирование с использованием линейной модели ARMA (с добавленной обработкой смещений) несколько сработало, и результаты показали высокую неопределенность по сравнению с фактическими популяциями в течение 12-20 лет. Это указывает на то, что динамика изменения населения может быть нелинейной.
Поместите нелинейную модель черного ящика в оценочные данные. Предварительные знания о уравнениях, регулирующих данные оценки, не требуются. Будет оценена линейно-в-регрессорная форма нелинейной модели ARX.
Создайте нелинейную модель ARX с 2 выходами и без входов.
sysNLARX = idnlarx([1 1;1 1],[],'Ts',0.1,'TimeUnit','years');
sysNLARX - нелинейная модель ARX первого порядка, которая не использует нелинейную функцию; он предсказывает отклик модели, используя взвешенную сумму его регрессоров первого порядка.
getreg(sysNLARX)
ans = 2x1 cell
{'y1(t-1)'}
{'y2(t-1)'}
Чтобы ввести нелинейную функцию, добавьте в модель полиномиальные регрессоры.
Создайте регрессоры до мощности 2 и включите перекрестные термины (продукты стандартных регрессоров, перечисленных выше). Добавьте эти регрессоры в модель в качестве пользовательских регрессоров.
R = polyreg(sysNLARX,'MaxPower',2,'CrossTerm','on'); sysNLARX.CustomRegressors = R; getreg(sysNLARX)
ans = 5x1 cell
{'y1(t-1)' }
{'y2(t-1)' }
{'y1(t-1).^2' }
{'y1(t-1).*y2(t-1)'}
{'y2(t-1).^2' }
Оцените коэффициенты (веса регрессора и смещение) модели, используя оценочные данные, ze.
sysNLARX = nlarx(ze,sysNLARX)
sysNLARX = Nonlinear multivariate time series model with 2 outputs Outputs: y1, y2 Regressors: 1. Linear regressors in variables y1, y2 2. Custom regressor: y1(t-1).^2 3. Custom regressor: y1(t-1).*y2(t-1) 4. Custom regressor: y2(t-1).^2 List of all regressors All model outputs are linear in their regressors. Sample time: 0.1 years Status: Estimated using NLARX on time domain data "ze". Fit to estimation data: [88.34;88.91]% (prediction focus) FPE: 3.265e-05, MSE: 0.01048
Вычислите прогнозируемый результат на 10 шагов вперед для проверки модели.
yhat2 = predict(sysNLARX,ze,10);
Спрогнозировать выходные данные модели на 100 шагов сверх оценочных данных.
yf2 = forecast(sysNLARX,ze,100);
Стандартные отклонения прогнозируемой реакции не вычисляются для нелинейных моделей ARX. Эти данные недоступны, поскольку информация о ковариации параметра не вычисляется во время оценки этих моделей.
Постройте график измеренных, прогнозируемых и прогнозируемых результатов.
t = yf2.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat2.y(:,1),... t,yf2.y(:,1),'m') xlabel('Time (year)'); ylabel('Predator population (thousands)'); subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat2.y(:,2),... t,yf2.y(:,2),'m') legend('Measured','Predicted','Forecasted') xlabel('Time (year)'); ylabel('Prey population (thousands)');
Графики показывают, что прогнозирование с использованием нелинейной модели ARX дало лучшие результаты прогнозирования, чем использование линейной модели. Нелинейное моделирование черного ящика не требовало предварительных знаний о уравнениях, регулирующих данные.
Обратите внимание, что для уменьшения числа регрессоров можно выбрать оптимальное подмножество (преобразованных) переменных с помощью анализа основных компонентов (см. pcaили выбор элемента (см. sequentialfs) в Toolbox™ статистики и машинного обучения.
При наличии знаний о динамике системы можно подогнать оценочные данные с помощью нелинейной модели «серого ящика».
Знание характера динамики может помочь улучшить качество модели и, таким образом, точность прогнозирования. Для динамики хищник-добыча изменения в хищнике (y1) и добыча (y2) население может быть представлено как:
y2 (t) * y1 (t)
(t) -p5 * y2 (t) 2
Дополнительные сведения о уравнениях см. в разделе Три системы экологического населения: MATLAB и C MEX-File Modeling of Time-Series.
Создайте нелинейную серую модель на основе этих уравнений.
Укажите файл, описывающий структуру модели для системы хищник-добыча. Файл определяет производные состояния и выходные данные модели как функцию времени, состояний, входных данных и параметров модели. Два выхода (популяции хищников и добычи) выбираются в качестве состояний для получения нелинейного описания состояния и пространства динамики.
FileName = 'predprey2_m';Укажите порядки моделей (количество выходов, входов и состояний)
Order = [2 0 2];
Укажите начальные значения параметров , , , и и укажите, что все параметры должны быть оценены. Обратите внимание, что при оценке моделей черного ящика не существовало требования указывать начальные догадки для параметров sysARMA и sysNLARX.
Parameters = struct('Name',{'Survival factor, predators' 'Death factor, predators' ... 'Survival factor, preys' 'Death factor, preys' ... 'Crowding factor, preys'}, ... 'Unit',{'1/year' '1/year' '1/year' '1/year' '1/year'}, ... 'Value',{-1.1 0.9 1.1 0.9 0.2}, ... 'Minimum',{-Inf -Inf -Inf -Inf -Inf}, ... 'Maximum',{Inf Inf Inf Inf Inf}, ... 'Fixed',{false false false false false});
Аналогично, укажите начальные состояния модели и укажите, что оба начальных состояния должны быть оценены.
InitialStates = struct('Name',{'Predator population' 'Prey population'}, ... 'Unit',{'Size (thousands)' 'Size (thousands)'}, ... 'Value',{1.8 1.8}, ... 'Minimum', {0 0}, ... 'Maximum',{Inf Inf}, ... 'Fixed',{false false});
Укажите модель как систему непрерывного времени.
Ts = 0;
Создайте нелинейную серую модель с заданными структурой, параметрами и состояниями.
sysGrey = idnlgrey(FileName,Order,Parameters,InitialStates,Ts,'TimeUnit','years');
Оцените параметры модели.
sysGrey = nlgreyest(ze,sysGrey); present(sysGrey)
sysGrey =
Continuous-time nonlinear grey-box model defined by 'predprey2_m' (MATLAB file):
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.
Вычислите прогнозируемый результат на 10 шагов вперед для проверки модели.
yhat3 = predict(sysGrey,ze,10);
Прогнозирование выходных данных модели на 100 шагов после данных оценки и вычисление стандартных отклонений выходных данных.
[yf3,x03,sysf3,ysd3] = forecast(sysGrey,ze,100);
Постройте график измеренных, прогнозируемых и прогнозируемых результатов.
t = yf3.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat3.y(:,1),... t,yf3.y(:,1),'m',... t,yf3.y(:,1)+ysd3(:,1),'k--', ... t,yf3.y(:,1)-ysd3(:,1),'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat3.y(:,2),... t,yf3.y(:,2),'m',... t,yf3.y(:,2)+ysd3(:,2),'k--', ... t,yf3.y(:,2)-ysd3(:,2),'k--') legend('Measured','Predicted','Forecasted','Forecast uncertainty (1 sd)') xlabel('Time (years)'); ylabel('Prey population (thousands)');
Графики показывают, что прогнозирование с использованием нелинейной серой модели дало хорошие результаты прогнозирования и низкую неопределенность результатов прогнозирования.
Сравните прогнозируемый ответ, полученный из идентифицированных моделей, sysARMA, sysNLARX, и sysGrey. Первые два являются дискретно-временными моделями и sysGrey является моделью непрерывного времени.
clf
plot(z,yf1,yf2,yf3)
legend({'Measured','Linear ARMA','Nonlinear AR','Nonlinear Grey-Box'})
title('Forecasted Responses')
Прогнозирование с помощью нелинейной модели ARX дало лучшие результаты, чем прогнозирование с помощью линейной модели. Включение знаний о динамике в нелинейную серую модель еще больше повысило надежность модели и, следовательно, точность прогнозирования.
Обратите внимание, что уравнения, используемые в серо-боксовом моделировании, тесно связаны с полиномиальными регрессорами, используемыми нелинейной моделью ARX. Если аппроксимировать производные в управляющих уравнениях разностями первого порядка, вы получите уравнения, подобные:
t-1) * y2 (t-1)
) -s5 * y2 (t-1) 2
Где, s5 - это некоторые параметры, связанные с исходными параметрами p5 и с временем выборки, используемым для дифференцирования. Эти уравнения предполагают, что только 2 регрессора необходимы для первого выхода (y1 (t-1) и * y1 (t-1) * y2 (t-1)) и 3 для второго выхода при построении нелинейной модели ARX. Даже в отсутствие такого предшествующего знания структуры модели линейного регресса, использующие полиномиальные регрессоры, остаются популярным выбором на практике.
Прогнозирование значений с использованием нелинейной серой модели в течение 200 лет.
[yf4,~,~,ysd4] = forecast(sysGrey, ze, 2000);
Постройте график последней части данных (с неопределенностью 1 sd)
t = yf4.SamplingInstants; N = 700:2000; subplot(1,2,1); plot(t(N), yf4.y(N,1), 'm',... t(N), yf4.y(N,1)+ysd4(N,1), 'k--', ... t(N), yf4.y(N,1)-ysd4(N,1), 'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); ax = gca; ax.YLim = [0.8 1]; ax.XLim = [82 212]; subplot(1,2,2); plot(t(N),yf4.y(N,2),'m',... t(N),yf4.y(N,2)+ysd4(N,2),'k--', ... t(N),yf4.y(N,2)-ysd4(N,2),'k--') legend('Forecasted population','Forecast uncertainty (1 sd)') xlabel('Time (years)'); ylabel('Prey population (thousands)'); ax = gca; ax.YLim = [0.9 1.1]; ax.XLim = [82 212];
Сюжет показывает, что популяция хищников, согласно прогнозам, достигнет устойчивого состояния приблизительно 890, а популяция добычи, согласно прогнозам, достигнет 990.