Остаточный анализ с автокорреляцией
Этот пример показывает, как использовать автокорреляцию с доверительным интервалом, чтобы анализировать невязки подгонки наименьших квадратов к шумным данным. Невязки являются различиями между подобранной моделью и данными. В модели шума сигнала-плюс-белый, если у вас есть подходящий вариант для сигнала, невязки должны быть белым шумом.
Создайте шумный набор данных, состоящий из полинома 1-го порядка (прямая линия) в аддитивном белом Гауссовом шуме. Аддитивный шум является последовательностью некоррелированых случайных переменных после N (0,1) распределение. Это означает, что все случайные переменные имеют средний нуль и модульное отклонение. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.
x = -3:0.01:3;
rng default
y = 2*x+randn(size(x));
plot(x,y)
Используйте polyfit, чтобы найти строку наименьших квадратов для шумных данных. Отобразите исходные данные на графике наряду с подгонкой наименьших квадратов.
coeffs = polyfit(x,y,1); yfit = coeffs(2)+coeffs(1)*x; plot(x,y) hold on plot(x,yfit,'linewidth',2)
Найдите невязки. Получите последовательность автокорреляции невязок, чтобы отстать 50.
residuals = y - yfit;
[xc,lags] = xcorr(residuals,50,'coeff');Когда вы осматриваете последовательность автокорреляции, вы хотите определить, существует ли доказательство автокорреляции. Другими словами, вы хотите определить, похожа ли демонстрационная последовательность автокорреляции на последовательность автокорреляции белого шума. Если последовательность автокорреляции невязок похожа на автокорреляцию белого шумового процесса, вы уверены, что ни один из сигнала не вышел из вашей подгонки и закончился в невязках. В этом примере используйте 99%-доверительных-интервалов. Чтобы создать доверительный интервал, необходимо знать распределение демонстрационных значений автокорреляции. Также необходимо найти критические значения на соответствующем распределении, между которым находятся 0.99 из вероятности. Поскольку распределение в этом случае является Гауссовым, можно использовать дополнительную обратную функцию ошибок, erfcinv. Отношение между этой функцией и инверсией Гауссовой кумулятивной функции распределения описано на странице с описанием для erfcinv.
Найдите критическое значение для 99%-доверительных-интервалов. Используйте критическое значение, чтобы создать более низкие и верхние доверительные границы.
conf99 = sqrt(2)*erfcinv(2*.01/2); lconf = -conf99/sqrt(length(x)); upconf = conf99/sqrt(length(x));
Постройте последовательность автокорреляции наряду с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.03 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
Кроме в нулевой задержке, демонстрационные значения автокорреляции лежат в этих 99%-доверительных-границах для автокорреляции белой шумовой последовательности. От этого можно прийти к заключению, что невязки являются белым шумом. А именно, вы не можете отклонить это, невязки являются реализацией белого шумового процесса.
Создайте сигнал, состоящий из синусоиды плюс шум. Данные выбираются на уровне 1 кГц. Частота синусоиды составляет 100 Гц. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
rng default
x = cos(2*pi*100*t)+randn(size(t));Используйте дискретное преобразование Фурье (DFT), чтобы получить подгонку наименьших квадратов к синусоиде на уровне 100 Гц. Оценка наименьших квадратов амплитуды равняется 2 / N времена коэффициент ДПФ, соответствующий 100 Гц, где N является длиной сигнала. Действительная часть является амплитудой косинуса на уровне 100 Гц, и мнимая часть является амплитудой синуса на уровне 100 Гц. Подгонка наименьших квадратов является суммой косинуса и синуса с правильной амплитудой. В этом примере интервал ДПФ 101 соответствует 100 Гц.
xdft = fft(x); ampest = 2/length(x)*xdft(101); xfit = real(ampest)*cos(2*pi*100*t)+imag(ampest)*sin(2*pi*100*t); figure plot(t,x) hold on plot(t,xfit,'linewidth',2) axis([0 0.30 -4 4]) xlabel('Seconds') ylabel('Amplitude')
Найдите невязки и определите демонстрационную последовательность автокорреляции, чтобы отстать 50.
residuals = x-xfit;
[xc,lags] = xcorr(residuals,50,'coeff');Постройте последовательность автокорреляции с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.03 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
Снова, вы видите, что кроме в нулевой задержке, демонстрационные значения автокорреляции лежат в этих 99%-доверительных-границах для автокорреляции белой шумовой последовательности. От этого можно прийти к заключению, что невязки являются белым шумом. А именно, вы не можете отклонить это, невязки являются реализацией белого шумового процесса.
Наконец, добавьте другую синусоиду с частотой 200 Гц и амплитудой 3/4. Соответствуйте только синусоиде на уровне 100 Гц и найдите демонстрационную автокорреляцию невязок.
x = x+3/4*sin(2*pi*200*t);
xdft = fft(x);
ampest = 2/length(x)*xdft(101);
xfit = real(ampest)*cos(2*pi*100*t)+imag(ampest)*sin(2*pi*100*t);
residuals = x-xfit;
[xc,lags] = xcorr(residuals,50,'coeff');Постройте демонстрационную автокорреляцию наряду с этими 99%-доверительными-интервалами.
figure stem(lags,xc,'filled') ylim([lconf-0.12 1.05]) hold on plot(lags,lconf*ones(size(lags)),'r','linewidth',2) plot(lags,upconf*ones(size(lags)),'r','linewidth',2) title('Sample Autocorrelation with 99% Confidence Intervals')
В этом случае значения автокорреляции ясно превышают эти 99%-доверительных-границ для белой шумовой автокорреляции во многих задержках. Здесь можно отклонить гипотезу, что невязки являются белой шумовой последовательностью. Значение - то, что модель не составляла весь сигнал, и поэтому невязки состоят из сигнала плюс шум.