В этом примере показано, как выполнить многомерное прогнозирование временных рядов данных, измеренных от хищника и популяций добычи в сценарии давки добычи. Движущие силы изменения населения добычи хищника моделируются с помощью линейных и нелинейных моделей временных рядов. Прогнозирование эффективности этих моделей сравнено.
Данные являются двумерными временными рядами, состоящими из популяций с 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)
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 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.