基于泰勒模式自动微分的阿基米德Copula推断

# 通过泰勒模式自动微分的阿基米德 Copula 推断

来源：https://arxiv.org/html/2605.23134

**杨剑桥** 独立研究员，[email protected] & **李冬冬** 哈佛医学院，[email protected]

###### 摘要

目前没有任何现成的嵌套阿基米德 copula 工具能够同时处理以下三个问题：(a) 生存分析中每个变量任意的右删失，(b) 任意的嵌套树结构，以及 (c) 精确的参数梯度。现有的实现只能处理二元问题、低维（即 d≤10）情形、两层嵌套，或仅限于手工推导的 copula 嵌套。我们提出 `acopula`，一个 JAX 原生框架，给定任何阿基米德生成元——经典的或神经网络的——它都能在多项式时间内计算任意删失掩码下的精确嵌套 copula 似然和参数梯度。其机制是对泰勒模式自动微分输出进行多项式幂运算，用一个可微分的计算替代了每个族手工推导的部分 Bell 多项式表，任何用户定义的生成元都可以驱动该计算。我们进行了广泛的模拟验证 `acopula` 的正确性。然后我们展示了：(a) 在 MIMIC-IV 的 85,229 条 ICU 住院记录上进行高维（d=53）的每个变量删失分析，分别使用经典阿基米德族和嵌套神经网络阿基米德 copula 进行拟合；(b) 在 S&P 500 日收益率上的一个 11 部门层次模型（d=98）；(c) 在一项视网膜病变研究中，跨越十个族（其中五个族之前没有实现）的与族无关的删失 MLE 方法；以及 (d) 在 d=35 时，相比 R 语言的 `nacLL` 实现了约 650 倍的密度计算加速，并且在 d=8,000 时仍保持二次缩放。

## 1 引言

在现代生存分析、金融和临床风险环境中，部分观测的分层多元数据是常态。一项关于配对器官的生存研究中的患者可能在随访期间从未经历事件（Huster 等人，1989）。股票收益率在行业板块内聚类，并在危机期间产生尾部联动。重症监护病房（ICU）的实验室测量值按临床面板分组，分析人员独立地对每个实验室指标进行出院右删失（Johnson 等人，2023）。每种情景都需要一个模型，该模型能够以可解释的参数——尾部依赖性和一致性——捕捉层次结构，并容纳每个观测值上的任意删失模式。

##### 阿基米德 copula。

一个 copula C: [0,1]^d → [0,1] 是单位超立方体上的一个多元分布，具有 d 个均匀边缘分布；根据 Sklar 定理（Sklar, 1959），每个关于 u = (u_1, ..., u_d) ∈ [0,1]^d 的联合分布都可以唯一地分解为其一维边缘分布和一个捕捉它们依赖关系的 copula。这种分解允许分析人员将特定领域的边缘模型——例如，用于生存时间的 Weibull 分布（Weibull, 1951），用于收益率的 GARCH 模型（Bollerslev, 1986）——与一个可解释的依赖结构配对，该结构的参数直接编码尾部行为和相关一致性。

嵌套阿基米德 copula（McNeil, 2008；Hofert 和 Pham, 2013）通过在树中组合生成元函数将 copula 扩展到层次依赖：同一内部节点下的变量比跨分支的变量具有更强的依赖关系。

##### 推断障碍。

拟合一个嵌套 copula 需要其密度——copula 累积分布函数 (CDF) 的混合偏导数，每个观测变量一个微分——通过一个由组合生成元构成的树（Hofert 和 Pham, 2013）；三个障碍阻碍了现有工具。

*   **每个变量的删失。** 我们考虑生存设置下的每个变量右删失。每个被删失变量的子集会改变定义似然的 copula CDF 的混合偏导数，因此一个包含数千种不同删失模式的 d=53 ICU 数据集需要框架能够统一处理每个观测值的掩码 δ ∈ {0,1}^d，其中 δ_j=1 表示变量 j 被观测。现有工具涵盖二元问题（Sun 和 Ding, 2020；Emura 和 Chen, 2018）或共享删失（Liu 等人, 2025b）；没有一个能达到对每个变量进行 d>30 的删失。
*   **任意嵌套和任意生成元。** R 语言的 `copula` 包（Hofert 等人, 2012）硬编码了五个族的 Stirling 数表，并且只支持两层嵌套；添加第六个族——或任何基于神经网络的生成元，例如深度阿基米德 copula 族（Ling 等人, 2020；Liu 等人, 2025b）——需要重新进行手工推导。
*   **精确的自动参数梯度。** 没有自动梯度，拟合一个 d=98 的模型就要依赖缩放性差的无导数或有限差分方法，或者蒙特卡洛近似器（Hofert 等人, 2013；Ng 等人, 2021），这些方法虽然一致但在有限样本下有偏。

##### 我们的方法。

我们提出 `acopula`，一个 JAX 原生框架，解决了所有三个障碍。用户至少指定一个 copula 生成元——经典的或神经网络的——然后 `acopula` 为任何树形状和任何删失掩码计算嵌套 copula 似然，并自动提供精确的参数梯度。其机制是泰勒系数的多项式幂运算：泰勒模式自动微分（Griewank 和 Walther, 2008；Radul 等人, 2023）在一次前向传播中返回每个生成元组合的高阶导数，通过卷积对这些系数向量进行连续幂运算，恢复了 Hofert 和 Pham (2013) 密度公式所需的部分 Bell 多项式，无需为每个族进行手工推导，适用于任何满足标准 d_c-单调嵌套条件（McNeil, 2008）的生成元。

##### 贡献。

我们做出如下贡献：

1.  **一种与族无关的、基于梯度的嵌套阿基米德 copula 似然推断算法**：只需给出一个经典或神经网络生成元（例如，深度阿基米德 copula，Ling 等人, 2020；Liu 等人, 2025b），该框架自动推导密度和梯度。
2.  **将每个变量删失作为一等操作**：似然算法通过其核心递推直接处理任意观测掩码 δ ∈ {0,1}^d，突破了先前删失感知工具仅限二元或单一共享时间的限制。
3.  **一个复杂度结果**：对于固定深度树，最坏情况为 O(d^3)；对于 √d-分解的固定深度树为 O(d^2)，其中每个非根子树最多持有 O(√d) 个未删失叶子；经验证 d 可达 8,000。
4.  **一个开源 JAX 实现 `acopula`**，提供简单的 API 用于嵌套 copula 密度和梯度，以及用于嵌套有效性条件的每个边 d_c-单调性诊断。
5.  **作为副产品**，这些 copula 的条件分布抽样（Hofert, 2008）变得可行，此前被高阶导数所阻碍；由于篇幅限制，我们将其推迟到附录 F。

##### 相关工作。

我们将 `acopula` 与四个相关的工作方向进行对比，每个方向总结如下；表 1 以表格形式汇总了由此产生的能力差距。

**表 1：** 嵌套 copula 和 copula-生存工具的能力比较（R 语言 `copula`（Hofert 等人, 2012, 2013），`HAC`（Okhrin 和 Ristig, 2014），`Copulas.jl`（Laverny 和 Jimenez, 2024），`HACSurv†`（Liu 等人, 2025b），`DCSurvival`（Foomani 等人, 2023），`CopulaCenR`（Sun 和 Ding, 2020），以及藤结构实现 `rvinecopulib`（Nagler 和 Vatter, 2023））。测试的 d 报告了每个工具在嵌套真实数据示例中成功尝试的最大值，包括本工作。

| 功能 | R copula | HAC | Copulas.jl | HACSurv† | DCSurvival | CopulaCenR | Vine (rvinecopulib) | 我们的工作 |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| 用户定义经典生成元 | × | × | × | ✓ | N/A | × | × | N/A | ✓ |
| 深度阿基米德生成元 | × | × | × | × | ✓ | ✓ | × | N/A | ✓ |
| 抽样 | ✓ | ✓ | ✓ | N/A | ✓ | ✓ | ✓ | ✓ |
| 嵌套深度 (精确密度) | 2 | 任意 | N/A | 2 | N/A | N/A | N/A | 任意 |
| 自动密度梯度 | × | × | × | ✓ | ✓ | × | × | ✓ |
| 与族无关的密度 | × | × | × | × | ✓ | N/A | × | ✓ |
| 每个变量删失 | 无 | 无 | 无 | 竞争风险 | 二元 | 二元 | 无 | 每个变量 |
| 测试的 d (嵌套真实数据) | 35 | 35 | ∼10 | N/A | 6 | 6 | 22 | 98 |

*   **经典 R 包生态系统。** R 语言的 `copula` 包（Hofert 等人, 2012, 2013）和 `HAC` 包（Okhrin 和 Ristig, 2014）通过五个族的硬编码闭合形式计算生成元导数，将每条边限制为同族对，并且不提供解析参数梯度。模拟最大似然估计器（Hofert 等人, 2013）依赖于拉普拉斯变换脆弱性抽样，其蒙特卡洛方差随维数增长；我们的精确似然方法在匹配精度下实现了约 100 倍的加速（附录 Q），并且我们对 R 的 `nacLL` 的确定性密度比较在 d=35 时达到约 650 倍，随后 R 中止（第 5.2 节）。
*   **感知删失的 copula 模型。** `CopulaCenR`（Sun 和 Ding, 2020）拟合删失的二元和三元 copula，包括通过生存 copula 的 2^d 角包含-排除法进行区间删失；Li 等人（2020, 2023）使用暴力符号微分处理右删失。两者都局限于低环境维数。`HACSurv`（Liu 等人, 2025b）拟合用于竞争风险的层次阿基米德 copula，其中每个观测值缩减为一个事件时间，似然是嵌套 CDF 的一个单一的一阶偏导数。我们处理生存意义下每个变量右删失产生的高阶混合偏导数，并在第 6 节讨论区间删失。
*   **深度和概率性 ML copula。** 最近的神经网络 copula 工作（Foomani 等人, 2023；Liu 等人, 2025a）学习用于二元依赖删失和离散扩散联合抽样的灵活生成元；我们是补充性的，目标是在中等至高维 d 下使用每个变量删失进行精确嵌套似然计算。藤 copula（Aas 等人, 2009；Nagler 和 Vatter, 2023；Stöber 等人, 2013）和层次 Kendall copula（Brechmann, 2014）是密度优先的：删失似然需要对删失块进行每个观测值的多元积分，计算量随其大小指数增长。METIC（Chen 等人, 2025）通过随机数值方案在小的 d 下处理此问题；没有藤结构实现在 d=53 时达到每个变量删失（附录 U 和 V）。
*   **高阶自动微分和可微分概率系统。** 我们的截断泰勒级数计算使用了 JAX 的泰勒模式自动微分原语（Bettencourt 等人, 2019），也在 Betancourt (2018) 的 Weil 代数框架中形式化了：在截断多项式代数上的单次前向传播中计算 d 阶导数（第 3.2 节）。Lin (2024) 将 JVP-规则协议扩展到泛函导数，与我们在此使用的泰勒模式 AD 互补。NumPyro（Phan 等人, 2019）公开了一阶梯度，但不是嵌套 copula 所需的高阶梯度；Gen-AC（Ng 等人, 2021）学习拉普拉斯变换，其蒙特卡洛方差随导数阶数增长。

##### 路线图。

第 2、3 和 4 节公式化问题，提出关键算法，并描述 JAX 实现。第 5 节通过模拟验证正确性，描述相对于 R 的 `nacLL` 的运行时间缩放（在 d=35 时达到约 650 倍加速，随后 R 中止），并进行四个真实数据实验：在 d=53 的 MIMIC-IV 上使用经典和神经网络阿基米德生成元进行每个变量删失；在 d=98 的 S&P 500 上——这是 R 无法企及的嵌套似然；跨越十个阿基米德族的糖尿病视网膜病变分析，其中五个族没有任何先前的删失实现；以及 HACSurv 的 Framingham 竞争风险管道的端到端重现。第 6 节讨论了混合族的有效性和跨观测值的摊销 JIT 编译成本。

## 2 背景与符号

![图1：六个叶子上的嵌套阿基米德 copula。节点携带生成元，形成其子节点的 copula CDF；内部子树捕捉组内依赖关系。叶子是具有边缘密度的变量；阴影叶子 u3 被删失。](https://arxiv.org/html/2605.23134/fig_arxiv/dag.png)

**图 1：** 六个叶子上的嵌套阿基米德 copula。节点携带生成元，形成其子节点的 copula CDF；内部子树捕捉组内依赖关系。叶子是具有边缘密度的变量；阴影叶子 u3 被删失。

一个 copula C: [0,1]^d → [0,1] 捕捉 d 个均匀随机变量之间的依赖关系，独立于它们的边缘分布（Sklar, 1959; Embrechts 等人, 2002; Nelsen, 2006）。一个**阿基米德 copula** 通过单一的生成元 ψ: [0,∞) → [0,1] 来编码这种依赖关系——一个连续的、严格递减且满足 ψ(0)=1 的函数，它是 d-单调的，对于 d=2 等价于凸函数，对于任意 d 则通过完全单调性保证（McNeil 和 Nešlehová, 2009），产生