文章目录

系统模型
似然函数
辨识过程
关于参数的梯度与海森矩阵
- 梯度
- 海森矩阵

受控自回归滑动平均模型 (Controled Auto Regression and Moving Average model, CARMA)，亦称带外部输入的自回归滑动平均模型 (Auto Regression and Moving Average model with eXogenous input, ARMAX)是应用非常广泛的线性系统模型，本文介绍该模型的一种系统辨识方法：最大似然法。

系统模型

yk+a1yk−1+⋯anyk−n=b1uk−1+bnuk−n+ek+c1ek−1+⋯+cnek−ny_k + a_1y_{k-1} + \cdots a_n y_{k-n} = \\ b_1u_{k-1}+ b_n u_{k-n} + e_{k} + c_1 e_{k-1} + \cdots + c_n e_{k-n} yk+a1yk−1+⋯anyk−n=b1uk−1+bnuk−n+ek+c1ek−1+⋯+cnek−n
即
yk=−∑i=1naiyk−i+∑i=1nbiuk−i+∑i=1nciek−i+eky_k = -\sum_{i=1}^n a_iy_{k-i} + \sum_{i=1}^n b_i u_{k-i}+ \sum_{i=1}^n c_i e_{k-i} + e_{k} yk=−i=1∑naiyk−i+i=1∑nbiuk−i+i=1∑nciek−i+ek
该模型在 ARMA 的基础上考虑了外部输入 u(k)u(k)u(k) 对输出的影响，为方便讨论，在此设定 a,b,ca,b,ca,b,c 的下标都是从 111 到 nnn，实际上不必如此。

似然函数

在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。
给定输出x时，关于参数θ的似然函数L(θ|x)（在数值上）等于给定参数θ后变量X的概率：

L(θ∣X)=P(X∣θ)L(\theta | X) = P(X|\theta)L(θ∣X)=P(X∣θ)

在统计学中，“似然性”和“概率”有明确的区分。概率用于在已知一些参数的情况下，预测接下来的观测所得到的结果，而似然性则是用于在已知某些观测所得到的结果时，对有关事物的性质的参数进行估计。

假定误差 e(k)e(k)e(k) 服从 000 均值、方差为 σ2\sigma^2σ2 的高斯分布，则有似然函数为：
P(YN∣UN,Θ)=P(yN,…∣uN,…,Θ)=P(y0)Πk=1NP(yk∣Yk−1,UN,Θ)=P(y0)Πk=1NP(ek)=P(y0)(2π)−N/2σ−Nexp⁡[−(1/2σ2)∑k=1Nek2](1)\begin{array}{ll} P(Y_N | U_N,\Theta) &= P(y_N, \ldots | u_N, \ldots, \Theta) \\\\ &= P(y_0)\Pi _{k=1}^NP(y_k| Y_{k-1},U_N,\Theta) \\\\ &= P(y_0)\Pi _{k=1}^NP(e_k) \\\\ &= P(y_0)(2\pi)^{-N/2}\sigma^{-N} \exp \left[ -(1/2\sigma^2)\sum_{k=1}^N e_k^2\right] \tag{1} \end{array} P(YN∣UN,Θ)=P(yN,…∣uN,…,Θ)=P(y0)Πk=1NP(yk∣Yk−1,UN,Θ)=P(y0)Πk=1NP(ek)=P(y0)(2π)−N/2σ−Nexp[−(1/2σ2)∑k=1Nek2](1)
其中
ek=yk−ϕiΘϕk=[−yk−1,…,−yk−n∣uk−1,…,uk−n∣ek−1,…,ek−n]⊤Θ=[a1,…,an∣bi,…,bn∣c1,…,cn]⊤\begin{array}{ll} e_k &= y_k - \phi_i \Theta \\\\ \phi_k & = [-y_{k-1}, \ldots, -y_{k-n} | u_{k-1}, \ldots, u_{k-n} | e_{k-1}, \ldots, e_{k-n}]^\top \\\\ \Theta &= [a_1, \ldots, a_n | b_i, \ldots, b_n | c_1, \dots, c_n ]^\top \end{array} ekϕkΘ=yk−ϕiΘ=[−yk−1,…,−yk−n∣uk−1,…,uk−n∣ek−1,…,ek−n]⊤=[a1,…,an∣bi,…,bn∣c1,…,cn]⊤
最大化似然(1)等价于最小化负对数似然(2)
J(σ,Θ)=Nln⁡σ+12σ2∑kNek2(2)J(\sigma, \Theta) = N \ln \sigma + \frac{1}{2\sigma^2}\sum_k^N e_k^2 \tag{2} J(σ,Θ)=Nlnσ+2σ21k∑Nek2(2)

辨识过程

在实际辨识过程中，由于参数 Θ\ThetaΘ 未知，误差 eke_kek 不能精确获得，所以计算过程中用其估计值 νk≃ek\nu_k \simeq e_kνk≃ek 来替代。

代价函数：
J(σν,Θ)=Nln⁡σν+12σν2∑kNνk2(2)J(\sigma_\nu, \Theta) = N \ln \sigma_{\nu} + \frac{1}{2\sigma_\nu^2}\sum_k^N \nu_k^2 \tag{2} J(σν,Θ)=Nlnσν+2σν21k∑Nνk2(2)

采集 NNN 组数据，估计一个初始 Θ0\Theta_0Θ0，比如先假定误差项系数 ci=0c_i=0ci=0，用最小二乘法求解ai,bia_i, b_iai,bi；
k=0k = 0k=0
固定 Θ\ThetaΘ（即固定ν\nuν），更新 σν\sigma_\nuσν:
σν=arg min⁡σνJ=∑kνk2/N\sigma_\nu = \argmin_{\sigma_\nu} J = \sqrt{\sum_k \nu_k^2/N} σν=σνargminJ=k∑νk2/N
固定σν\sigma_\nuσν，更新 Θ\ThetaΘ，即最小化
L=∑kνk2L = \sum_k \nu_k^2 L=k∑νk2
采用牛顿法更新参数：
Θt+1=Θt−H−1∇ΘL\Theta_{t+1} = \Theta_t - H^{-1} \nabla_{\Theta} L Θt+1=Θt−H−1∇ΘL
k=k+1k = k +1k=k+1，重复以上交替优化过程直到
σt2−σt−12σt−12<10−4\frac{\sigma_t^2 - \sigma_{t-1}^2}{\sigma_{t-1}^2} < 10^{-4} σt−12σt2−σt−12<10−4

关于参数的梯度与海森矩阵

梯度

∂L∂Θ=2∑νk∂νk∂Θ(3)\frac{\partial L}{\partial \Theta} = 2\sum \nu_k \frac{\partial \nu_k}{\partial \Theta} \tag{3} ∂Θ∂L=2∑νk∂Θ∂νk(3)
这里稍稍有一点复杂，因为：
νk=yk+∑i=1naiyk−i−∑i=1nbiuk−i−∑i=1nciνk−i=yk−ϕk⊤Θ\nu_k = y_k + \sum_{i=1}^n a_iy_{k-i} - \sum_{i=1}^n b_i u_{k-i} - \sum_{i=1}^n c_i \nu_{k-i} = y_k - \phi_k ^\top\Theta νk=yk+i=1∑naiyk−i−i=1∑nbiuk−i−i=1∑nciνk−i=yk−ϕk⊤Θ
所以
∂νk∂ai=yk−i−∑l=1ncl∂νk−l∂ai∂νk∂bi=−uk−i−∑l=1ncl∂νk−l∂bi∂νk∂ci=−νk−i−∑l=1ncl∂νk−l∂ci\frac{\partial \nu_k}{\partial a_i} = y_{k-i} - \sum_{l=1}^n c_l \frac{\partial \nu_{k-l}}{\partial a_i} \\ \frac{\partial \nu_k}{\partial b_i} = -u_{k-i} - \sum_{l=1}^n c_l \frac{\partial \nu_{k-l}}{\partial b_i} \\ \frac{\partial \nu_k}{\partial c_i} = -\nu_{k-i}- \sum_{l=1}^n c_l \frac{\partial \nu_{k-l}}{\partial c_i} ∂ai∂νk=yk−i−l=1∑ncl∂ai∂νk−l∂bi∂νk=−uk−i−l=1∑ncl∂bi∂νk−l∂ci∂νk=−νk−i−l=1∑ncl∂ci∂νk−l
可以合并成
∂νk∂Θ=−ϕk⊤−∑l=1ncl∂νk−l∂Θ\frac{\partial \nu_k}{\partial \Theta} = -\phi_{k}^\top - \sum_{l=1}^n c_l \frac{\partial \nu_{k-l}}{\partial \Theta} ∂Θ∂νk=−ϕk⊤−l=1∑ncl∂Θ∂νk−l
所以∂L/∂Θ\partial L/\partial \Theta∂L/∂Θ 已求得！

海森矩阵

假定 ∂νk/∂Θ\partial \nu_k/\partial \Theta∂νk/∂Θ 为行向量
H=∂2L∂Θ2=∂∂Θ(2∑νk∂νk∂Θ)=2(∑k(∂νk∂Θ)⊤(∂νk∂Θ)+νk∂2νk∂Θ2)H = \frac{\partial ^2L}{\partial \Theta^2} = \frac{\partial }{\partial \Theta}\left(2\sum \nu_k \frac{\partial \nu_k}{\partial \Theta} \right) \\ = 2\left( \sum_k \left(\frac{\partial \nu_k}{\partial \Theta}\right)^\top \left( \frac{\partial \nu_k}{\partial \Theta}\right) + \nu_k \frac{\partial^2 \nu_k}{\partial \Theta^2}\right) H=∂Θ2∂2L=∂Θ∂(2∑νk∂Θ∂νk)=2(k∑(∂Θ∂νk)⊤(∂Θ∂νk)+νk∂Θ2∂2νk)
主要需要考虑 ∂2ν/∂Θ2\partial^2 \nu/ \partial \Theta^2∂2ν/∂Θ2:
∂2νk∂ai∂aj=−∑l=1ncl∂2νk−l∂ai∂aj∂2νk∂ai∂bj=−∑l=1ncl∂2νk−l∂ai∂bj∂2νk∂ai∂cj=−∂νk−j∂ai−∑l=1ncl∂2νk−l∂ai∂cj∂2νk∂bi∂cj=−∂νk−j∂bi−∑l=1ncl∂2νk−l∂ai∂cj∂2νk∂ci∂cj=−2∂νk−j∂ci−∑l=1ncl∂2νk−l∂ci∂cj\frac{\partial^2 \nu_k}{ \partial a_i \partial a_j} = - \sum_{l=1}^n c_l \frac{\partial^2 \nu_{k-l}}{\partial a_i\partial a_j} \\ \frac{\partial^2 \nu_k}{ \partial a_i \partial b_j} = - \sum_{l=1}^n c_l \frac{\partial^2 \nu_{k-l}}{\partial a_i\partial b_j} \\ \frac{\partial^2 \nu_k}{ \partial a_i \partial c_j} = -\frac{\partial \nu_{k-j}}{\partial a_i} - \sum_{l=1}^n c_l \frac{\partial^2 \nu_{k-l}}{\partial a_i\partial c_j} \\ \frac{\partial^2 \nu_k}{ \partial b_i \partial c_j} = -\frac{\partial \nu_{k-j}}{\partial b_i} - \sum_{l=1}^n c_l \frac{\partial^2 \nu_{k-l}}{\partial a_i\partial c_j} \\ \frac{\partial^2 \nu_k}{ \partial c_i \partial c_j} = -2\frac{\partial \nu_{k-j}}{\partial c_i} - \sum_{l=1}^n c_l \frac{\partial^2 \nu_{k-l}}{\partial c_i\partial c_j} \\ ∂ai∂aj∂2νk=−l=1∑ncl∂ai∂aj∂2νk−l∂ai∂bj∂2νk=−l=1∑ncl∂ai∂bj∂2νk−l∂ai∂cj∂2νk=−∂ai∂νk−j−l=1∑ncl∂ai∂cj∂2νk−l∂bi∂cj∂2νk=−∂bi∂νk−j−l=1∑ncl∂ai∂cj∂2νk−l∂ci∂cj∂2νk=−2∂ci∂νk−j−l=1∑ncl∂ci∂cj∂2νk−l
写成矩阵形式：
Hk=−∑l=1nclHk−l−Gk−Gk⊤Gk=[0…0∣0…0∣(∂νk−1∂Θ)⊤…(∂νk−n∂Θ)⊤]3n×3nH_k = -\sum_{l=1}^n c_lH_{k-l} - G_k -G_k^\top \\ G_k = \left[ 0 \ldots 0 \bigg| 0 \ldots 0 \bigg| \left(\frac{\partial \nu_{k-1}}{\partial \Theta}\right)^\top \ldots \left(\frac{\partial \nu_{k-n}}{\partial \Theta}\right)^\top \right]_{3n\times 3n} Hk=−l=1∑nclHk−l−Gk−Gk⊤Gk=[0…0∣∣∣∣0…0∣∣∣∣(∂Θ∂νk−1)⊤…(∂Θ∂νk−n)⊤]3n×3n
再次强调一下：∂νk/∂Θ\partial \nu_{k}/\partial \Theta∂νk/∂Θ 是行向量，因为 Θ\ThetaΘ 是列向量！

到此，完整的计算过程应该心中有数了吧！

受控自回归滑动平均模型（ARMAX）的系统辨识相关推荐

运用自回归滑动平均模型、灰色预测模型、BP神经网络三种模型分别预测全球平均气温，并进行预测精度对比（附代码、数据）
大家好,我是带我去滑雪,每天教你一个小技巧!全球变暖是近十年来,人们关注度最高的话题.2022年夏天,蔓延全球40℃以上的极端天气不断刷新人们对于高温的认知,人们再也不会像从前那样认为全球变暖离我们遥 ...
基于自回归整合滑动平均模型(ARIMA)的时间序列预测
基于自回归整合滑动平均模型(ARIMA)的时间序列预测 ID:8629644191951129誩宝
TensorFlow模型保存和提取方法（含滑动平均模型）
一.TensorFlow模型保存和提取方法 1. TensorFlow通过tf.train.Saver类实现神经网络模型的保存和提取.tf.train.Saver对象saver的save方法将Tens ...
tensorflow 滑动平均模型 ExponentialMovingAverage
____tz_zs学习笔记滑动平均模型对于采用GradientDescent或Momentum训练的神经网络的表现都有一定程度上的提升. 原理:在训练神经网络时,不断保持和更新每个参数的滑动平均值, ...
tensorflow中滑动平均模型的说明
作用:使用滑动平均模型可以使模型在测试数据上更健壮.即如果在测试过程中,出现了一些噪声数据,滑动平均模型可以很好地应对这些数据,使这些噪声数据不会对模型的变量造成太大的影响. 1.滑动平均模型原理: ...
深入解析TensorFlow中滑动平均模型与代码实现
因为本人是自学深度学习的,有什么说的不对的地方望大神指出指数加权平均算法的原理 TensorFlow中的滑动平均模型使用的是滑动平均(Moving Average)算法,又称为指数加权移动平均算法( ...
滑动平均模型（MA）—tensorflow
在采用梯度下降的方式训练神经网络的时候,我们使用滑动平均模型会在一定的程度上提高最终模型在测试集上的表现. 在TensorFlow中提供了tf.train.ExponentialMovingAvera ...
tensorflow tf.train.ExponentialMovingAverage() （滑动平均模型）（移动平均法 Moving average，MA）（用于平滑数据波动对预测结果的影响）
tf.train.ExponentialMovingAverage 函数定义 tensorflow中提供了tf.train.ExponentialMovingAverage来实现滑动平均模型,他使用指 ...
深度学习中滑动平均模型的作用、计算方法及tensorflow代码示例
滑动平均模型: 用途:用于控制变量的更新幅度,使得模型在训练初期参数更新较快,在接近最优值处参数更新较慢,幅度较小方式:主要通过不断更新衰减率来控制变量的更新幅度衰减率计算公式 : dec ...

受控自回归滑动平均模型（ARMAX）的系统辨识