LiNGAM与非高斯：线性非高斯模型LiNGAM

可视化演示

动态展示LiNGAM模型如何利用非高斯性识别因果方向：源变量混合产生高斯化观测，ICA解混恢复非高斯源，进而确定唯一因果顺序与有向边。

模型设定：一张机器的装配图

老师傅修机器时，不会凭空猜测，而是先画出零件之间的传动关系。如果齿轮 x1 带动了齿轮 x2，他就写下“x2 的转动量 = 0.7 倍的 x1 转动量 + x2 自己的咔哒扰动”。这句话有两层意思：一是 x2 确实由 x1 影响，二是 x2 身上还有一股不属于 x1 的“原生活力”。LiNGAM 的模型设定正是如此——每个变量都等于它所有父节点的线性组合，再加上一个独立的噪声项。

用数学写出来，对第 i 个变量有：

其中 \(b_{ij}\neq 0\) 表示存在从 \(x_j\) 到 \(x_i\) 的因果边；所有 \(b_{ij}\) 拼在一起就是连接强度矩阵 \(B\)。同时要求这些变量之间有向无环，也就是不能出现“x1 带 x2、x2 带 x3、x3 又带回 x1”这种卡死的循环。这个 DAG 假设保证了因果方向能排出先后顺序。

回到车间视角：只要记录下每个轴转了多少、机器内部怎么响，老师傅就能用这张表反推传动关系。接下来最关键的问题，就是为什么“咔哒声”必须非高斯。

非高斯核心：齿轮的独特声音

如果所有齿轮都发出同一种均匀、对称、钟形的嗡嗡声，老师傅会非常头疼。因为多个这样的声音叠加后还是相同的嗡鸣，就像高斯噪声叠加后仍服从同一个分布，谁也没法分辨声源。相反，如果齿轮发出的“咔哒”是有边界的均匀抖动，或者是斜向一侧的指数式短促冲击，那么这些声音即使混在一起，也有各自鲜明的“签名”。

LiNGAM 要求每个外生项 \(e_i\) 连续、相互独立，并且是非高斯分布。这样，观测变量的联合概率密度就可以写成每个噪声密度的乘积，再除以一个由 \(B\) 决定的伸缩因子：

这个公式的本质是：先从观测 \(x\) 中解出每个“咔哒” \(e_i\)，然后利用独立性把联合概率拆成多个单零件概率的乘积。如果噪声中有两个以上高斯源，旋转后的分布仍然像同一批声音，上面的分解就会出现多解；最多只能有一个高斯源，才能保证声音签名彼此可分。

于是，独特音色让每个齿轮在混合声中保持可分辨。下一步，我们要把这种直觉变成严格的“唯一性”。

可识别性：为什么装配图只有一张

同一台机器，从外面看也许有几种内部装配都能产生相同的转动表现。如果只有高斯嗡嗡声，这些装配方案根本无法区分。但有了非高斯咔哒声，情况就不同了：只有真实的传动顺序，才能让每个齿轮的噪声源在分离后真正独立；其他顺序都会留下“残留的相互牵扯”。

LiNGAM 的核心可识别性定理说：在线性、DAG、外生项独立非高斯（最多一个高斯）的假设下，观测分布能唯一确定因果顺序和连接矩阵 \(B\)。把模型改写成：

这里 \(A\) 是“混合矩阵”，表示每个噪声源如何混合成观测变量。非高斯保证 \(A\) 可以被唯一识别；如果是高斯噪声，\(A\) 可以任意乘一个正交旋转矩阵 \(Q\)，得到同样的观测分布，因果方向就说不清了。

唯一的 \(A\) 意味着能还原唯一的 \(B\)。那么，具体怎么把独立声源从混合观测中提出来？这就轮到 ICA 登场。

ICA方法：用盲源分离听声辨位

老师傅在机器外壳上放了几只麦克风，每只麦克风录到的是多个齿轮声音的混合。要从混合录音里把每个齿轮的咔哒声分别拎出来，需要一种“盲源分离”技术。ICA（独立成分分析）做的就是这件事：给定观测到的混合信号，在源信号非高斯且独立的假设下，估计出混合矩阵和独立成分。

在 LiNGAM 中，观测 \(x = A e\)，所以可以直接对 \(x\) 跑 ICA，估计 \(\hat A\) 与独立成分 \(\hat e\)。令 \(W = \hat A^{-1}\)，则可以恢复连接矩阵：

这里 \(P\) 是置换矩阵，\(L\) 是严格下三角矩阵。为什么需要 \(P\)？因为 ICA 的输出有“先按哪个变量排序”的模糊性，而 DAG 要求因果顺序对应一个拓扑排序，恰好可以用下三角结构来消除这种排列模糊。

不过，借助 ICA 只是其中一条路线。另一位老师傅更喜欢更直接的听诊法，这就是 DirectLiNGAM。

DirectLiNGAM：逐步剥离法

老师傅戴上听诊器，不去做完整的盲源分离，而是先找到最上游的那个齿轮。最上游的齿轮有个特点：它的咔哒声不会被其他齿轮影响，因此当它去解释任何一个其他零件时，解释完剩下的残差和它自己是独立的。反过来，如果一个齿轮下游还有别的零件，那么用下游零件去回归它，残差会带着它的痕迹。

DirectLiNGAM 的迭代思路就是反复找最外生变量：对所有剩余变量计算“我回归你、你回归我”的残差独立性，把独立性最强的那位拎出来作为当前最上游变量；然后用回归剥离它对其他变量的影响，再对剩下的子机器重复。每一步的回归斜率是：

这个公式计算的是：把 \(x_i\) 中能被 \(x_j\) 解释的线性部分去掉，剩下 \(r_i^{(j)}\)。如果 \(x_j\) 是最外生的原因，那么 \(r_i^{(j)}\) 与 \(x_j\) 独立；如果 \(x_j\) 其实是结果，残差中会残留依赖。

剥完所有层之后，因果顺序已经确定。接下来要理解，这个顺序如何体现在最终的 DAG 和矩阵结构里。

因果顺序：流水线上的下三角

一旦确定谁在上游、谁在下游，整张装配图就成了一条不循环的流水线。数学上，DAG 一定存在拓扑排序：把变量重新编号后，所有因果箭头只从编号小的指向编号大的，连接矩阵 \(B\) 就变成严格下三角矩阵。

设因果顺序为 \(k_1, k_2, \dots, k_p\)，则重新排序后的矩阵满足：

这个条件非常漂亮：它同时表达了“有因果边，则原因一定排在结果之前”和“对角线为零，所以没有变量是自己影响自己”。没有环，机器才不会卡住。

现在我们已经有了完整的“装配图”。它和传统方法最大的区别，到底是什么？

方法对比：约束法只知相关不知方向

过去，很多因果发现方法是“约束法”：只看条件独立性，判断哪些变量之间有边。这就像只看两个齿轮的振动是否相关，能感觉到它们连在一起，却看不出谁带动谁。对于一条链 \(x_1 \to x_2\) 和反向 \(x_2 \to x_1\)，它们在条件独立性上完全等价，约束法只能把它们归入同一个马尔可夫等价类：

PC 等算法最后只能输出一个 CPDAG（部分有向图），其中某些边保持双向/无向。LiNGAM 则不同：它利用非高斯噪声的独立性分解，直接从观测分布中唯一挑选出真实 DAG。在满足假设时，它能确定每一条边的方向。

然而，这套听诊法并非万能。在真实车间的嘈杂环境中，它也有严格的适用边界。

应用与局限：在真实车间里

LiNGAM 已经在量化金融、营销归因、基因调控等领域得到应用。比如分析几支基金的日收益谁先驱动谁，或者用户从曝光到点击再到购买的传导路径。它最大的吸引力是：不需要做昂贵的随机实验，只用观测数据就能给出完整的 DAG。

但必须清醒地看到局限：LiNGAM 假设关系是线性的、因果图是无环的、没有隐藏变量、外生项独立且非高斯。如果实际机制是非线性或存在反馈，模型就会给出误导性的“最像因果”答案。样本量也很重要，独立性检验需要足够样本才能稳定。常用统计量近似为：

这个式子提醒我们：残差独立性是用样本估计的，样本越大，微弱依赖越容易被检测出来，因果判断也越可靠。

AI 应用场景

Python 演示

import numpy as np  def dvar(u):    n = len(u)    a = np.abs(u[:, None] - u[None, :])    A = a - a.mean(0) - a.mean(1)[:, None] + a.mean()    return np.sqrt(np.maximum((A * A).sum() / n**2, 0.0))  def dcor(u, v):    n = len(u)    a = np.abs(u[:, None] - u[None, :])    b = np.abs(v[:, None] - v[None, :])    A = a - a.mean(0) - a.mean(1)[:, None] + a.mean()    B = b - b.mean(0) - b.mean(1)[:, None] + b.mean()    dcov = np.sqrt(np.maximum((A * B).sum() / n**2, 0.0))    denom = np.sqrt(dvar(u) * dvar(v))    return 0.0 if denom == 0 else dcov / denom  def direct_lingam(X):    X = np.asarray(X, float)    p = X.shape[1]    order = []    active = np.arange(p)    cur = X.copy()     while len(active) > 1:        scores = []        for j in range(len(active)):            dep = 0.0            for k in range(len(active)):                if k == j:                    continue                xj = cur[:, j]                xk = cur[:, k]                beta = xj.dot(xk) / xj.dot(xj)                resid = xk - beta * xj                dep = max(dep, dcor(xj, resid))            scores.append(dep)         idx = int(np.argmin(scores))        order.append(active[idx])         xj = cur[:, idx]        cols = np.delete(np.arange(len(active)), idx)        next_cur = cur[:, cols].copy()        for t in range(next_cur.shape[1]):            beta = xj.dot(next_cur[:, t]) / xj.dot(xj)            next_cur[:, t] -= beta * xj         cur = next_cur        active = active[cols]     order.append(active[0])     order_map = {v: i for i, v in enumerate(order)}    B = np.zeros((p, p))    for j in range(p):        for i in range(j + 1, p):            cause = order[j]            effect = order[i]            xc = X[:, cause]            xe = X[:, effect]            B[effect, cause] = xc.dot(xe) / xc.dot(xc)     return np.array(order), B  def generate_data(n=800):    rng = np.random.default_rng(0)    e1 = rng.uniform(-1, 1, n)    e2 = rng.exponential(1, n)    e3 = rng.uniform(-1, 1, n)    x1 = e1    x2 = 0.8 * x1 + e2    x3 = 0.4 * x1 + 0.6 * x2 + e3    return np.column_stack([x1, x2, x3])  if __name__ == "__main__":    X = generate_data()    order, B = direct_lingam(X)    print("因果顺序:", order)    print("连接矩阵 B:\n", B)

1. dep = max(dep, dcor(xj, resid))计算当前变量与所有回归残差的最大依赖度；依赖越小，越可能是外生变量。
2. idx = int(np.argmin(scores))选出依赖度最小的变量作为本轮最外生变量，并追加到因果顺序末尾。
3. next_cur[:, t] -= beta * xj用回归剥离该变量对剩余变量的线性影响，缩小“机器”规模。
4. B[effect, cause] = xc.dot(xe) / xc.dot(xc)确定因果顺序后，对每个下游变量回归所有上游变量，得到连接强度矩阵 B 的非对角元。

常见误区

面试题