在R中使用Copula计算条件分布F(Y|X)与F(X|Y)

已知三个随机变量X、Z、Y的边缘分布函数F(X)、F(Z)、F(Y)，但变量间依赖结构未知，希望计算条件分布F(Y|X)和F(X|Y)。目前已使用R语言的copula和gamlss.dist包构建了theta=2的Gumbel copula三维联合分布，代码如下：

pacman::p_load(copula, gamlss.dist)

# Generate Gumbel copula with theta=2
thetaVal <- 2
copula <- gumbelCopula(thetaVal, dim = 3)

copPre <- mvdc(copula, c("SHASH", "norm", "norm"),
               list(
                 list(mu = 46, sigma = 5, nu = 2, tau = 1),  # X SHASH
                 list(mean = 46, sd = 1),                      # Y normal
                 list(mean = 23, sd = 1)))                     # Z normal

# 计算(X ≤ 47, Y ≤ 45, Z ≤ 22)的联合概率
pMvdc(c(47, 45, 22), copPre)

现寻求基于该copula模型计算目标条件分布的方法。

解决方案 #

核心原理 #

对于连续随机变量，条件分布F(Y|X=x)表示给定X=x时Y≤y的概率。利用copula的性质：

1. 三维联合分布可分解为边缘分布+copula：(F(X,Y,Z) = C(F_X(X), F_Y(Y), F_Z(Z)))
2. 边缘化Z后，(X,Y)的二维联合copula为原三维copula在(F_Z(∞)=1)时的结果，对于Gumbel copula，二维边缘copula与三维copula参数一致（均为theta=2）
3. 二维copula的条件分布公式：(F(Y≤y | X=x) = \frac{\partial C_{X,Y}(F_X(x), F_Y(y))}{\partial u})（其中(u=F_X(x))，(v=F_Y(y))），该偏导数可通过copula包直接计算。

代码实现 #

1. 加载依赖包

pacman::p_load(copula, gamlss.dist)

2. 定义基础模型（复用已有的三维copula定义）

thetaVal <- 2
# 三维Gumbel copula
copula <- gumbelCopula(thetaVal, dim = 3)
copPre <- mvdc(copula, c("SHASH", "norm", "norm"),
               list(
                 list(mu = 46, sigma = 5, nu = 2, tau = 1),  # X的SHASH分布参数
                 list(mean = 46, sd = 1),                      # Y的正态分布参数
                 list(mean = 23, sd = 1)))                     # Z的正态分布参数

# 构造(X,Y)的二维边缘Gumbel copula
bicop_xy <- gumbelCopula(thetaVal, dim = 2)

3. 计算F(Y ≤ y | X = x)

定义函数实现条件分布计算：

# 给定X=x，计算Y≤y的条件概率
cond_y_given_x <- function(y, x) {
  # 计算X=x对应的copula变量u=F_X(x)
  u <- pshash(x, mu = 46, sigma = 5, nu = 2, tau = 1)
  # 计算Y=y对应的copula变量v=F_Y(y)
  v <- pnorm(y, mean = 46, sd = 1)
  # 利用copula包的pcondCopula计算条件分布
  pcondCopula(v, bicop_xy, u, side = "left")
}

# 示例：计算X=47时，Y≤45的条件概率
cond_y_given_x(y = 45, x = 47)

4. 计算F(X ≤ x | Y = y)

同理，定义对应的条件分布函数：

# 给定Y=y，计算X≤x的条件概率
cond_x_given_y <- function(x, y) {
  # 计算Y=y对应的copula变量v=F_Y(y)
  v <- pnorm(y, mean = 46, sd = 1)
  # 计算X=x对应的copula变量u=F_X(x)
  u <- pshash(x, mu = 46, sigma = 5, nu = 2, tau = 1)
  # 计算条件分布
  pcondCopula(u, bicop_xy, v, side = "right")
}

# 示例：计算Y=45时，X≤47的条件概率
cond_x_given_y(x = 47, y = 45)

验证说明 #

若需要验证条件概率的合理性，可通过三维联合分布计算近似条件概率（给定X∈[x-ε, x+ε]）：

epsilon <- 0.01
# 计算P(Y≤45, X∈[47-ε,47+ε])
joint_prob <- pMvdc(c(47+epsilon, 45, Inf), copPre) - pMvdc(c(47-epsilon, 45, Inf), copPre)
# 计算P(X∈[47-ε,47+ε])
marginal_prob <- pMvdc(c(47+epsilon, Inf, Inf), copPre) - pMvdc(c(47-epsilon, Inf, Inf), copPre)
# 近似条件概率
approx_cond_prob <- joint_prob / marginal_prob

该结果会趋近于cond_y_given_x(45,47)的计算值（ε越小越接近）。

内容的提问来源于stack exchange，提问作者user17326436