随机信号的参数建模法

随机信号的参数建模法

在对语音信号进行编码时，对语言产生的物理过程进行建模，编码器根据输入信号计算模型参数，然后对模型参数进行编码。当解码器收到模型参数后，再利用数学模型重建原始数据。这种利用数学模型和参数计算的编码方式称为信源模型编码技术。可以用下图表示过程：

对于平稳随机信号，有三种常用的线性模型分别是：

• AR模型（自回归模型Auto-regression model）
• MA模型（滑动平均模型Moving average model）
• ARMA模型（自回归滑移平均模型Auto-regression-Moving average model）

由于任何平稳的ARMA模型或者MA模型均可以用无限阶或阶数足够的AR模型去近似。因此本文着重介绍AR模型的基本原理和方法。

MA模型

随机信号$x(n)$由当前的激励$w(n)$和若干次过去的激励$w(n-k)$线性组合产生:
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)$$
对上述表达式进行拉普拉斯变换可以得到该模型的系统函数为：
$$H(z)=\frac{X(z)}{W(z)}=\sum _{k=0}^q{b}_k{z}^{-k}$$

q表示系统阶数，系统函数只有零点，没有极点，所以该系统一定是稳定的系统，也称为全零点模型，用$MA(q)$来表示。

AR模型

随机信号$x(n)$由本身的若干次过去值$x(n-k)$和当前的激励值$w(n)$线性组合产生：
$$x(n)=w(n)-\sum _{k=1}^p{a}_{k}x(n-k)$$

该模型的系统函数为：
$$H(z)=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$

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

ARMA模型

ARMA是AR与MA模型的结合：
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)-\sum _{k=1}^p{a}_{k}x(n-k)$$
该模型的系统函数是：
$$H(z)=\frac{{\sum }_{k=0}^q{b}_k{z}^{-k}}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$

它既有零点又有极点，所以也称极零点模型，要考虑极零点的分布位置，保证系统的稳定，用$ARMR(p，q)$表示。

在随机信号时域分析中，提出了许多数学模型用来由已知在最大不确定原则下预测将来值，其优点是只需要很少的已知值。但是它不能用在信号是确定性的场合，在确定信号的情况下，信号是由确定的数学方程预测的

AR模型参数的估计

随机信号的建模法最近在信号处理中应用相当普遍。。实际中选用哪一种模型就要考虑到节约和计算量，选定模型后，剩下的任务就是用适当的算法估计模型参数$(a_k,b_k,p,q)$，以便用模型对随机信号进行预测。

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

对于$x(n)=w(n)-{\sum }_{k=1}^p{a}_{k}x(n-k)$，两边同时乘以$x(n-m)$然后求均值：

$$E[x(n)x(n-m)]=E\left[w(n)x(n-m)-\sum _{k=1}^p{a}_{k}x(n-k)x(n-m)\right]$$
由自相关函数的对称性，可以将原式化为：
$$R_{xx}(m)=R_{xw}(m)-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)$$
而系统的单位冲激相应$h(n)$是因果的，所以输出的平稳随机信号和输入的白噪声之间的互相关函数有下列推导：

$R_{xw}(m)$的定义为
$$R_{xw}(m)=E[x(n)w(n+m)]$$

而随机信号的当前值又是随机信号的当前激励通过滤波器产生的，所以有：
$$x(n)=w(n)\ast h(n)=\sum _{k=0}^{\infty }h(k)w(n-k)$$
由于均值是对随机变量来求的，所以均值时可以把$h(n)$当作常量。
$$\begin{gathered} R_{xw}(m)=E\left[\sum _{k=0}^{\infty }h(k)w(n-k)w(n+m)\right]=\sum _{k=0}^{\infty }h(k)E[w(n-k)w(n+m)]] \\ =\sum _{k=0}^{\infty }h(k)R_{ww}(m+k)=\sum _{k=0}^{\infty }h(k){\sigma }_w^2\delta (m+k)={\sigma }_w^{2}h(-m) \end{gathered}$$
所以有：
$$R_{xw}(m)=\{\begin{array}{cc}0& m>0\\ {\sigma }_w^{2}h(-m)& m\le 0\end{array}$$
综上，有：

$$R_{xx}(m)=\{\begin{array}{cc}-{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)& m>0\\ -{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)+h(0){\sigma }_w^2& m=0\\ R_{xx}(-m)& m<0\end{array}$$

由$$H(z)=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$转化为时域可以得到
$$h(n)+\sum _{k=1}^p{a}_{k}h(n-k)=\delta (n)$$
所以有$h(0)=\delta (0)=1$
根据$R_{xx}$的表的式可以看出，AR模型输出信号的自相关函数具有递推的性质：
$$R_{xx}(m)=-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)m>0$$
这就是著名的Yule-Walker方程，将上式变换：
$$0=R_{xx}(m)+\sum _{k=1}^p{a}_k{R}_{xx}(m-k)m>0$$
可以求得输入的白噪声方差为：
$${\sigma }_w^2=R_{xx}(0)+\sum _{k=1}^p{a}_k{R}_{xx}(-k)m=0$$
最后，结果为：
$$\left[\begin{array}{cccc}R(0)& R(-1)& \cdots & R(-p)\\ R(1)& R(0)& \cdots & R(-p+1)\\ ⋮& ⋮& \cdots & ⋮\\ R(p)& R(p-1)& \cdots & R(0)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}{\sigma }_w^2\\ 0\\ ⋮\\ 0\end{array}\right]$$

方程组的系数都是自相关矩阵$[R{]}_{p+1}$，由于自相关函数是偶对称函数：$R_{xx}(m)=R_{xx}(-m)$，因而自相关矩阵是对称矩阵，与主对角线平行的斜对角线上的元素都是相同的，是的维托毕兹（Toeplitz）矩阵，所以存在高效算法，其中应用广泛的有Levinson－Durbin（L-D）算法。Yule－Walker（Y-W）方程表明，只要已知输出平稳随机信号的自相关函数，就能求出AR模型中的参数$a_k$，并且需要的观测数据较少。

实例

自回归信号模型$AR(3)$为：
$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)$，式中$w(n)$式具有方差${\sigma }_w^2=1$的平稳白噪声，求：

1. 自相关序列$R_{xx}(m)$,m=1,2,3,4,5,6,72. 用求出来的自相关序列估计$AR(3)$的参数$a_k$以及输入白噪声方差$\widehat{{\sigma }_w}$的大小。
2. 用a求出的自相关序列来估计AR（3）的参数${\hat{a}}_k$，以及输入白噪声的方差$\widehat{{\sigma }_w}$的大小
3. 利用AR模型，仿真给出的32点观测值x(n)＝[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],用观测值的自相关序列来估计$AR(3)$的参数$\widehat{{\sigma }_w^2}$,以及输入白噪声的$\widehat{{\sigma }_w}$。

求解

1. 题干给出了该信号本身的表达式，可以得到$a_1=-\frac{14}{24},a_2=-\frac{9}{24},a_3=\frac{1}{24}$那么可以直接使用公式：
   $$\left[\begin{array}{cccc}R(0)& R(1)& R(2)& R(3)\\ R(1)& R(0)& R(1)& R(2)\\ R(2)& R(1)& R(0)& R(1)\\ R(3)& R(2)& R(1)& R(0)\end{array}\right]\left[\begin{array}{c}1\\ -14/24\\ -9/24\\ 1/24\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]$$
   利用matlab解方程：

   	a = [-14/24,-9/24,1/24 ];
   	
   	A = [ 1, a(1), a(2), a(3);  
   	    a(1), 1 + a(2), a(3), 0;
   	    a(2), a(1)+a(3), 1, 0; 
   	    a(3), a(2), a(1), 1]
   	B = [1;0;0;0];
   	Rxx = A \ B
   

   可以得到

   Rxx =
   4.9377
   4.3287
   4.1964
   3.8654
   

   此后可以使用公式$$R_{xx}(m)=-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)m>0$$来计算R(4), R(5)等

   for i = 5:7
       Rxx(i) = 0;
       for k = 1:3
           Rxx(i) = Rxx(i) - a(k) * Rxx(i-k);
       end
   end
   Rxx
   

   结果如下：

   Rxx =
   
   4.9377
   4.3287
   4.1964
   3.8654
   3.6481
   3.4027
   3.1919
   

2. 利用$$\left[\begin{array}{llll}R(0)& R(1)& R(2)& R(3)\\ R(1)& R(0)& R(1)& R(2)\\ R(2)& R(1)& R(0)& R(1)\\ R(3)& R(2)& R(1)& R(0)\end{array}\right]\left[\begin{array}{c}1\\ {\hat{a}}_1\\ {\hat{a}}_2\\ {\hat{a}}_3\end{array}\right]=\left[\begin{array}{c}{\hat{\sigma }}_w^2\\ 0\\ 0\\ 0\end{array}\right]$$
   如同上一问，直接解方程即可。解得${\hat{a}}_1=-\frac{14}{24},{\hat{a}}_2=-\frac{9}{24},{\hat{a}}_3=\frac{1}{24},{\hat{\sigma }}_w^2=1$
   可以发现对AR模型参数是无失真的估计，因为已知AR模型，我们可以得到完全的输出观测值，因而求得的自相关函数没有失真，当然也就可以不失真的估计。
3. 利用公式$$R_{xx}(m)=\frac{1}{n}\sum _{i=1}^n{x}_i{x}_{i+m}m>0$$可以求得自相关系数。

   x_n = [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];
   
   R_xx = zeros( 1, length(x_n) );
   
   for m  = 0 : length(x_n) - 1 
       j = 1; 
       while j + m <= length(x_n)
           R_xx(m+1) = R_xx(m+1) + x_n(j) * x_n(j + m );
           j = j + 1; 
       end
       R_xx(m+1) =  R_xx(m+1) / 32;
   end
   R_xx
   

   结果如下：

   R_xx =
   
     1 至 11 列
   
       1.9271    1.6618    1.5381    1.3545    1.1349    0.9060    0.8673    0.7520    0.7637    0.8058    0.8497
   
     12 至 22 列
   
       0.8761    0.9608    0.8859    0.7868    0.7445    0.6830    0.5808    0.5622    0.5134    0.4301    0.3998
   
     23 至 32 列
   
       0.3050    0.2550    0.1997    0.1282    0.0637    0.0329   -0.0015   -0.0089   -0.0143   -0.0083
   

   带入矩阵$$\left[\begin{array}{llll}R(0)& R(1)& R(2)& R(3)\\ R(1)& R(0)& R(1)& R(2)\\ R(2)& R(1)& R(0)& R(1)\\ R(3)& R(2)& R(1)& R(0)\end{array}\right]\left[\begin{array}{c}1\\ {\hat{a}}_1\\ {\hat{a}}_2\\ {\hat{a}}_3\end{array}\right]=\left[\begin{array}{c}{\hat{\sigma }}_w^2\\ 0\\ 0\\ 0\end{array}\right]$$可以解得
   ${\hat{a}}_1=－0.6984,{\hat{a}}_2=＝－0.2748,{\hat{a}}_3=0.0915,{\hat{\sigma }}_w^2=.04678$
   与真实AR模型参数误差为：$e_1＝0.1151，e_2＝0.1002，e_3＝0.0498$，原因在于我们只有一部分的观测数据，使得自相关序列值与理想的完全不同。输入信号的方差误差比较大：$e_{\sigma }＝0.5322$，造成的原因比较多，计算机仿真的白噪声由于只有32点长，32点序列的方差不可能刚好等于1。给出一段观测值求AR模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

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

在上述实例中，要得到更精确的估计值，就要建立更高阶的AR模型，直接用观测值的自相关序列来求解Y-W方程计算量太大。因此把AR模型和预测系统联系起来，换个方法来估计AR参数。
从AR模型的时域表达式$$x(n)=w(n)-\sum _{k=1}^p{a}_{k}x(n-k)$$可以知道模型的当前输出值与它过去的输出值有关。预测是推断一个给定序列的未来值，即利用信号前后的相关性来估计未来的信号值。

若序列的模型已知而用过去观测的数据来推求现在和将来的数据称为前向预测器，表示为：$$\hat{x}(n)=-\sum _{k=1}^m{a}_m(k)x(n-k)$$
式中$a(k),k=1,2,\cdots ,m$代表m阶预测器的预测系数，负号是为了与技术文献保持一致。显然预测出来的结果与真实的结果存在预测误差或前向预测误差，设误差为$e(n)$:
$$e(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^m{a}_m(k)x(n-k)$$
把$e(n)$看成系统的输出，$x(n)$看成系统的输入，得到系统函数：
$$\frac{E(z)}{X(z)}=1+\sum _{k=1}^m{a}_m(k)z^{-k}$$

假如m=p，且预测系数和AR参数模型参数相同，(如上图所示，左为预测误差系统，右为AR模型)那么就有$w(n)=e(n)$，即前向预测误差系统中的输入为$x(n)$，输出为预测误差$e(n)$等于白噪声。也就是说前向预测误差系统对观测信号起了白化的作用。由于AR模型和前向预测误差系统有着密切的关系，两者的系统函数互为倒数，所以求AR模型参数就可以通过求预测误差系统的预测系数来实现。
对$$e(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^m{a}_m(k)x(n-k)$$
中的预测误差求均方值，可得
$$\begin{gathered} E[e^2(n)]=E[x(n)+\sum _{k=1}^m{a}_m(k)R_{xx}(k)+\sum _{k=1}^m\sum _{l=1}^m{a}_m(l)a_m(k)R_{xx}(l-k)] \\ =R_{xx}(0)+2\left[\sum _{k=1}^m{a}_m(k)R_{xx}(k)\right]+\sum _{k=1}^m\sum _{l=1}^m{a}_m(l)a_m(k)R_{xx}(l-k) \end{gathered}$$
要使得均方误差最小，将上式右边对预测系数求偏导并且等于零，得到m个等式：
$$R_{xx}(l)=-\sum _{k=1}^m{a}_m(k)R_{xx}(l-k)l=1,2,\cdots m$$
带入均方误差公式$R_{xx}(0)+2\left[{\sum }_{k=1}^m{a}_m(k)R_{xx}(k)\right]+{\sum }_{k=1}^m{\sum }_{l=1}^m{a}_m(l)a_m(k)R_{xx}(l-k)$，可以得到最小均方误差：
$$E_m\left[e^2(n)\right]=R_{xx}(0)+\sum _{k=1}^m{a}_m(k)R_{xx}(k)$$
或
$$\begin{array}{c}{E}_p\left[e^2(n)\right]=R_{xx}(0)+{\sum }_{k=1}^{T-}{a}_k{R}_{xx}(k)\\ a_k=a_m(k)m=p\end{array}$$
也就是p阶预测器的预测系数等于p阶AR模型的参数，由于$e(n)=w(n)$，所以最小均方预测误差等于白噪声方差，即$E_p[e^2(n)]={\sigma }_m^2$

有了上面的知识后，我们回来看怎样估计AR模型参数，也即要估计参数$a_k,p,{\sigma }_w^2$，这里介绍应用广泛的L-D算法。L-D算法的基本思想就是根据Y-W方程式或上面介绍的两个式子，自相关序列具有递推的性质，L-D递推算法是模型阶数逐渐加大的一种算法，先计算阶次m＝1时的预测系数 { a m ( k ) } ∣ = a 1 ( 1 ) \{a_m(k)\}|=a_1(1) {am​(k)}∣=a1​(1)和${\sigma }_{w1}^2$然后计算m＝2时的 { a m ( k ) } ∣ = a 2 ( 1 ) \{a_m(k)\}|=a_2(1) {am​(k)}∣=a2​(1),$a_2(2)$和${\sigma }_{w1}^2$一直计算到m=p阶时，$a_p(1),a_p(2),\cdots ,a_p(p)$以及${\sigma }_{wp}^2$。这种递推算法的特点是，每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。
最终可以得到：
$$\{\begin{array}{c}{a}_m(k)=a_{m-1}(k)+a_m(m)a_{m-1}(m-k)\\ a_m(m)=-\frac{R(m)+{\sum }_{k=1}^{m-1}{a}_{m-1}(k)R(m-k)}{E_{m-1}}\\ E_m={\sigma }_{wm}^2=\left[1-a_m^2(m)\right]E_{m-1}=R(0){\prod }_{k=1}^m\left[1-a_k^2(k)\right]\end{array}$$
其中$a_m(m)$称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值$E_0=R(0),a_0(0)=1$以及AR模型的阶数，就可以按照图下图所示流程图进行估计啦。

L-D算法的优缺点

• 优点
  • 计算速度快
  • 求得的AR模型必定稳定
  • 均方预测误差随着阶次的增加而减小
• 缺点
  • 由于在求自相关序列时，是假设除了观测值之外的数据都为零，必然会引入较大误差。

实例

自回归信号模型$AR(3)$为：
$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)$，式中$w(n)$式具有方差${\sigma }_w^2=1$的平稳白噪声.
利用AR模型，仿真给出的32点观测值x(n)＝[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],用观测值的自相关序列来估计$AR(3)$的参数$\widehat{{\sigma }_w^2}$,以及输入白噪声的$\widehat{{\sigma }_w}$。

对于该问题，采用L-D算法进行运算。

首先使用相同的做法算出$R_{xx}(m)$。之后利用L-D算法。
由于matlab下标从1开始，所以对于下标需要进行一些特殊处理。

x_n = [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];

R_xx = zeros( 1, length(x_n) );

for m  = 0 : length(x_n) - 1 
    j = 1; 
    while j + m <= length(x_n)
        R_xx(m+1) = R_xx(m+1) + x_n(j) * x_n(j + m );
        j = j + 1;
    end
    R_xx(m+1) =  R_xx(m+1) / 32;
end

%初始化
p = 4;
E = zeros(1, p); E(1) = R_xx(1);
a = zeros(p, 3); a(1 , 1) = 1;

for m = 1:3
    a( m + 1, m ) = R_xx(m + 1);
    for k = 1:( m - 1 )
        a( m + 1, m ) = a( m + 1, m ) + a( m, k ) * R_xx( m + 1 - k );
    end
    a( m + 1, m ) = - a( m + 1, m ) / E(m);
    for k = 1:(m - 1)
       a( m + 1, k ) = a( m, k ) + a( m + 1, m ) * a( m, m - k );
    end 
    E( m + 1 ) = ( 1 - a( m + 1, m ).^2 ) * E(m); 
end
a
E

最终结果为：

a =

    1.0000         0         0
   -0.8623         0         0
   -0.6789   -0.2127         0
   -0.6984   -0.2748    0.0915


E =

    1.9271    0.4941    0.4718    0.4678

此外，Matlab有专门的函数实现L-D算法的AR模型参数估计：[a, E]＝aryule（x，p），输入x表示观测信号，输入p表示要求的阶数，输出a表示估计的模型参数，e表示噪声信号的方差估计。例如本题用该函数计算：

[a, E] = aryule(x_n,3);
a
E

结果为：

a =

    1.0000   -0.6984   -0.2748    0.0915


E =

    0.4678

可以看到结果是一样的。
如果使用更高阶的AR模型，可以看到均方误差会越来越小。
给定观测序列后我们很容易用L-D算法来进行AR模型参数的估计，但是如何确定阶数以及L-D算法存在的缺点如何克服呢？

AR模型参数估计的各种算法的比较和阶数的选择

为了克服L-D算法导致的误差，1968年Burg提出了Burg算法，其基本思想是对观测的数据进行前向和后向预测，然后让两者的均方误差之和为最小作为估计的准则估计处反射系数，进而通过L-D算法的递推公式求出AR模型参数。Burg算法的优点是，求得的AR模型是稳定的，较高的计算效率，但递推还是用的L-D算法，因此仍然存在明显的缺点。
Matlab中有专门的函数实现Burg算法的AR模型参数估计：[a, E]＝arburg（x，p）。例如上面的例子用Burg算法计算结果为：

[a, E] = arburg(x_n,p);
a
E

结果：

a =

    1.0000   -0.6852   -0.3088   -0.0489    0.1759


E =

    0.4426

结果与L-D算法结果略有不同，${\hat{a}}_3$估计得更精确些。

1980年Marple在前人的基础上提出一种高效算法，Marple算法或者称不受约束的最小二乘法（LS），该算法的思想是，让每一个预测系数的确定直接与前向、后向预测的总的平方误差最小，这样预测系数就不能由低一阶的系数递推确定了，所以不能用L-D算法求解，实践表明该算法比L-D、Burg算法优越。由于该算法是从整体上选择所有的模型参数达到总的均方误差最小，与自适应算法类似，不足是该算法不能保证AR模型的稳定性。
AR模型的阶数选择不同得到的模型不同，效果相差较大，因而如何选择阶数很重要。因此国内外学者在这方面都做了许多研究工作，其中基于均方误差最小的最终预测误差（FPE：final predidyion error）准则是确定AR模型阶次比较有效的准则。
最终预测误差（FPE：final predidyion error）准则定义：给定观测长度为N，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差，该误差在某个阶数时为最小，其表达式为：
$$FPE(p)={\hat{\sigma }}_{wp}^2\left(\frac{N+P+1}{N-P-1}\right)$$
估计的方差随着阶数的增加而减小，而括号内的值随着p的增加而增加，因而能找到一最佳的p，使得FPE最小。
除了用上式进行阶数的选择，还可以采用试验的方法，如果某一个p值满足预先规定的某些要求，则认为该p值是最优的。在短数据情况下，根据经验，AR模型的阶次选在N/3～N/2的范围内较好。