Проблемы с привязкой данных к линейной модели

Question

Предполагая бесшумный процесс AR (1) y(t)= a*y(t-1). У меня есть следующие концептуальные вопросы и буду рад разъяснению.

Q1 - Несоответствие между математической формулировкой и реализацией. Математическая формулировка модели AR представлена ​​в виде y(t) = - summmation over i=1 to p[a*y(t-p)] + eta(t) где p = порядок модели и eta(t) - белый гауссовский шум. Но при оценке коэффициентов с использованием любого метода, такого как arburg() или наименьшего квадрата, мы просто вызываем эту функцию. Я не знаю, будет ли явно добавлен белый гауссовский шум. Затем, когда мы разрешаем уравнение AR с оцененными коэффициентами, я видел, что отрицательный знак не учитывается, и добавленный шум не добавляется.

Какое правильное представление модели AR и как найти средние коэффициенты по k числу испытаний, когда у меня есть только один образец из 1000 точек данных?

Q2 - Кодирование проблемы в Как имитировать fitted_data для к числу испытаний, а затем найти остатки - я приспособил данные «данные», полученных от неизвестной системы и получил коэффициент по

load('data.txt'); 

for trials = 1:10 

    model = ar(data,1,'ls'); 
    original_data=data; 

    fitted_data(i)=coeff1*data(i-1); % **OR** 
    data(i)=coeff1*data(i-1); 

    fitted_data=data; 

    residual= original_data - fitted_data; 
    plot(original_data,'r'); hold on; plot(fitted_data); 

end 

При расчете остаточного - это данные, полученные как описано выше, путем решения уравнения AR с полученными коэффициентами? У Matlab есть функция для этого, но я хотел сделать свой собственный. Итак, после определения коэффициентов из исходных данных, как мне решить? Указанное выше кодирование неверно. Прикреплен график исходных данных и встроенных_данных.

Answer 1

Модели AR-типа могут служить целым рядом целей, включая линейное предсказание, линейное предсказательное кодирование, фильтрующий шум. Eta (t) не то, что мы заинтересованы в сохранении, а часть точки алгоритмов заключается в том, чтобы удалить их влияние в любой степени, если искать постоянные шаблоны в данных.

У меня есть учебники, которые в контексте линейного предсказания не включают отрицательный знак, включенный в ваше выражение до суммы. С другой стороны, функция от Matlab lpc делает:

Xp(n) = -A(2)*X(n-1) - A(3)*X(n-2) - ... - A(N+1)*X(n-N) 

Я рекомендую вам взглянуть на функции lpc, если вы еще не сделали, и в examples из документации, такой как:

randn('state',0); 
noise = randn(50000,1); % Normalized white Gaussian noise 
x = filter(1,[1 1/2 1/3 1/4],noise); 
x = x(45904:50000); 
% Compute the predictor coefficients, estimated signal, prediction error, and autocorrelation sequence of the prediction error: 
p = lpc(x,3); 
est_x = filter([0 -p(2:end)],1,x); % Estimated signal 
e = x - est_x;      % Prediction error 
[acs,lags] = xcorr(e,'coeff');  % ACS of prediction error 

В Оценочная оценка x рассчитывается как est_x. Обратите внимание, как пример использует filter. снова Цитируя Matlab док, «является„Прямая форма II Транспозиция“реализация стандартного разностного уравнения:

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) 
         - a(2)*y(n-1) - ... - a(na+1)*y(n-na) 

, что означает, что в предыдущем примере est_x(n) вычисляется как

est_x(n) = -p(2)*x(n-1) -p(3)*x(n-2) -p(4)*x(n-3) 

который что вы ожидаете

Edit:

Как Rega rds функция ar, matlab documentation объясняет, что выходные коэффициенты имеют то же значение, что и в описанном выше сценарии lp.

Правильный способ оценить выход модели AR является вычисление

data_armod(i)= -coeff(2)*data(i-1) -coeff(3)*data(i-2) -coeff(4)*data(i-3) 

где коэфф является матрица коэффициентов возвращается с

model = ar(data,3,'ls'); 
coeff = model.a; 

Answer 2

Если модель просто y(n)= a*y(n-1) со скалярным a, то вот решение.

y = randn(10, 1); 
a = y(1 : end - 1) \ y(2 : end); 

y_estim = y * a; 
residual = y - y_estim; 

Конечно, вы должны разделить данные на поезд-тест, а также применять a на тестовых данных. Вы можете обобщить этот подход к y(n)= a*y(n-1) + b*y(n-2) и т.д.

Обратите внимание, что \ представляет mldivide() функцию: mldivide

Edit:

% model: y[n] = c + a*y(n-1) + b*y(n-2) +...+z*y(n-n_order) 
n_order = 3; 
allow_offset = true; % alows c in the model 

% train 
y_train = randn(20,1); % from your data 
[y_in, y_out] = shifted_input(y_train, n_order, allow_offset); 
a = y_in \ y_out; 

% now test 
y_test = randn(20,1); % from your data 
[y_in, y_out] = shifted_input(y_test, n_order, allow_offset); 
y_estim = y_in * a; % same a 
residual = y_out - y_estim; 

здесь является shifted_input():

function [y_in, y_out] = shifted_input(y, n_order, allow_offset) 
y_out = y(n_order + 1 : end); 
n_rows = size(y, 1) - n_order; 
y_in = nan(n_rows, n_order); 
for k = 1 : n_order  
    y_in(:, k) = y(1 : n_rows); 
    y = circshift(y, -1); 
end 
if allow_offset 
    y_in = [y_in, ones(n_rows, 1)]; 
end 
return