将数据拟合到线性模型的问题

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 模型的正确表示？当我只有一个包含 1000 个数据点的样本时，如何找到 k 次试验的平均系数？

Q2 - 如何为 k 次试验模拟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 方程获得了上述 fit_data？ Matlab 有这样做的功能，但我想自己做。那么，在从原始数据中找到系数后，我该如何解决？上面的编码不正确。附上原始数据和fitted_data的图。

Answer 1

如果您的模型只是 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;

当然，你应该将数据分成train-test，并在测试数据上应用a。您可以将此方法推广到y(n)= a*y(n-1) + b*y(n-2)等。

注意\代表mldivide()函数：mldivide

编辑：

% 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

Answer 2

AR 类型的模型可以用于多种用途，包括线性预测、线性预测编码、过滤噪声。 eta(t) 不是我们有兴趣保留的东西，而是算法的一部分是通过寻找数据中的持久模式来尽可能地消除它们的影响。

我有一些教科书，在线性预测的上下文中，不包括在求和之前的表达式中包含的负号。另一方面Matlab的函数lpcdoes：

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 文档，filter(b,a,x)“是标准差分方程的“直接形式 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)

这是你所期望的！

编辑：

关于函数ar，matlab documentation 解释了输出系数与上面讨论的lp 场景中的含义相同。

评估AR模型输出的正确方法是计算

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

其中 coeff 是使用

返回的系数矩阵

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