因果语言模型和半因果语言模型 因果关系模型有哪些

因果模型二：线性非高斯无环模型

• 一、前提条件
• 二、方程形式
• 三、算法过程
• 四、检验

• 1. 剪枝检验（检验一）
• 2. 模型拟合度检验（检验二）
• 3. 嵌套模型的差异卡方检验（检验三）
• 4. 综合三种检验方法的剪枝检验算法

• 五、算法效果验证

• 1. 实验数据检验算法准确性
• 2. 真实数据检验算法适用性

本篇主要介绍一种基于贝叶斯网络的、具体化的求解因果关系的模型：线性非高斯无环模型（LiNGAM）。主要是让我们对如何把因果关系研究抽象化为一个数学模型，以及如何求解它先有一个初步的认识。

一、前提条件

线性非高斯无环模型要求三个基本的前提条件：

1. 观测变量之间是存在因果顺序的，后续变量不会导致前序变量；这些变量可以用一个DAG（有向无环图）表示。
2. 变量间的因果关系是线性的，可以用如下公式描述：
   $x_{i}=\sum_{k(j)<k(i)}b_{i,j}x_{j}+e_{i}+c_{i}$
3. 扰动项e服从方差非零的非高斯分布，且互相独立。

二、方程形式

$X = BX+e$

其中X，B，e都是矩阵，分别表示样本矩阵、系数矩阵和扰动项矩阵。如果我们能够按变量之间的因果顺序将X进行排序，则系数矩阵B可以排列成为一个严格意义上的下三角矩阵（对角线全为0值）。我们以一个具体的例子简单阐述：

假设一个DAG如下图所示：

其变量间的因果关系都是线性的，则整个DAG的因果关系方程如下：

$x_{2} = b_{21} * x_{1} + e_{2}\\ x_{3} = b_{32} * x_{2} + e_{3}\\ x_{4} = b_{42} * x_{2} + b_{43} * x_{3} + e_{4}$

将上述方程写成矩阵的形式，即

即$X = BX+e$的具体展开形式。

对方程中的X进行求解，可得到如$X=Ae$的形式，其中$A = (I-B)^{-1}$。在得知了变量X的因果顺序之后，B可以排列成一个下三角矩阵，那A同样也可以排列成一个下三角矩阵。（参考下图，A矩阵中右上三角部分的元素一定全是0）

当方程被抽象到这一步的时候，因果关系求解的大致脉络也就清晰了，我们只需要找到一种算法，首先把X矩阵分解成A和e两部分，然后再找到一种算法，吧A矩阵排序成一个下三角矩阵，就可以把B求解出来了，整个因果关系图也就得到了。接下来我们看一下算法步骤：

三、算法过程

1. 样本矩阵X（m X n, m<<n），每一列代表一个样本向量，首先每一行去均值，然后使用ICA算法（比如fastICA）得到样本矩阵X的一个分解式：X = AS，其中S和A同大小，每一行存放独立项。

• 这里需要解释一下什么是ICA算法：ICA算法常用于信号分离，其作用就是将一个混杂了多种信号的信息分解为一个个独立的信号。以FastICA算法为例，该算法寻找方向使得随机变量的某种“非高斯性”(non-Gaussianity)的度量最大化。一种常用的非高斯性的度量是四阶矩E[(WTX)4]。类似PCA的流程，我们首先找w1使得E[(w1TX)4]最大；然后在与w1正交的空间里找w2，使得E[(wTX)4]最大，以此类推直到找到所有的w1,w2,w3,…,wn. 可以证明，用这种方法得到的是相互独立的w1Tz,w2Tz,w3Tz,…,wnTz。这里可以证明，只要源信号是满足非高斯分布的，那么这种分解就是唯一的。（参考这篇博客和这篇文献）。这也是为什么LiNGAM模型要求数据要服从非高斯分布。

2. 在$W（W = A^{-1}）$所有排序矩阵中，找到主对角线上全是非零值的那个排序矩阵（唯一）$\tilde{W}$，实际运算中误差会导致W的所有值都是非零的，所以实际算法中寻找$\sum_{i}1/|\tilde{W}_{ii}|$取最小值时的排列矩阵。
3. 对$\tilde{W}$中的每一行，除以相应行上对角线元素的值，生成一个对角线上全是1的新矩阵$\tilde{W}^{\prime}$。
4. 利用如下公式计算矩阵B的一个估计值：$\hat{B}=I-\tilde{W}^{\prime}$
5. 最后，为了找到因果顺序，需要找到一个变换矩阵P（对行和列进行同样的变换），使得B在经过$\tilde{B}=P\hat{B}P^{T}$变换之后能够非常接近于一个下三角矩阵。（可以使用公式$\sum_{i<=j}{\tilde{B}_{i,j}^{2}}$来检验）。

• 第5步这里应该如何对B进行排序呢？如果算法精准，所有该等于0的地方都等于0，则可以直接采用如下算法（算法B1）：
  1）初始化一个存放排序的列表p。
  2）重复如下步骤直到B不再包含任何元素：
  （i）找到全是0的第i行，如果没有则返回False；
  （ii）Append i到列表p中；
  （iii）移除B的第i行和第i列。
  3） 返回True并输出排序列表p。
  但实际算法中不可能如此精准，该为0的地方不能够全都精确为0，则采用如下算法：
  1）将矩阵B中绝对值最小的m(m+1)/2个元素全都置为0。
  2） 循环如下步骤：
  （i）使用B1算法检验B是否能够被排序成一个严格的下三角矩阵，如果能，停止循环并返回B；
  （ii）把B中下一个绝对值最小的元素也置为0；

四、检验

1. 剪枝检验（检验一）

算法完成后，矩阵B中的上三角部分都被置零；但下三角部分基本全都是非0，即预测出的基本是一个全连接网络。剪枝即通过Wald检验，把B矩阵中不显著的系数置为0，从而简化网络连接。
检验方式：零假设（null hypothesis）：H0：bij = 0等价于H0：wij = 0
可以通过以下统计量进行检验：
$\tilde{W}^{2}/avar(\tilde{W}^{2})$
其中avar(*)代表渐进方差。这样的统计量渐进近似于一个自由度为1的卡方统计量。通过与经验值相比较，我们就可以对H0是否成立做出判断，即bij是否显著为0，如果不显著，则拒绝H0假设，即认为bij不显著为0，可以保留该系数。

2. 模型拟合度检验（检验二）

首先我们引入一些运算符号：

其中$vec^{+}(X)$表示这样一种运算，即将矩阵X的协方差矩阵中所有重复元素删除掉后，再纵向排列为一个向量，具体举例如下：

令x1,…xn为LiNGAM模型中随机采样的样本，定义：

当样本数n足够大的时候，m2就可以用来估计$\sigma_{2}(\tau_{0})$，通过衡量m2和$\sigma_{2}(\tau_{0})$之间的距离，可以判断模型拟合的好坏，距离越大，说明拟合的越差。

依然采用零假设的检验方式H0 : $E(m_{2}) = σ2(τ_{0})$，构造统计量：$T_{1}=n\times F(\hat{\tau})$，其中：

$T_{1}$统计量渐进近似于一个自由度为u-v的卡方变量，其中u为所使用的不同矩的数目，v为表示σ2(τ)的二阶矩所使用的参数的个数。$T_{1}$统计量就可以用来对H0进行检验。然而，T1需要大量的样本才能表现得像一个卡方变量，所以通常使用T2统计量：

3. 嵌套模型的差异卡方检验（检验三）

嵌套模型：假设模型1和模型2分别拥有q和q-1条边，即模型2是模型1剪掉一条边后的简化版，则模型1和2就被成为嵌套模型。
检验方式：用$T_{2}(q)$和$T_{2}(q-1)$分别表示模型1和模型2的拟合统计量，则可以通过检验$T_{2}(q)-T_{2}(q-1)$的显著性来表明模型1和模型2的拟合效果是否有显著差异，如果没有显著差异，则表明边可以被修剪掉。

4. 综合三种检验方法的剪枝检验算法

我们通过以下算法过程来对LiNGAM模型中的每个边的显著性进行一个系统的检验：

1. 设定一个显著性水平a（比如0.05）；
2. 对每条边进行Wald检验（检验一），找出所有不显著的边；
3. 把严格下三角矩阵B中最不显著的边所对应的系数元素置为0；
4. 重复以下步骤，直到所有不显著的边都被检验过：
   (a) 使用差异卡方检验方法（检验三），检验剪掉一条边后的模型和不剪这条边的模型之间拟合度是否有显著差异；使用卡方检验方法（检验二），检验当前模型的拟合效果是否显著；如果两个零检验都通过（即剪掉这条边不会带来显著效果差异，而且剪掉后模型依然有良好的拟合效果），则剪掉这个边，保留当前模型。
   (b) 进一步把矩阵B中下一个最不显著的边对应的元素置为0。
5. 返回剪枝后的矩阵B。

五、算法效果验证

1. 实验数据检验算法准确性

首先，我们通过人为构造出完全服从非高斯分布的数据，并检验该算法能否正确拟合出相应的系数。数据构造及检验过程如下：

1. 首先，随机构造严格的下三角矩阵B，涵盖不同维度（3，5，10，20，100），涵盖不同稀疏度（全连接，稀疏矩阵）；随机设定扰动项的方差和常数项ci。
2. 从一个高斯分布中对$e_{i}$进行独立采样，然后对$e_{i}$进行一个幂的非线性化（符号不变，绝对值进行一个区间在 [0.5, 0.8]或[1.2, 2.0]之间的指数变换）以使e非高斯化，最后调整e的尺度使e达到我们设定的方差。有了e之后，就可以构造出样本数据X。
3. 打乱样本数据X的行顺序。
4. 把数据X喂给上述算法，得到$\hat{B}$，比较人为构造的B和算法计算出来的$\hat{B}$之间的差异。
   验证的结果可用下图表示：
   
   其中横轴代表数据的样本量，纵轴代表变量数，可以看到数据维度在10维以下时，200个及以上的样本量就足够LiNGAM算法正确拟合出系数B的值，而当变量维度为20时，1000个变量才能让算法达到令人满意的精度，当维度变为100时，则至少需要5000个变量。总之，从构造的实验数据来看，算法对系数b的拟合准确度还是可以的。
   
   上面的表格是在实验数据上对剪枝效果的一个验证结果，显著性水平设置为5%，其中True Positive表示应保留的边被正确识别，False Negative表示应保留的边被错误剪掉，True Negative表示应剪掉的边被正确剪掉，False Positive表示应剪掉的边没有被剪掉。从表格结果可以看到，只要样本量充足，剪枝整体的正确率都在90%以上。

2. 真实数据检验算法适用性

本篇主要是针对几类不同的时序序列，检验了LiNGAM算法的适用性。首先我们需要了解几个概念和问题：

1）什么是p阶自回归过程？

一个严格意义上的p阶自回归过程的时间点数据都是通过如下公式构造出来的：

如果要用LiNGAM方程去模拟一个p阶自回归过程，则要求白噪声at必须满足非高斯分布。

2）如何把自回归时间序列数据构造成LiNGAM模型的训练数据？

固定时间窗口滑动采样。

3）固定时间窗口滑动采样对于模拟1阶以上的自回归过程会有什么问题？

潜变量（Confounder）问题，我们以一个2阶自回归过程为例，滑动时间窗口长度选为3，如下图所示：

对于一个时间窗口中的$X_{t+2}$而言，预测它的信息是充分的，但对于$X_{t+1}$而言，预测它的信息是不充分的，因为时间窗口的限制导致预测$X_{t+1}$所必须的$X_{t-1}$数据点并不在样本当中，所以系数$b_{21}$理论上讲是永远预测不准的。

4）如何得知一个真实的时间序列数据是几阶自回归过程呢？

借助acf（自回归函数）和pacf（偏自回归函数）分析，简单来讲acf衡量了$X_{t}$和$X_{t-k}$之间的关系，而pacf是在去除了$X_{t-1},X_{t-2},......X_{t-k+1}$的影响后，衡量$X_{t}$和$X_{t-k}$之间的关系。（具体可以看这篇博客）

了解了以上几个问题之后，我们首先来看一下典型的四类时序序列的特征及LiNGAM算法对它们的模型效果：

第一个时序序列为美国东港的月降水量序列，从其pacf分析中可以看出，滞后1阶和2阶的偏自相关系数是在95%显著性水平上是显著的，表明这样一个时序序列可以近似用一个2阶自回归过程来表述，除此之外，滞后11阶偏自相关系数也显著，表明其带有一定的年际振荡特征。LiNGAM算法能够准确计算出这样一个时序的正确因果关系.

第二个时序序列为IBM的股票数据，从pacf分析中可以看出，其只有滞后1阶偏自相关系数是显著的，表明这样一个时间序列基本就是一个随机游走的过程，LiNGAM算法找到了与之相反的因果关系。

第三个序列为惠拉维德利地区的日降水量，从pacf分析中看出它没有明显显著的偏自相关系数，LiNGAM算法对这样的时间序列进行估计，得出的系数矩阵B中的元素基本上都为0或者不显著的十分小的值。

第四个时间序列为夏威夷莫那罗亚地区的月而养活他排放数据。LiNGAM算法并没有正确拟合出它的时序因果关系，很有可能是LiNGAM的基本前提假设并不满足，包括：一，这个时间序列间的前后因果关系可能是非线性的；二，该时序的变化特征可能受某些外部自然环境（潜变量）所影响，并不能够完全通过时序本身推导出前后因果关系；三，时间样本并不满足服从非高斯分布的条件。

---

参考文章：

1. Shimizu S , Hoyer P O , Hyvrinen A , et al. A Linear Non-Gaussian Acyclic Model for Causal Discovery[J]. Journal of Machine Learning Research, 2006, 7(4):2003-2030.
2. http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=A71A4C9E6102B30C77E7C5CD40948C82?doi=10.1.1.105.5884&rep=rep1&type=pdf