Наносим линейный тренд на график временного хода в MATLAB

Сегодня мы нанесём линию тренда с доверительными интервалами на график временного ряда.

Работать мы будем с теми же данными (среднемесячные значения индекса NAO), с которыми работали в предыдущей заметке про добавление дат на график временного хода.

nao = load('norm.nao.monthly.b5001.current.ascii.txt');

Для того, чтобы тренд был явно выражен мы выбрали только последние несколько лет (да-да, это тот самый случай, когда учёные манипулируют с данными, никогда так не делайте в своей научной работе):

dates = datetime(nao(750:end,1),nao(750:end,2),15);
nao=nao(750:end,3);

Построим график для того, чтобы представлять, с чем имеем дело (параметр 'k' делает линию чёрной, а 'LineWidth' устанавливает толщину линии):

plot(dates, nao,'k','LineWidth',2)
Индекс NAO

Тренд здесь хорошо заметен. Если вам нужно удалить линейный тренд из данных, то сделать это можно при помощи простой функции detrend:

plot(dates, detrend(nao),'k','LineWidth',2)
удалить тренд MATLAB

Если же вы хотите не удалить тренд, а наоборот выделить его из данных, то нужно будет использовать другую функцию, или функции, так как выделить тренд в MATLAB можно разными способами. 

Первый способ - через функцию polyfit, которая имеет три входных аргумента:

p = polyfit(x, y, n)

Эта функция аппроксимирует функцию y(x) полиномом p(x) степени n, с помощью метода наименьших квадратов. Т.е. y - это значения нашего индекса NAO, которые изменяются в зависимости от времени (x). Поскольку нас интересует линейный тренд, то нам нужен полином первой степени, значит n = 1. Насколько я понял, функция polyfit пока не умеет работать с временными переменными (datetime), поэтому в качестве аргумента x мы подставим вектор с номерами месяцев: [1:length(nao)]. Что ж, давайте пробовать:

p = polyfit([1:length(nao)],nao,1)
Error using polyfit (line 44)
X and Y vectors must be the same size.

Как видите, возникла ошибка из-за того, что функция polyfit не любит, когда одна переменная имеет размерность 78х1, а другая 1х78. Не страшно, транспонируем одну из переменных и у нас всё получится:

p = polyfit([1:length(nao)]',nao,1)

p =

    0.0162   -0.4251

На выходе мы получили два числа, которые являются коэффициентами полинома и расположены в убывающем порядке, т.е. грубо говоря, уравнение нашей прямой выглядит так:

y  = p(2)+p(1)x

Значит наш индекс NAO в выбранном промежутке времени растёт в среднем на 0.0162 в месяц. 

По полученному уравнению не сложно рассчитать значения для построения линии тренда, однако трудится не стоит, в MATLAB для этого есть специальная функция polyval:

plot(dates,nao,'k','LineWidth',2)
hold on
plot(dates,polyval(p,[1:length(nao)]),':k','LineWidth',2)
Построить тренд MATLAB

В целом, мы справились с задачей, однако скорее всего вам понадобится дополнить этот график доверительными интервалами. И для этого нам нужно указать ещё один выходной аргумент S в функции polyfit, для того, чтобы в него записалась информация, необходимая для дальнейшего расчёта доверительных интервалов:

[p,S] = polyfit([1:length(nao)]',nao,1);

Теперь осталось рассчитать значения линий доверительных интервалов при помощи функции polyconf:

[Y,DELTA] = polyconf(p,[1:length(nao)]',S,'predopt','curve');

Не забудьте выставить параметр 'predopt' на значение 'curve', для того, чтобы доверительные интервалы рассчитывались именно для нашей кривой (по умолчанию рассчитываются спрогнозированные интервалы для новых значений x). Кроме того, следует знать, что по умолчанию рассчитывается 95% доверительный интервал и, если вы хотите его изменить, то следует поменять параметр 'alpha', например на 0.01 для 99% доверительного интервала (т.е. 100*(1-alpha))

На выходе функции polyconf мы будем иметь две переменных. В переменной Y будут содержаться значения для прямой тренда (кстати говоря, мы можем их использовать вместо функции polyval при построении), а в переменной DELTA - отклонения от этой прямой. Имея эти данные, построить доверительные интервалы не составит никакого труда:

plot(dates,nao,'k','LineWidth',2)
hold on
plot(dates,polyval(p,[1:length(nao)]),'--k','LineWidth',2)
hold on
plot(dates,Y+DELTA,':k','LineWidth',2)
hold on
plot(dates,Y-DELTA,':k','LineWidth',2)
тренд с доверительными интервалами при помощи polyfit MATLAB

Другой распространённый способ выделения тренда - при помощи регрессии. Для построения регрессии в MATLAB в разное время создавали разные функции и, например, можно использовать функцию с говорящим названием regression. Однако, начиная с 2013 года появилась более удобная для использования, функция fitlm, которой мы и воспользуемся. 

mdl = fitlm([1:length(nao)],nao)

mdl = 


Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate       SE         tStat      pValue  
                   ________    _________    _______    _________

    (Intercept)    -0.42515      0.25759    -1.6505      0.10297
    x1             0.016168    0.0056655     2.8537    0.0055653


Number of observations: 78, Error degrees of freedom: 76
Root Mean Squared Error: 1.13
R-squared: 0.0968,  Adjusted R-Squared 0.0849
F-statistic vs. constant model: 8.14, p-value = 0.00557

Как видно, мы получили те же самые значения для линии тренда. Однако, следует заметить, что функция polyfit сделает это быстрее. Кстати говоря, если вы не знаете, увидеть время для выполнения функции можно, написав перед ней @(), а затем использовать функцию timeit:

mdl = @() fitlm([1:length(nao)],nao);
timeit(mdl)

ans =

    0.0079

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = @() polyfit([1:length(nao)]',nao,1);
>> timeit(p)

ans =

   4.5270e-04
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timeit(p)<timeit(mdl)

ans =

  logical

   1

Но вернёмся к результатам нашей регрессии. Для того, чтобы получить доверительные интервалы для коэффициентов нашего тренда, нужно воспользоваться функцией coefCI:

ci = coefCI(mdl)

ci =

   -0.9382    0.0879
    0.0049    0.0275

В итоге мы получили две строки, в первой строке содержатся доверительные интервалы для свободного члена уравнения регрессии, а вторая строка - интервалы для трендовой составляющей.

Если мы хотим построить доверительные интервалы для линии тренда, чтобы построить их на графике, как мы это делали ранее, то воспользуемся функцией predict, которая по нашему уравнению регрессии mdl вычисляет спрогнозированные значения ypred и значения для доверительных интервалов yci (два вектроа - для нижней и верхней границы). В нашем случае в качестве входных значений мы используем те же самые значения, по которым расчитывалось уравнение регрессии:

[ypred,yci] = predict(mdl,[1:length(nao)]');

Теперь построим:

plot(dates,nao,'k','LineWidth',2)
hold on
plot(dates,ypred,'--k','LineWidth',2)
hold on
plot(dates,yci(:,1),':k','LineWidth',2)
hold on
plot(dates,yci(:,2),':k','LineWidth',2)
тренд с доверительными интервалами при помощи fitlm MATLAB

Кстати говоря, вам не обязательно строить всё вручную. Если вы хотите быстро проверить результат и вас устроит график, который сделает MATLAB  по умолчанию, то можно просто написать:

plot(mdl)

и получить такую картинку:

fitlm регрессия в MATLAB
 telegram

Чтобы не пропустить новые материалы с рецептами по работе с океанологическими данными, подпишитесь на канал в Telegram: https://t.me/koldunovaleksey

P.S. Пока писал пост для Яндекс Дзена, сочинил универсальную функцию, на случай, если вам очень быстро нужно построить тренд на графике:

function plot_data_and_trend (data)
plot([1:length(data)],polyval(polyfit([1:length(data)]',data,1),[1:length(data)]),[1:length(data)],data);