小样本马尔可夫链独立性检验
Published:
本文主要包括以下内容：
- 案例构建和样本采集
- 小样本量下的独立性检验
一、案例构建和样本采集
假定有两个独立的骰子 $X$ 和 $Y$，每个骰子有6个面。独立重复 $N_{\text{trails}}$ 次掷骰子试验，每次试验持续较少的 $N_{\text{steps}}$ 步（即小样本量）。分别记录第i次试验的第j步骰子结果为 $x_{i,j}$ 和 $y_{i,j}$。最终获得样本数组 $X_{N_{\text{trails}}×N_{\text{steps}}}$ 和 $Y_{N_{\text{trails}}×N_{\text{steps}}}$ 用于分析。
在本案例中，使用状态转移矩阵 $\pi$ 控制骰子 $X$ 和 $Y$ 各自的点数变化。
状态转移矩阵性质：
- 每个元素 $π_{i,j} \geq 0$ 表示从当前状态 $i$ 转移到下个状态 $j$ 的概率
- $\boldsymbol{\pi}$ 的每一行之和必为1
- 骰子点数为1至6，状态空间维数为6
- $\boldsymbol{\pi}$ 为 $6 \times 6$ 的矩阵
独立投掷的状态转移矩阵：
\[\boldsymbol{\pi} = \frac{1}{6} \cdot \mathbf{1}_{6 \times 6}\]1阶投掷的状态转移矩阵：
\[\boldsymbol{\pi} = \begin{bmatrix} 0.5 & 0.25 & 0 & 0 & 0 & 0.25 \\ 0.25 & 0.5 & 0.25 & 0 & 0 & 0 \\ 0 & 0.25 & 0.5 & 0.25 & 0 & 0 \\ 0 & 0 & 0.25 & 0.5 & 0.25 & 0 \\ 0 & 0 & 0 & 0.25 & 0.5 & 0.25 \\ 0.25 & 0 & 0 & 0 & 0.25 & 0.5 \end{bmatrix}\]2阶投掷的状态转移矩阵：
对于任意历史状态 $(s_{t-2}, s_{t-1}) \in \mathcal{S}^2$：
\[\pi(X_t \mid X_{t-2}, X_{t-1}) = \begin{cases} [0.7, 0.2, 0.1] & \text{if } (0,0) \\ [0.1, 0.6, 0.3] & \text{if } (0,1) \\ [0.2, 0.2, 0.6] & \text{if } (0,2) \\ [0.3, 0.4, 0.3] & \text{if } (1,0) \\ [0.1, 0.8, 0.1] & \text{if } (1,1) \\ [0.0, 0.1, 0.9] & \text{if } (1,2) \\ [0.5, 0.5, 0.0] & \text{if } (2,0) \\ [0.2, 0.3, 0.5] & \text{if } (2,1) \\ [0.1, 0.1, 0.8] & \text{if } (2,2) \end{cases}\]对应样本变化曲线如下：
可见，随着阶数的增加，样本序列的变化趋势逐渐平缓。因此，可以通过设计状态转移的方式对每次试验的点数变化规律进行控制。需注意，由于骰子 $X$ 和 $Y$ 相互独立，不论状态转移矩阵如何设置，$X$ 和 $Y$ 试验所得样本都应无关。
本文研究：如何使用小样本，准确检验 $X$ 和 $Y$ 的独立性，并确定马尔可夫链的阶数。
二、小样本量下的独立性检验
2.1 非马尔可夫链的独立性检验方法
如果 $X$ 和 $Y$ 的样本序列不具有马尔可夫性，则理论上在任意试验 $i$ 中，$X$ 和 $Y$ 应始终相互独立。因此，对每组样本 $x_{i}=\left[ x_{i,1},\cdots,x_{i,N_{\text{steps}}} \right]$ 和 $y_{i}=\left[ y_{i,1},\cdots,y_{i,N_{\text{steps}}} \right]$ 进行独立性检验时，结果应显示二者独立。需要注意的是，由于每次试验的步数 $N_{\text{steps}}$ 较少（实际应用中常见），必须确保在如此小样本量下，独立性检验结果依然具有可靠性。为此，置换检验和 Bootstrap 检验等非参数方法特别适合用于小样本量下复杂分布数据的独立性检验。具体而言，对于某一 $X$-$Y$ 联合分布（其来源不一定为马尔可夫链），通过独立同分布采样获得的样本 $x_{i}$ 和 $y_{i}$，可先计算互信息 $I(x_{i};y_{i})$ 以衡量其关联度，再通过置换检验进行判断：若 $I(x_{i};y_{i})$ 未显著大于 0（零假设 $H_{0}$），则接受 $X$ 与 $Y$ 无关的结论；反之，若显著大于 0（备择假设 $H_{1}$），则认为两者存在关联。
互信息是一种衡量随机变量间线性和非线性关联程度的指标，定义为：
\[I(X;Y) = \sum_{x \in X} \sum_{y \in Y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)} \tag{1}\]2.1.1 排列置换检验
在零假设 $H_0$ 成立条件下（即 $X$ 和 $Y$ 样本来自同一总体），其组别标签可随机交换而不影响统计量分布。具体实施时：首先将原始样本合并为总集 $Z = X \cup Y = {z_1,\cdots,z_{N_X+N_Y}}$，随后进行 $N_{\text{perm}}$ 轮置换抽样，每轮从 $Z$ 中无放回抽取 $\lfloor(N_X+N_Y)/2\rfloor$ 个样本作为置换组 $X_k^{\text{perm}}$，剩余样本作为 $Y_k^{\text{perm}}$，并计算每轮的关联系数$I(X_k^{\text{perm}};Y_k^{\text{perm}})$；通过重复该过程构建经验分布函数
\[F(t) = \frac{1}{N_{\text{perm}}}\sum_{k=1}^{N_{\text{perm}}} \mathbb{1}\left(I(X_k^{\text{perm}};Y_k^{\text{perm}}) \leq t\right) \tag{2}\]最终，结合样本实际观测值计算 $p$ 值，并与单边检验显著性水平 $\alpha$ 比较：
\[p = \frac{\sum_{k=1}^{N_{\text{perm}}} \mathbb{1}\left(I(X_k^{\text{perm}};Y_k^{\text{perm}}) \geq I_{\text{obs}}\right)}{N_{\text{perm}}} \tag{3}\]若 $p \leq \alpha$，则拒绝零假设 $H_0$ 接受 $H_1$，认为 $X$ 和 $Y$ 之间存在显著关联。
2.1.2 蒙特卡洛置换检验
蒙特卡洛置换检验提供了一种非参数的关联性分析方法，其流程与排列置换检验类似，但在各变量内部进行置换操作。对于一组 样本 $X$ 和 $Y$，首先计算其互信息 $I(X;Y)$。然后进行 $N_{\text{perm}}$ 次随机抽样，每次随机打乱 $X$ 顺序得到置换样本 $X_k^{\text{perm}}$，并计算互信息 $I(X_k^{\text{perm}};Y)$。最终通过多轮置换获得经验分布，计算 $p$ 值和显著性。
\[F(t) = \frac{1}{N_{\text{perm}}}\sum_{k=1}^{N_{\text{perm}}} \mathbb{1}\left(I(X_k^{\text{perm}};Y) \leq t\right) \quad \tag{4}\] \[p = \frac{\sum_{k=1}^{N_{\text{perm}}} \mathbb{1}\left(I(X_k^{\text{perm}};Y) \geq I_{\text{obs}}\right)}{N_{\text{perm}}} \tag{5}\]2.1.3 Bootstrap 检验
Bootstrap 检验是一种基于重采样的非参数方法，适用于小样本量下的独立性检验。其基本思想是通过对原始样本进行有放回抽样，构建多个 Bootstrap 样本集，从而估计统计量的分布。对于一组样本 $X$ 和 $Y$，首先计算其互信息 $I(X;Y)$。然后进行 $N_{\text{bootstrap}}$ 次随机抽样，每次从 $X$ 中有放回地抽取 $N_X$ 个样本，得到 Bootstrap 样本 $X_k^{\text{bootstrap}}$。最终通过多轮 Bootstrap 抽样获得经验分布，计算 $p$ 值和显著性。
\[F(t) = \frac{1}{N_{\text{bootstrap}}}\sum_{k=1}^{N_{\text{bootstrap}}} \mathbb{1}\left(I(X_k^{\text{bootstrap}};Y) \leq t\right) \quad \tag{6}\] \[p = \frac{\sum_{k=1}^{N_{\text{bootstrap}}} \mathbb{1}\left(I(X_k^{\text{bootstrap}};Y) \geq I_{\text{obs}}\right)}{N_{\text{bootstrap}}} \tag{7}\]2.2 马尔可夫链的独立性检验方法
2.1节中介绍的3种方法均适用于非马尔可夫链的独立性检验。然而，当样本序列具有马尔可夫性时，序列具有时序性和马尔可夫性，直接随机打乱样本会破坏内部的时序结构和依赖关系，丢失本底关联信息，进而造成检验结果失真。因此，针对马尔可夫链的独立性检验，需要采用更复杂的检验方法。
2.1节方法用于马尔可夫链检验的问题主要在于零假设分布样本构建失真，需要采用如下的保序替代样本（order-preserving surrogates）构建方法：
- 设定变量数据的马尔可夫阶数 $k$，并基于原始序列 $X$ 估计状态转移矩阵$\hat{\boldsymbol{\pi}}$；
- 利用估计的转移概率矩阵$\hat{\boldsymbol{\pi}}$，从 $X$ 的初始状态出发，随机生成保序替代样本 $X_{k}^{\text{surrog}}$，并计算其互信息 $I(X_{k}^{\text{surrog}};Y)$；
- 重复步骤2，进行 $N_{\text{surrogates}}$ 次生成计算，得到保序替代样本互信息值的经验分布；
- 计算原始样本的互信息 $I(X;Y)$，并与保序替代样本的经验分布进行比较，计算 $p$ 值和显著性。
三、案例分析
3.1 独立投掷样本序列
下图展示了 $N_{\text{trials}}=1000$ 组独立投掷样本序列 $X$ 和 $Y$ 所得到的真实互信息分布，以及仅通过第1组100个样本在 $N_{\text{resample}}=1000$ 次重采样参数下通过置换检验（perm）、蒙特卡洛置换检验（MC_perm）、Bootstrap 检验（Bootstrap）和马尔可夫链检验（Markov_chain）获得的零假设互信息分布。
可以看出，所有检验方法得到的零假设分布均与真实分布高度一致，说明这四种方法均能有效地检验非马尔可夫数据的独立性。
3.2 一阶投掷样本序列
下图继续展示了四种方法在1阶投掷样本序列 $X$ 和 $Y$ 上的零假设互信息分布与真实分布对比。
这次情况发生了变化，置换检验、蒙特卡洛置换检验和 Bootstrap 检验的零假设分布与真实分布存在明显负向偏差，说明这三种方法在1阶马尔可夫链数据上失效，且更容易拒绝掉实际无关的MI值，接受 $H_1$ 备择假设，导致更高的I类错误率。相比之下，只有马尔可夫链检验所得的零假设分布与真实分布高度一致，说明其能够有效地检验1阶马尔可夫数据的独立性。
综上，只有马尔可夫链检验能够同时胜任非马尔可夫链和马尔可夫链数据的独立性检验，而其他三种方法在马尔可夫链数据上会出现较高的I类错误率。在算法实现时，可将非马尔可夫链数据的阶数设为0代入计算，提升通用性。