Этот пример показывает, как выполнить многомерное прогнозирование временных рядов данных, измеренных от хищника и населения добычи в сценарии давки добычи. Движущие силы изменения генеральной совокупности добычи хищника моделируются с помощью линейных и нелинейных моделей временных рядов. Прогнозирование производительности этих моделей сравнено.
Данные являются двумерными временными рядами, состоящими из населения с 1 добычей с 1 хищником (в тысячах), собрался 10 раз в год в течение 20 лет. Для получения дополнительной информации о данных, см. Три Экологических Системы Генеральной совокупности: MATLAB и Моделирование Файла MEX на C Timeseries.
Загрузите данные временных рядов.
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, но находится в форме пространства состояний. Симуляция 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 лет. Это указывает, что динамика изменения генеральной совокупности может быть нелинейной.
Соответствуйте нелинейной модели черного ящика к данным об оценке. Вы не требуете предварительных знаний об уравнениях, управляющих данными об оценке. Будет оценена линейная в регрессоре форма модели Nonlinear ARX.
Создайте нелинейную модель ARX с 2 выходными параметрами и никакими входными параметрами.
sysNLARX = idnlarx([1 1;1 1],[],'Ts',0.1,'TimeUnit','years');
sysNLARX является первым порядком нелинейная модель ARX, которая не использует нелинейной функции; это предсказывает образцовый ответ с помощью взвешенной суммы его регрессоров первого порядка.
getreg(sysNLARX)
Regressors:
For output 1:
y1(t-1)
y2(t-1)
For output 2:
y1(t-1)
y2(t-1)
Чтобы ввести функцию нелинейности, добавьте полиномиальные регрессоры в модель.
Создайте регрессоры до степени 2 и включайте перекрестные условия (продукты стандартных упомянутых выше регрессоров). Добавьте те регрессоры в модель как пользовательские регрессоры.
R = polyreg(sysNLARX,'MaxPower',2,'CrossTerm','on'); sysNLARX.CustomRegressors = R; getreg(sysNLARX)
Regressors:
For output 1:
y1(t-1)
y2(t-1)
y1(t-1).^2
y1(t-1).*y2(t-1)
y2(t-1).^2
For output 2:
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 ARX model with 2 outputs and 0 inputs
Inputs:
Outputs: y1, y2
Standard regressors corresponding to the orders:
na = [1 1; 1 1]
nb = []
nk = []
Custom regressors:
For output 1:
y1(t-1).^2
y1(t-1).*y2(t-1)
y2(t-1).^2
For output 2:
y1(t-1).^2
y1(t-1).*y2(t-1)
y2(t-1).^2
Nonlinear regressors:
For output 1:
none
For output 2:
none
Model output is 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 Timeseries.
Создайте нелинейное серое поле, основанное на модели на этих уравнениях.
Задайте файл, описывающий образцовую структуру для системы добычи хищника. Файл задает производные состояния и образцовые выходные параметры как функция времени, состояний, входных параметров и параметров модели. Эти два выходных параметров (хищник и население добычи) выбраны в качестве состояний, чтобы вывести нелинейное описание пространства состояний динамики.
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 дало лучшие результаты, чем прогнозирование с линейной моделью. Включение знания динамики в нелинейной модели серого поля далее улучшило надежность модели и поэтому точности прогнозирования.
Обратите внимание на то, что уравнения, используемые в моделировании серого поля, тесно связаны с полиномиальными регрессорами, используемыми моделью Nonlinear ARX. Если вы аппроксимируете производные в управляющих уравнениях различиями первого порядка, вы получите уравнения, подобные:
Где, некоторые параметры, связанные с исходными параметрами и к шагу расчета используется для дифференцирования. Эти уравнения предполагают, что только 2 регрессора необходимы для первого вывода (y1 (t-1) и *y1 (t-1) *y2 (t-1)) и 3 для второго вывода при построении модели Nonlinear 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.