简介：围绕Yule-Walker方程求解与AR模型建立，这份实验报告PDF系统整理了生物医学信号处理中的关键方法。内容从随机信号的自回归模型出发，讲解Yule-Walker方程的推导、自相关矩阵构造及L-D快速算法，并给出完整的Matlab实现流程。实验部分以心电、脑电等实际生理信号为对象，完成AR建模、系数求解、白噪声驱动仿真和功率谱对比，同时通过最小均方误差、预测误差及FPE指标评估模型精度，可与Matlab内置aryule函数结果互相验证。报告包含实验目的、原理、步骤、结果图表与程序代码，单文件PDF约847KB。目前已有182人学习下载，适合生物医学工程、信号处理相关专业学生用于课程实验、复习或自学参考。

1. YuleWalker方程.pdf：AR模型参数估计为什么绕不开这张纸

点开“YuleWalker方程.pdf”的人，大多数不是来欣赏推导过程的，而是手里攒了一列时序数据——风速、脑电、设备振动、量化收益——想用AR模型估一组系数，结果发现最短路不是直接调statsmodels，而是先弄明白这一组方程在干什么。Yule-Walker方程干的事用一个公式就能说清：把“自相关序列等于AR系数与过去自相关的线性组合”写成矩阵方程，解方程就得到AR系数。这一步是线性预测、LPC语音编码、AR功率谱估计的地基。适合谁？做时间序列预测、信号特征提取、用Python科学计算栈的从业者。看完这篇，你能自己写出求解函数，并知道什么时候该信它、什么时候该换Burg或最小二乘。

2. 从自相关到Toeplitz矩阵：Yule-Walker方程的推导与结构

2.1 AR(p)模型下自相关满足的约束关系

假设零均值平稳序列x_t满足AR(p)模型：

x_t = φ_1 x_{t-1} + φ_2 x_{t-2} + ... + φ_p x_{t-p} + ε_t

对等式两边同时乘以x_{t-k}再取期望。因为k >= 1时ε_t与过去的x不相关，最后一项直接消失，剩下：

γ_k = φ_1 γ_{k-1} + φ_2 γ_{k-2} + ... + φ_p γ_{k-p}, k = 1, 2, ..., p

这里γ_k就是滞后k步的自协方差。把k=1到p逐个写出来，就是Yule-Walker方程。注意两个前提：序列零均值、协方差平稳。如果序列带趋势，比如股价原始收盘价，γ_k随时刻变化，这个等式从第一步就不成立——这是后文会单独讲的坑。

实际工程里常用样本自相关系数替代理论自相关。计算时有一个选择：分母除以N还是N-k。这两个版本后面结果差异很大，我会在第4章专门讲，这里先按教科书惯例用除以N的“有偏估计”，因为它在短序列下能保住Toeplitz矩阵的正定性，递推不容易发疯。

2.2 方程组的矩阵形式：Toeplitz结构为什么值得专门讲

把上面的p个等式写成矩阵：

R φ = ρ

其中R是一个p×p矩阵，第i行第j列等于γ_{|i-j|}，右边的ρ是[γ_1, γ_2, ..., γ_p]的转置。把自相关写开你马上能看到R每条对角线的值都一样：

[[γ_0, γ_1, γ_2, ..., γ_{p-1}], [γ_1, γ_0, γ_1, ..., γ_{p-2}], [γ_2, γ_1, γ_0, ..., γ_{p-3}], ...

这条对角线相等的性质叫Toeplitz。第一次接触的人会觉得它就是“矩阵长得整齐一点”，实际价值在于求解速度。通用高斯消元法解这个方程组是O(p³)，而利用Toeplitz结构做Levinson-Durbin递推只需要O(p²)次乘加，p到50以上差距就非常明显。在线估计场景里，每次新数据进来都要解一次方程，这个复杂度差距基本决定了能不能实时跑。

我通常会在解方程前先检查R是否满足Toeplitz，方法很简单：对比R[i][j]和R[i+1][j+1]是否相等。由于浮点误差和样本自相关的计算顺序，这里偶尔会差出10⁻⁶量级，不影响求解，但如果你用无偏自相关而且滞后段取得很大，对角线差异可能会被放大到足以让矩阵失去正定性。

2.3 两种求解路径：直接解与Levinson-Durbin递推

第一种路径是把矩阵R显式构造出来，用NumPy的np.linalg.solve直接解。好处是代码直观、不容易写错；坏处是当p超过100时内存和耗时都开始吃紧，而且你手里其实有更强的工具——scipy.linalg.solve_toeplitz，后面会写。

第二种路径是Levinson-Durbin递推，这是工程上真正在用的算法。它从一阶开始，逐阶把模型从AR(1)升级到AR(p)，每一阶只用到前一阶的系数和一个反射系数k_m。反射系数在信号处理里也叫PARCOR系数，它天然落在(-1, 1)区间内，这是AR模型稳定的充要条件。递推公式长这样：

k_m = (γ_m - Σ_{j=1}^{m-1} a_{m-1,j} γ_{m-j}) / E_{m-1} a_{m,m} = k_m a_{m,j} = a_{m-1,j} - k_m * a_{m-1,m-j} (j = 1..m-1) E_m = E_{m-1} * (1 - k_m²)

其中E_m是m阶模型的前向预测误差方差，也就是残差方差。这套递推每算一阶还顺带给了你一个非常有用的副产品：k_m本身就是滞m处的偏自相关函数值。于是选阶数时不用再单独调包计算PACF，递推过程里已经全有了。

3. 用NumPy和SciPy求解Yule-Walker方程：三步拿到AR系数

3.1 造一段AR(2)仿真数据作为测试集

手边没有合适的实测数据时，先生成一段已知系数的AR(2)序列最靠谱，这样能直接对比解出来的系数和真实值的偏差。下面这段是按φ=[0.6, -0.4]生成的：

import numpy as np from scipy.linalg import solve_toeplitz import matplotlib.pyplot as plt np.random.seed(42) N = 500 phi_true = np.array([0.6, -0.4]) x = np.zeros(N) eps = np.random.randn(N) * 0.5 for t in range(2, N): x[t] = phi_true[0] * x[t-1] + phi_true[1] * x[t-2] + eps[t]

逻辑说明：自回归生成必须循环，不能直接用np.convolve套白噪声，因为卷积假设系统初始条件为零且输入无限长，边界效应会污染前几十个点。phi_true里索引0对应滞后1步，索引1对应滞后2步。eps方差取0.5是让信噪比低一点，检验求解在噪声较大时是否还稳。这个仿真序列后面所有代码都用同一份。

参数说明：N=500足够让前几个滞后的样本自相关收敛到理论值；如果你想刻意演示短序列的不稳定性，可以把N改成25再对比一次，这正是第4章第一个坑的实验环境。种子固定是为了结果可复现。

3.2 计算样本自相关：优先选有偏估计

这里不要直接用np.correlate，它的返回长度和模式容易把人绕晕。一个明确的循环更不容易错：

def autocorr_biased(x, max_lag): N = len(x) x = x - np.mean(x) r = np.zeros(max_lag + 1) for k in range(max_lag + 1): r[k] = np.dot(x[:N-k], x[k:]) / N return r p = 3 r = autocorr_biased(x, p) print("r0 =", r[0], " r1 =", r[1], " r2 =", r[2], " r3 =", r[3])

逻辑说明：先减去均值满足零均值假设。np.dot(x[:N-k], x[k:])计算的是滞后k步的乘积和，除以N（不是N-k），得到的就是有偏自协方差估计。有偏估计的方差比无偏估计小，而且能保证后面构造出的Toeplitz矩阵正定，这两点对Y-W求解比“无偏”这个优点重要得多。

参数说明：max_lag至少要等于模型阶数p。如果你后面要用Levinson-Durbin递推往上算到p，这里多留几个滞后没坏处，代价只是多几个点乘。对500个点的序列，求到20阶滞后也没有性能压力。

3.3 用SciPy的solve_toeplitz解矩阵方程

SciPy里linalg.solve_toeplitz专门解Toeplitz矩阵方程，比np.linalg.solve快一个量级。它的参数是(第一列, 第一行)和右侧向量。因为自相关对称，这里第一列和第一行相等：

def yule_walker_solve(x, p): N = len(x) r = autocorr_biased(x, p) r0 = r[0] r_n = r / r0 # solve_toeplitz要求c[0] == r[0]，归一化后必然满足 phi = solve_toeplitz((r_n[:p], r_n[:p]), r_n[1:p+1]) sigma2 = r0 * (1 - np.dot(phi, r_n[1:p+1])) return phi, sigma2 phi, sigma2 = yule_walker_solve(x, 2) print("phi =", phi, " sigma2 =", sigma2)

逻辑说明：solve_toeplitz第一个元组里，c代表矩阵的第一列，r代表第一行。因为R[i][j] = γ_{|i-j|}，第一列是[γ_0, γ_1, ..., γ_{p-1}]，第一行完全一样，所以传两遍r_n[:p]。函数要求c[0]必须等于r[0]，归一化后两个都是1，这就是先除r0的原因。右侧传r_n[1:p+1]，对应上面推到过的[γ_1, ..., γ_p]。

参数说明：p=2时理论上解就是[0.6, -0.4]附近；p=3时第三个系数应该接近0，代表过拟合不吸收额外系数。残差方差sigma2用γ_0 - Σφ_i γ_i这个公式，它等价于Levinson递推里的E_p，只是少了一次递推的积累误差，实测两者在小数点后两位数内一致。如果你要用AR系数做功率谱估计，sigma2记得乘回r0，因为r_n已经是归一化的自相关，直接拿归一化结果算谱会把能量尺度弄丢。

3.4 手动实现Levinson-Durbin递推

上面那个算法内部其实已经用了Levinson结构，但自己写一遍递推能帮你彻底弄懂第2.3节的公式，调试时也看得见每一步的中间量：

def levinson_durbin(r, p): a = np.zeros(p + 1) a[0] = 1.0 E = float(r[0]) for m in range(1, p + 1): s = float(r[m]) for j in range(1, m): s -= a[j] * r[m - j] k = s / E a_prev = a.copy() for j in range(1, m): a[j] = a_prev[j] - k * a_prev[m - j] a[m] = k E *= (1.0 - k * k) return a[1:], E r_full = autocorr_biased(x, 2) phi_ld, E_ld = levinson_durbin(r_full, 2) print("Levinson phi =", phi_ld, " E =", E_ld)

逻辑说明：内层第一个循环计算γ_m - Σ a_j γ_{m-j}，这正是反射系数k_m的分子。第二个循环做系数反向更新，核心是两行之间的对称关系：AR系数在阶数从m-1升到m时，旧系数会被k_m与反序旧系数的组合修正。E每阶乘以(1-k²)，只要|k| < 1，预测误差就逐阶单调下降，这是判断递推是否数值健康的天然指标。

参数说明：r的第0个元素r[0]不能为0，全零序列会在这里直接报错。a数组长度比阶数多1，a[0]固定为1.0不参与最终输出，返回时从a[1:]取值。对比3.3节的solve_toeplitz结果，两个phi在浮点精度上应当一致；如果出现明显偏差，大概率是你手写递推里索引写错，可以在m=1时先打印k对照r[1]/r[0]，m=2时再打印一次，逐阶定位。

3.5 残差方差与BIC选阶：把方程结果用起来

Y-W方程本身不负责回答“p选几”，需要外在准则。常见做法是遍历p=1..max_p，计算每个阶数下的BIC：

max_p = 10 bic_list = [] for p in range(1, max_p + 1): r = autocorr_biased(x, p) phi_p, E_p = levinson_durbin(r, p) bic = N * np.log(E_p) + p * np.log(N) bic_list.append(bic) best_p = int(np.argmin(bic_list)) + 1 print("best_p =", best_p)

逻辑说明：BIC第一项是残差方差的负对数似然，阶数增加时它会下降；第二项p * log(N)是惩罚项，防止系数越多越好。Y-W解出的E_p拿来直接算BIC非常省事，不用重新拟合。注意BIC对样本量敏感，N只有30时惩罚项几乎失效，这时更建议用AICc或交叉验证，但N在几百以上BIC表现更稳。

参数说明：max_p别拍脑袋取大。经验上限是min(N/10, 50)，对500点数据取10已经含了余量。阶数超过这个上限后，样本自相关的尾部噪声会开始支配E_p，BIC曲线会出现十几阶后继续单调下行的假象。

4. Yule-Walker方程求解避坑：五个真实翻车点

4.1 短序列下系数被夸大：同一个AR信号，N=20和N=500解出两套结果

现象：用同一段AR(2)仿真，把N改成20，解出来的系数可能变成[0.82, -0.65]，而真实值是[0.6, -0.4]。继续缩短到N=15，系数甚至会冒出模大于1的情况。原因：Y-W方程用的是样本自相关替代理论自相关，短序列时滞后2步、3步的自相关估计方差大得惊人，而高阶滞后在方程里的权重又不低，一点误差就被放大进系数。解决：N小于30时不要单独信Y-W，优先用Burg法或其他最小二乘类估计；如果一定要用Y-W，把阶数上限压到3以内，并用多段子序列交叉验证系数稳定区间。

4.2 solve_toeplitz报错或者解出对不上的值：第一列和第一行传反了

现象：调用solve_toeplitz((r[:-1], r[1:]), r[1:])时，要么长度对不上直接抛异常，要么解出来的系数和np.linalg.solve结果差很多。原因：solve_toeplitz的参数约定是(c, r)分别代表Toeplitz矩阵的第一列和第一行，改成其他形式时函数内部按T[i][j]的列行索引去取，结果自然错位。解决：自相关对称且第一列第一行完全一致，永远传(r[:p], r[:p])，右侧固定传r[1:p+1]；调用前加一行归一化r_n = r / r[0]，避免c[0] != r[0]的运行时检查报错。

4.3 高阶系数震荡：BIC还在降，预测误差先炸了

现象：p从5往上加，BIC一路变小，看起来模型越来越好；但看系数时发现φ_5=0.44、φ_6=-0.38，符号交替，用这套参数做一步预测的均方误差反而比p=3时大了30%。原因：阶数升高后，Y-W方程右端的自相关向量开始进入尾部噪声区，这些噪声被当成真实周期成分吸收进系数，模型在训练段过拟合了。解决：看BIC的同时加一个条件——最大阶数不得超过min(N/10, 50)；再交叉验证一次，把预测误差作为最终判据，BIC只用于初筛。遇到系数符号震荡，直接砍半阶数。

4.4 非平稳序列直接解：φ_1被估成0.99

现象：拿一段带趋势的原始序列（比如逐日累计值）直接跑Y-W，p取5，解出φ_1≈0.99，其余系数都很小，残差方差并没有随阶数下降。原因：非平稳序列的自相关不衰减，Y-W方程从推导阶段就不成立。φ_1≈1只是方程强行拟合出的结果，不是数据里有长记忆结构。解决：先做ADF检验或直接看自相关图，r_k衰减到0很慢就说明要先差分。差分后重新求自相关，通常一阶差分后就能看到快速衰减，再跑Y-W就正常了。

4.5 Levinson递推里反射系数大于1：无偏自相关在捣乱

现象：手写Levinson-Durbin递推时，算到m=3发现k>1，E变成负值，程序输出一堆nan。原因：递推里k的分母是上一阶的误差方差E_{m-1}，它必须是正数；当样本自相关用的是除以N-k的无偏估计时，尾部滞后对应的方差可能把Toeplitz矩阵推得不正定，E被算成负数。解决：切回有偏估计，r[k] = Σx_t x_{t-k} / N，不要除以N-k。这个选择不是“精度”而是“稳定性”的权衡——有偏估计虽然在小滞后上有微小偏差，但能保证递推不翻车，实践中利远大于弊。

5. 用Yule-Walker结果验证模型：残差检验与AR谱估计

跑完Y-W方程拿到系数只是第一步，验证得到的模型真的白化数据才是关键。我习惯依次做三件事。

先看残差是否还是白噪声。残差resid = x_t - Σφ_i x_{t-i}，对一段500点的数据，用Ljung-Box Q统计量检验前10阶自相关：

resid = np.zeros(N) for t in range(2, N): resid[t] = x[t] - phi[0]*x[t-1] - phi[1]*x[t-2] h = 10 acf = np.array([np.dot(resid[:N-k], resid[k:]) / np.dot(resid, resid) for k in range(1, h+1)]) Q = N * (N+2) * np.sum(acf**2 / (N - np.arange(1, h+1))) df = h - 2 # 减去AR阶数 from scipy.stats import chi2 print("Q =", Q, " p =", 1 - chi2.cdf(Q, df))

p大于0.05说明残差里没有显著剩余自相关，模型可以接受。注意自由度要减去p，不然Q检验偏高，容易拒绝本该接受的模型。

再看偏自相关函数。Levinson递推里每一阶的k_m其实就是PACF在滞后m处的值，所以直接复用3.4节的递推结果，找到k_m落在两倍标准误带之外的最后一个位置，那个位置就是推荐阶数。这里有一个常见陷阱：Y-W解出的PACF在N较小时尾部会有少量假显著点，别看到第7阶超过带宽就激动，优先落在低阶位置才可信。

最后用Y-W系数做AR谱估计，这是线性预测之外最有用的落地方式。AR(p)谱密度公式为S(f) = σ² / |1 - Σ_{i=1}^p φ_i e^{-j2πfi}|²，频点一块儿算：

f = np.linspace(0, 0.5, 400) omega = 2 * np.pi * f den = np.ones(len(f), dtype=complex) for i in range(len(phi)): den -= phi[i] * np.exp(-1j * omega * (i + 1)) S = sigma2 / np.abs(den)**2 plt.plot(f, S) plt.xlabel("frequency") plt.show()

AR谱估计比直接做FFT平滑，短数据下更稳，特别适合EEG频带分析和机械振动特征提取。公式里的σ²要用带尺度的残差方差，也就是3.3节里没有归一化的sigma2。

现在拿到一段新序列，我的流程固定了：先画自相关图判断平稳性，差分到自相关快速衰减；再用Levinson递推一路算反射系数和BIC，确定阶数；最后用Ljung-Box确认残差白噪声。整套下来Y-W方程只是第一站，但这一站做扎实了，后面无论接谱估计还是预测都省心很多。希望帮到你。

本文还有配套的精品资源，点击获取