R语言ar函数预测 r语言ar模型参数估计

AR模型的参数的估计

一、AR模型

随机信号 x(n) 由本身的若干次过去值 x(n-k) 和当前的激励值 w(n)线性组合产生：

该模型的系统函数是：

p 是系统阶数，系统函数中只有极点，无零点，也称为全极点模型，系统由于极点的原因，要考虑到系统的稳定性，因而要注意极点的分布位置，用 AR( p ) 来表示。

二、AR模型参数和自相关函数的关系

1、推导过程

2、举例

已知自回归信号模型 **AR（3）**为：

式中 w( n ) 是具有方差 σ² =1 的平稳白噪声，求

a. 求自相关序列

b. 利用a求出的自相关序列估计参数a_k以及输入的白噪声的方差 σ² 大小

可以发现对 AR 模型参数是无失真的估计，因为已知 AR 模型，我们可以得到完全的输出观测值，因而求得的自相关函数没有失真，当然也就可以不失真的估计。

c. 用仿真出的32点观测值 x(n) 直接估计参数a_k以及输入的白噪声的方差 σ² 大小

先求自相关序列

把头4个相关序列值代入矩阵求得估计值：

与真实 AR 模型参数误差为： e₁ ＝ 0.1151， e₂ ＝ 0.1002， e₃ ＝ 0.0498，原因在于我们只有一部分的观测数据，使得自相关序列值与理想的不同。造成的原因比较多，计算机仿真的白噪声由于只有 32 点长，32 点序列的方差不可能刚好等于1。给出一段观测值求 AR 模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

三、Y-W方程的解法——L-D算法

要得到更精确的估计值，就要建立更高阶的 AR 模型，直接用观测 值的自相关序列来求解 Y-W 方程计算量太大。因此把 AR 模型和预测系统联系起来，换个方法来估计 AR 参数。

前向预测误差系统对观测信号起了白化的作用。由于 AR 模型和前向预测误差系统有着密切的关系，两者的系统函数互为倒数， 所以求 AR 模型参数就可以通过求预测误差系统的预测系数来实现。

求预测误差均方值：

代入求得最小均方误差：

也就是 p 阶预测器的预测系数等于 p 阶 AR 模型的参数，由于 e( n) ＝ w( n ) ，所以最小均方预测误差等于白噪声方差，即 E_p [e² (n)] = σ² 。

已知自相关函数， 给定初始值 E 0 = R(0) ，a 0 (0) = 1 ，以及 AR 模型的阶数 p ，就可以按照图 7.7 所示流程图进行估计。

直接利用MATLAB中 [a E]＝aryule(x,p) 函数实现.输入 p 表示要求的阶数，输出 a 表示估计的模型参数，e 表示噪声信号的方差估计。

clear;
close all;
clc;
xn=[0.4282,1.1454,1.5597,1.8994,1.6854,2.3075,2.4679,1.9790,1.6063,1.2804,-0.2083,0.0577,0.0206,0.3572,1.6572,0.7488,1.6666,1.9830,2.6914,1.2521 ,1.8691, 1.6855,0.6242,0.1763,1.3490,0.6955,1.2941,1.0475,0.4319,0.0312,0.5802,-0.6177];
[a E]=aryule(xn,3)

输出

a =

    1.0000   -0.6984   -0.2748    0.0915


E =

    0.4678

假如我们用更高阶的 AR 模型来估计：

[a E]=aryule(xn,8)

a =

  1 至 7 列

    1.0000   -0.6878   -0.3225   -0.1095    0.1661    0.3815   -0.2619

  8 至 9 列

    0.0630   -0.1426


E =

    0.3852

[a E]=aryule(xn,12)

a =

  1 至 7 列

    1.0000   -0.6703   -0.3254   -0.0793    0.1407    0.3676   -0.2451

  8 至 13 列

    0.0483   -0.0912   -0.0522    0.0515    0.0186   -0.0955


E =

    0.3783

四、示例

给出一个信号，利用自相关矩阵求解和L-D算法分别求解估计参数a_k及白噪声方差σ²

clear;
close all;
clc;

[y,fs] = audioread("test.m4a");
y=y';
n=6;
%% ---------------------------------------------------
% 求自相关
r=zeros(1,n);
[c,lags] = xcorr(y);
c=c/length(y);
r=c(1,(length(c)+1)/2:(length(c)+1)/2+n-1);
% 生成自相关矩阵
T = toeplitz(r,r)
T_temp=T(2:end,2:end)
b=-r(1,2:end)';
a1=inv(T_temp)*b;
a1=[1,a1']
E1=r*a1'
%% ---------------------------------------------------
% L-D算法
[a2 E2]=aryule(y,n-1)

结果如下