Линейная регрессия частотного диапазона

В этом примере показано, как использовать дискретное преобразование Фурье, чтобы создать модель линейной регрессии какое-то время ряд. Временные ряды, используемые в этом примере, являются ежемесячным количеством смертей от несчастного случая в Соединенных Штатах от 1 973 до 1979. Данные публикуются в Броквелле и Дэвисе (2006). Первоисточником является американский Национальный совет по безопасности.

Введите данные. Скопируйте exdata матрица в рабочую область MATLAB®.

exdata = [
        9007        7750        8162        7717        7792        7836
        8106        6981        7306        7461        6957        6892
        8928        8038        8124        7776        7726        7791
        9137        8422        7870        7925        8106        8129
       10017        8714        9387        8634        8890        9115
       10826        9512        9556        8945        9299        9434
       11317       10120       10093       10078       10625       10484
       10744        9823        9620        9179        9302        9827
        9713        8743        8285        8037        8314        9110
        9938        9129        8433        8488        8850        9070
        9161        8710        8160        7874        8265        8633
        8927        8680        8034        8647        8796        9240];

exdata 12 6 матрица. Каждый столбец exdata содержит 12 месяцев данных. Первая строка каждого столбца содержит количество американских смертей от несчастного случая на январь соответствующего года. Последняя строка каждого столбца содержит количество американских смертей от несчастного случая на декабрь соответствующего года.

Измените матрицу данных в 72 1 временные ряды и отобразите данные на графике в течение лет 1973 - 1978.

ts = reshape(exdata,72,1);
years = linspace(1973,1979,72);

plot(years,ts,'o-','MarkerFaceColor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')

Визуальный осмотр данных указывает, что количество смертей от несчастного случая варьируется по периодическому способу. Период колебания, кажется, составляет примерно 1 год (12 месяцев). Периодическая природа данных предполагает, что соответствующая модель может быть

где $$\mu$$ полное среднее значение, $$N$$ длина временных рядов, и $$\epsilon (n)$$ белая шумовая последовательность независимого политика и тождественно распределенных Гауссовых случайных переменных (IID) с нулевым средним значением и некоторым отклонением. Аддитивный шумовой термин составляет случайность, свойственную от данных. Параметры модели являются полным средним значением и амплитудами косинусов и синусов. Модель линейна в параметрах.

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

Выполните спектральный анализ данных, чтобы показать, какие частоты значительно способствуют изменчивости в данных. Поскольку полное среднее значение сигнала - приблизительно 9 000 и пропорционально преобразованию Фурье на 0 частотах, вычтите среднее значение до спектрального анализа. Это уменьшает большую величину коэффициент Фурье на 0 частотах и делает любые значительные колебания легче обнаружить. Частоты в преобразовании Фурье расположены с интервалами в интервале, который является обратной величиной длины временных рядов, 1/72. Производя данные ежемесячно, самая высокая частота в спектральном анализе составляет 1 месяц цикла/2. В этом случае удобно посмотреть на спектральный анализ в терминах циклов/год, так масштабируйте частоты соответственно для визуализации.

tsdft = fft(ts-mean(ts));
freq = 0:1/72:1/2;

plot(freq.*12,abs(tsdft(1:length(ts)/2+1)),'o-', ...
    'MarkerFaceColor','auto')
xlabel('Cycles/Year')
ylabel('Magnitude')
ax = gca;
ax.XTick = [1/6 1 2 3 4 5 6];

На основе величин частота 1 месяца цикла/12 является старшим значащим колебанием в данных. Величина в 1 месяц цикла/12 является более двух раз столь же большой как любая другая величина. Однако спектральный анализ показывает, что существуют также другие периодические компоненты в данных. Например, кажется, существуют периодические компоненты в гармониках (целочисленные множители) 1 месяца цикла/12. Также, кажется, существует периодический компонент с периодом 1 месяца цикла/72.

На основе спектрального анализа данных подбирайте простую модель линейной регрессии использование косинуса и термина синуса с частотой старшего значащего компонента: 1 цикл/год (1 месяц цикла/12).

Определите интервал частоты в дискретном преобразовании Фурье, которое соответствует 1 месяцу цикла/12. Поскольку частоты расположены с интервалами в 1/72, и первый интервал соответствует 0 частотам, правильный интервал является 72/12+1. Это - интервал частоты положительной частоты. Необходимо также включать интервал частоты, соответствующий отрицательной частоте:-1 месяц цикла/12. С индексацией MATLAB интервал частоты отрицательной частоты является 72-72/12+1.

Создайте 72 1 нулевой вектор. Заполните соответствующие элементы вектора с коэффициентами Фурье, соответствующими положительной и отрицательной частоте 1 месяца цикла/12. Инвертируйте преобразование Фурье и добавьте полное среднее значение, чтобы получить подгонку к данным о смерти от несчастного случая.

freqbin = 72/12;
freqbins = [freqbin 72-freqbin]+1;
tsfit = zeros(72,1);
tsfit(freqbins) = tsdft(freqbins);
tsfit = ifft(tsfit);
mu = mean(ts);
tsfit = mu+tsfit;

Отобразите исходные данные на графике наряду с подходящим рядом с помощью двух коэффициентов Фурье.

plot(years,ts,'o-','MarkerFaceColor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')
hold on
plot(years,tsfit,'linewidth',2)
legend('Data','Fitted Model')
hold off

Подобранная модель, кажется, получает общую периодическую природу данных и поддерживает первоначальное заключение, что данные колеблются с циклом 1 года.

Чтобы оценить, как соответственно одна частота 1 месяца цикла/12 составляет наблюдаемые временные ряды, сформируйте остаточные значения. Если остаточные значения напоминают белую шумовую последовательность, простая линейная модель с одной частотой соответственно смоделировала временные ряды.

Чтобы оценить остаточные значения, используйте последовательность автокорреляции с 95%-доверительными-интервалами для белого шума.

resid = ts-tsfit;
[xc,lags] = xcorr(resid,50,'coeff');

stem(lags(51:end),xc(51:end),'filled')
hold on
lconf = -1.96*ones(51,1)/sqrt(72);
uconf = 1.96*ones(51,1)/sqrt(72);
plot(lags(51:end),lconf,'r')
plot(lags(51:end),uconf,'r')
xlabel('Lag')
ylabel('Correlation Coefficient')
title('Autocorrelation of Residuals')
hold off

Значения автокорреляции выходят за пределы 95% доверительных границ во многих задержках. Не кажется, что остаточные значения являются белым шумом. Заключение состоит в том, что простая линейная модель с одним синусоидальным компонентом не составляет все колебания в количестве смертей от несчастного случая. Это ожидается, потому что спектральный анализ показал дополнительные периодические компоненты в дополнение к доминирующему колебанию. Создание модели, которая включает дополнительные периодические условия, обозначенные спектральным анализом, улучшит подгонку и побелит остаточные значения.

Подберите модель, которая состоит из трех самых больших содействующих величин Фурье. Поскольку необходимо сохранить коэффициенты Фурье, соответствующие и отрицательным и положительным частотам, сохраните самые большие 6 индексов.

tsfit2dft = zeros(72,1);
[Y,I] = sort(abs(tsdft),'descend');
indices = I(1:6);
tsfit2dft(indices) = tsdft(indices);

Продемонстрируйте, что сохранение только 6 из 72 коэффициентов Фурье (3 частоты) сохраняет большую часть энергии сигнала. Во-первых, продемонстрируйте, что сохранение всех коэффициентов Фурье дает к энергетической эквивалентности между исходным сигналом и преобразованием Фурье.

norm(1/sqrt(72)*tsdft,2)/norm(ts-mean(ts),2)

ans = 1.0000

Отношение равняется 1. Теперь исследуйте энергетическое отношение, где только 3 частоты сохраняются.

norm(1/sqrt(72)*tsfit2dft,2)/norm(ts-mean(ts),2)

ans = 0.8991

Почти 90% энергии сохраняются. Эквивалентно, 90% отклонения временных рядов составляются 3 частотными составляющими.

Составьте мнение данных на основе 3 частотных составляющих. Сравните исходные данные, модель с одной частотой и модель с 3 частотами.

tsfit2 = mu+ifft(tsfit2dft,'symmetric');

plot(years,ts,'o-','markerfacecolor','auto')
xlabel('Year')
ylabel('Number of Accidental Deaths')
hold on
plot(years,tsfit,'linewidth',2)
plot(years,tsfit2,'linewidth',2)
legend('Data','1 Frequency','3 Frequencies')
hold off

Используя 3 частоты улучшил подгонку к исходному сигналу. Вы видите это путем исследования автокорреляции остаточных значений модели с 3 частотами.

resid = ts-tsfit2;
[xc,lags] = xcorr(resid,50,'coeff');

stem(lags(51:end),xc(51:end),'filled')
hold on
lconf = -1.96*ones(51,1)/sqrt(72);
uconf = 1.96*ones(51,1)/sqrt(72);
plot(lags(51:end),lconf,'r')
plot(lags(51:end),uconf,'r')
xlabel('Lag')
ylabel('Correlation Coefficient')
title('Autocorrelation of Residuals')
hold off

Используя 3 частоты привел к остаточным значениям, которые более тесно аппроксимируют белый шумовой процесс.

Продемонстрируйте, что значения параметров, полученные из преобразования Фурье, эквивалентны модели линейной регрессии временного интервала. Найдите оценки наименьших квадратов для полного среднего значения, амплитуд косинуса и амплитуд синуса для этих трех частот путем формирования матрицы проекта и решения нормальных уравнений. Сравните подходящие временные ряды с полученным из преобразования Фурье.

X = ones(72,7);
X(:,2) = cos(2*pi/72*(0:71))';
X(:,3) = sin(2*pi/72*(0:71))';
X(:,4) = cos(2*pi*6/72*(0:71))';
X(:,5) = sin(2*pi*6/72*(0:71))';
X(:,6) = cos(2*pi*12/72*(0:71))';
X(:,7) = sin(2*pi*12/72*(0:71))';
beta = X\ts;
tsfit_lm = X*beta;
max(abs(tsfit_lm-tsfit2))

ans = 1.2733e-11

Эти два метода приводят к идентичным результатам. Максимальное абсолютное значение различия между этими двумя формами волны находится порядка 10-12. В этом случае подход частотного диапазона был легче, чем эквивалентный подход временного интервала. Вы естественно используете спектральный анализ, чтобы визуально смотреть, какие колебания присутствуют в данных. От того шага это просто в использовании коэффициенты Фурье, чтобы создать модель для сигнала, состоящего из суммы косинусы и синусы.

Для получения дополнительной информации о спектральном анализе во временных рядах и эквивалентности с регрессией временного интервала смотрите (Шумвей и Стоффер, 2006).

В то время как спектральный анализ может ответить, какие периодические компоненты значительно способствуют изменчивости данных, это не объясняет, почему те компоненты присутствуют. Если вы исследуете эти данные тесно, вы видите, что минимальные значения в 12-месячном цикле имеют тенденцию происходить в феврале, в то время как максимальные значения происходят в июле. Вероятное объяснение этих данных состоит в том, что люди естественно более активны летом, чем зимой. К сожалению, в результате этого увеличенного действия, существует увеличенная вероятность вхождения несчастных случаев со смертельным исходом.

Ссылки

Броквелл, Питер Дж. и Ричард А. Дэвис. Временные ряды: теория и методы. Нью-Йорк: Спрингер, 2006.

Shumway, Роберт Х. и Дэвид С. Стоффер. Анализ временных рядов и его приложения с примерами R. Нью-Йорк: Спрингер, 2006.

Смотрите также

fft | ifft | xcorr