数学狂想曲（十一）——平稳离散时间随机过程, 功率谱, 高阶统计
2018-09-28平稳离散时间随机过程
Toeplitz矩阵
Toeplitz矩阵（diagonal-constant matrix），指矩阵中每条自左上至右下的斜线上的元素相同。
Otto Toeplitz，1881～1940，德国犹太裔数学家。University of Breslau博士（1905），先后执教于Göttingen University（在David Hilbert手下供职）、University of Kiel和Bonn University。1939年，为了躲避元首的迫害，逃亡耶路撒冷，次年去世。
广义平稳离散时间随机过程的相关矩阵是Hermite矩阵，也是Toeplitz矩阵。反之，如果相关矩阵是Toeplitz矩阵，则该离散时间随机过程，一定是广义平稳的。
离散时间随机过程的相关矩阵是非负定的，并且几乎总是正定的。（等于零，只有在无噪声且观测向量线性相关的情况下，才会出现。）
白噪声
\[E[v(n)v^*(n-k)]=\begin{cases} \sigma_v^2, & k = 0 \\ 0, & k \neq 0 \\ \end{cases}\]线性差分方程
时间随机过程本身是由时间序列组成的，因此也可以使用《机器学习（三十七）》中提到的ARIMA模型。该模型的关键是求解线性差分方程。这通常要使用“信号与系统”课程中的z变换（离散域的拉普拉斯变换）求解。考虑到“信号与系统”是一个很大的课程。这里仅对本人关心的要点，做一个简要记录。
绝对可积->收敛域
z变换：\(f(z)\to F(z)\)
z逆变换：\(F(z)\to f(z)\)
系统函数：\(H(z)=\frac{R(z)}{E(z)}\)。其中，E是激励信号，R是系统响应。
E的收敛域：\(\mid z\mid >1\)
差分算子->特征方程->特征根
H的平稳条件：H的特征根满足\(\mid z\mid \le 1\)
特征根是正实数，且\(\mid z \mid<1\)：自相关函数为阻尼曲线，仅有幅变。
特征根是负实数或者复数，且\(\mid z \mid<1\)：自相关函数为正弦阻尼曲线，不仅有幅变，还有相变。
选择ARIMA的阶数
如前所述，ARIMA(p,d,q)除了一些参数之外，还包括p，d，p这三个阶数的超参数。
AIC信息准则即Akaike information criterion，是衡量统计模型拟合优良性(Goodness of fit)的一种标准，由于它为日本统计学家赤池弘次创立和发展的，因此又称赤池信息量准则。AIC方法主要使用了KL散度。
MDL(minimum description length,最小描述长度) 原理是Rissane在研究通用编码时提出的。其基本原理是选择总描述长度最小的模型。
参考：
https://mp.weixin.qq.com/s/66lY17sOO83Q-xhvQi72dw
周期性时间序列的预测
功率谱
随机过程（设时间序列为\(u(n)\)）二阶统计：
时域——自相关函数：
\[r_N(n-k)=E[u_N(n)u_N^*(k)]\tag{1}\]其中，\(u_N^*(k)\)是\(u_N(k)\)的复共轭。
频域：
\[U_N(\omega)=\sum_{n=-N}^Nu_N(n)e^{-j\omega n}\tag{2}\] \[S(\omega)=\lim_{N\to\infty}\frac{1}{N}E[\mid U_N(\omega)\mid^2]=\sum_{l=-\infty}^{+\infty}r(l)e^{-j\omega l}\tag{3}\]其中，\(S(\omega)\)就是功率谱密度（power spectral density, PSD），也称为功率谱（power spectrum）。
自相关函数和功率谱密度组成了傅立叶变换对，这种关系又被称为EWK（Einstein-Wiener-Khintchine）关系。
Einstein最早提出idea，Wiener证明了一个特例，Khintchine做了扩展证明。
Aleksandr Yakovlevich Khinchin，1894～1959，苏联数学家。莫斯科州立大学毕业，并留校任教，直到去世。苏联概率学派的重要人物。苏联科学院院士。概率论中，著名的Khintchine inequality就是他的成果。
在频域上，我们有Nyquist频率，相应的在时域上，我们也有Nyquist间隔：在这个间隔之外，\(S(\omega)\)是周期性的。
离散时间随机过程的功率谱密度是非负实函数。
\[S_o(\omega)=\mid H(e^{j\omega})\mid S(\omega)\tag{4}\]其中，H为系统函数，\(S_o\)输出信号的功率谱密度。
功率谱密度的Cramér表示：
\[u(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{j\omega n}\mathrm{d}Z(\omega)\tag{5}\]其中，\(\mathrm{d}Z(\omega)\)被称为增量过程（increment process）。
Harald Cramér，1893～1985，瑞典数学家、统计学家。Stockholm University博士（1917）、教授、校长、瑞典高等教育系统大臣。被誉为“统计理论的巨人”。
由公式2和5，可得：
\[U_N(\omega)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\sum_{n=-N}^N e^{(-j(\omega-v) n)}\mathrm{d}Z(v)\tag{6}\]我们定义：
\[K_N(\omega)=\sum_{n=-N}^N e^{j\omega n}=\frac{\sin((2N+1)\omega/2)}{\sin(\omega/2)}\tag{7}\]则公式6可改写为：
\[U_N(\omega)=\frac{1}{2\pi}\int_{-\pi}^{\pi}K_N(\omega-v)\mathrm{d}Z(v)\tag{8}\]这里的K被称作Dirichlet Kernel。参见《数学狂想曲（一）》的相关章节。
一般来说，在公式8中，\(U_N(\omega)\)是已知的，而\(\mathrm{d}Z(\omega)\)是未知的。从数学上来说，这个积分方程可看做第一类Fredholm积分方程的一个例子。
Erik Ivar Fredholm，1866～1927，瑞典数学家。Uppsala University博士（1898）+Stockholm University教授。不知道是不是瑞典的保险业比较发达，他和Cramér居然都当过兼职的精算师。。。瑞典皇家科学院院士。
Uppsala University是瑞典，也是北欧，最古老的大学，始建于1477年。
功率谱密度的估计方法主要包括参数法和非参数法两大类。
参数法包括：
1.模型辨识法。基本就是上面提到的ARIMA或者其变种。
2.最小方差无失真响应法（MVDR）。
3.特征分解法。将相关矩阵R分解为两个子空间：信号子空间和噪声子空间。
非参数法包括：
1.周期图法。
2.多窗口法。
一般来说，随机过程的功率谱包含两个分量：确定性分量和连续分量。前者是增量过程\(\mathrm{d}Z(\omega)\)的一阶矩，后者是\(\mathrm{d}Z(\omega)\)的二阶中心矩。
参数法一般在知道相关物理规律时使用，它具有较高的精确度。而非参数法由于只依赖增量过程的一阶矩和二阶中心矩，因此适用范围更广泛，即使不知道系统的物理规律也可以使用。（有些类似万能拟合的GMM）
参考：
https://www.zhihu.com/question/29520851
功率谱密度如何理解？
https://mp.weixin.qq.com/s/i84yW8L7qGVpPtwpaJhhgw
深入浅出聊抖动（Jitter）
https://mp.weixin.qq.com/s/XZHh7O-aiTA4Z-jnz20HgA
功率谱密度
高阶统计
上面讨论的基本都是一阶和二阶统计量，实际上我们还可以使用更高阶的统计量。使用高阶统计量的学科，一般被称为高阶统计学（higher-order statistics）。
Moment
Moment（矩）的定义为：
\[\mu_n = \int_{-\infty}^\infty (x - c)^n\,f(x)\,\mathrm{d}x\]其中，当c=0时，被称作Raw Moment。当c为均值时，被称作Central Moment。如果用\(\mu_n/\sigma^n\)替换\(\mu_n\)，就是所谓的Normalised Moment（标准矩）了。
1阶Raw Moment，常称为Mean。2阶Central Moment，常称为Variance（方差）。3阶Normalised Moment，常称为Skewness（偏态）。4阶Normalised Moment，常称为Kurtosis（峰度）。
Cumulants
Cumulants（累积量）的思想最早是Thorvald Thiele提出的，后来被Ronald Fisher和John Wishart发扬光大。
Thorvald Nicolai Thiele，1838～1910，丹麦天文学家。哥本哈根大学博士。哥本哈根天文台台长（1978～1907）。曾研究过三体问题。被Ronald Fisher誉为“最伟大的统计学家”。
John Wishart，1898～1956，苏格兰数学家和农业统计学家。Edinburgh University本科+Cambridge University硕士+University College London博士。导师是Karl Pearson，和Ronald Fisher也有过合作。Royal Society of Edinburgh会员。Cambridge University统计实验室首任主任。
苏格兰人的自我意识真是强，足球有自己的协会，就连皇家学会也有自己的。
在介绍Cumulants之前，我们首先看一下Moment-generating function：
\[M_X(t) := \operatorname E \left[e^{tX}\right], \quad t \in \mathbb{R}\]可以看出，MGF和《数学狂想曲（二）》中提到的随机变量的特征函数（Characteristic function, CF）的形式非常类似。
而cumulant-generating function则是MGF的对数，即：
\[K(t)=\log\operatorname{E}\left[e^{tX}\right]\]对上式进行Maclaurin展开，可得：
\[K(t)=\sum_{n=1}^\infty \kappa_{n} \frac{t^{n}}{n!} = \mu t + \sigma^2 \frac{t^2}{2} + \cdots\]这里的\(\kappa_{n}\)就是Cumulants了。
由MGF和CF的关系易知，使用CF的对数的Maclaurin展开，也可以求Cumulants。
Cumulants有如下性质：
\[cum(\lambda_1x_1,\dots,\lambda_kx_k)=\sum_{i=1}^k\lambda_i cum(x_1,\dots,x_k)\tag{1}\]其中，\(\lambda_i\)为常数。
\[cum(x_1,\dots,x_k)=cum(x_{i_1},\dots,x_{i_k})\tag{2}\]其中，\((i_1,\dots,i_k)\)为\((1,\dots,k)\)的任意一种排列。
\[cum(x_0+y_0,z_1,\dots,z_k)=cum(x_0,z_1,\dots,z_k) + cum(y_0,z_1,\dots,z_k)\tag{3}\]如果\(\alpha\)为常数，则：
\[cum(\alpha+z_1,\dots,z_k)=cum(z_1,\dots,z_k)\tag{4}\]如果\(x_i\)与\(y_i\)相互独立，则：
\[cum(x_1+y_1,\dots,x_k+y_k)=cum(x_1,\dots,x_k)+cum(y_1,\dots,y_k)\tag{5}\]Polyspectrum
Polyspectrum一般翻译成高阶谱或者多谱。
吐槽一下，我最早看的一本书把Polyspectrum翻译成多谱，结果我以此为关键词搜索，基本一无所获。直到我发现它还有另一个中文名。。。
以上讨论的是Cumulants的一般定义，在处理离散时间随机过程的时候，我们往往采用如下形式的定义。
参考
https://www.zhihu.com/question/25344430
随机变量的矩和高阶矩有什么实在的含义？
https://www.zhihu.com/question/43469699
信号的矩和高阶累积量的定义是什么？
http://www.doc88.com/p-1127198771359.html
高阶累积量与高阶谱读书笔记
https://wenku.baidu.com/view/7c4931085727a5e9856a6139.html
高阶谱分析
https://wenku.baidu.com/view/136b666c561252d380eb6e8c.html
高阶统计量的定义与性质
https://mp.weixin.qq.com/s/SOpBf0fAoqau0ac5kXlQIQ
数学好吃吗？好吃！
