Этот пример показов, как выполнить многомерное прогнозирование временных рядов данных, измеренных от хищника и популяций жертв в сценарии скопления добычи. Динамика изменения популяции хищника-добычи моделируется с помощью линейных и нелинейных моделей временных рядов. Сравнивают прогнозирующую эффективность этих моделей.
Данные представляют собой двухмерные временные ряды, состоящий из 1-хищных 1-добычных популяций (в тысячах), собираемых 10 раз в год в течение 20 лет. Для получения дополнительной информации о данных смотрите Три Экологические Системы Населения: MATLAB и Файл MEX на C Modeling of Time-Series.
Загрузите данные временных рядов.
load PredPreyCrowdingData z = iddata(y,[],0.1,'TimeUnit','years','Tstart',0);
z является iddata объект, содержащий два выходных сигнала, y1 и y2, которые относятся к хищнику и популяциям жертв, соответственно. The 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);
The 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);The 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 но находится в форме пространства состояний. Симуляция 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) в Statistics and Machine Learning Toolbox™.
Если вы знаете ранее динамику системы, можно подогнать данные оценки с помощью нелинейной модели серого-прямоугольника.
Знания о природе динамики могут помочь улучшить качество модели и, следовательно, точность прогноза. Для динамики хищник-добыча, изменения у хищника (y1) и добыча (y2) население может быть представлена как:
Для получения дополнительной информации о уравнениях смотрите Три Экологические Системы Населения: MATLAB и Файл MEX на C 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. Если вы аппроксимируете производные в управляющих уравнениях различиями первого порядка, вы получите уравнения, подобные:
Где, являются некоторыми параметрами, связанными с исходными параметрами и к шагу расчета, используемой для дифференцирования. Эти уравнения предполагают, что только 2 регрессоры необходимы для первого выхода (y1 (t-1) и * y1 (t-1) * y2 (t-1)) и 3 для второго выхода при построении нелинейной модели ARX. Даже при отсутствии таких предшествующих знаний структуры модели линейного регрессора, использующие полиномиальные регрессоры, остаются популярным выбором на практике.
Прогнозируйте значения с помощью нелинейной модели серого ящика в течение 200 лет.
[yf4,~,~,ysd4] = forecast(sysGrey, ze, 2000);
Постройте график последней части данных (показав 1 неопределенность сд)
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.