使用Copula模拟自定义分布相关数据时相关性不符的问题排查

问题背景 #

我使用Copula三步法模拟相关数据，三步法为：

• 从相关矩阵模拟多元正态相关数据，边缘为标准正态分布；
• 用标准正态CDF将正态边缘转换为均匀分布；
• 用逆CDF将均匀边缘转换为目标分布。

已成功用R复现示例代码，代码如下：

library(MASS)
library(VGAM)  

set.seed(12345)
N <- 1000
# 定义相关矩阵Sigma
Sigma <- matrix(c(1.0, 0.5, 0.5, 1.0), nrow = 2)
# 1. 生成相关多元正态数据Z
Z <- mvrnorm(n = N, mu = c(0, 0), Sigma = Sigma)
# 2. 将边缘变量转换为U(0,1)分布
U <- pnorm(Z)
# 3. 构造目标边缘分布
gamma <- qgamma(U[, 1], shape = 4)
LN <- qlnorm(U[, 2], meanlog = 0.5, sdlog = 0.8)
X <- cbind(gamma, LN)

# 计算秩相关系数
rankCorrZ <- cor(Z, method = "spearman")[2]
rankCorrU <- cor(U, method = "spearman")[2]
rankCorrX <- cor(X, method = "spearman")[2]

# 计算Pearson相关系数
rhoZ <- cor(Z, method = "pearson")[2]
rhoX <- cor(X, method = "pearson")[2]

输出结果显示秩相关性保持一致：

> print(rankCorrZ)
#[1] 0.471292
> print(rankCorrU)
#[1] 0.471292
> print(rankCorrX)
#[1] 0.471292
> print(rhoZ)
#[1] 0.4930677
> print(rhoX)
#[1] 0.4636544

我希望生成符合以下自定义分布的x和y：

sigma_x <- 2
sigma_y <- 2
beta = c(4,1)
x <- rnorm(N, 0, sigma_x)
y <- beta[1] + beta[2]*x + sigma_y * rexp(N, rate = 1) - 1

沿用三步法实现后，得到的相关性远高于Z和U的结果，测试代码如下：

Sigma <- matrix(c(1.0, 0.5, 0.5, 1.0), nrow = 2)
Z <- rmvnorm(n = N, mean = c(0, 0), sigma = Sigma)
U <- pnorm(Z)
x <- qnorm(U[,1], mean = 0, sd = sigma_x)
y <- beta[1] + beta[2] * x + sigma_y * (qexp(U[,2], rate = 1) - 1)
DATA <- cbind(x, y)

rankCorrZ <- cor(Z, method = "spearman")[2]
rankCorrU <- cor(U, method = "spearman")[2]
rankCorrDATA <- cor(DATA, method = "spearman")[2]
rhoZ <- cor(Z, method = "pearson")[2]
rhoDATA <- cor(DATA, method = "pearson")[2]

输出结果：

> print(rankCorrZ)
#[1] 0.4609489
> print(rankCorrU)
#[1] 0.4609489
> print(rankCorrDATA)
#[1] 0.8821912
> print(rhoZ)
#[1] 0.4830919
> print(rhoDATA)
#[1] 0.8609049

请问为何会出现这种情况？我哪里操作错误？

错误原因分析 #

你的核心问题是混淆了Copula的适用场景：
Copula三步法的作用是生成边缘分布独立、但联合分布具有指定相关性的变量，每一个变量的边缘分布都是独立通过逆CDF从均匀变量转换而来的。

但你在构造y时，直接让y依赖于x（y <- beta[1] + beta[2] * x + ...），这相当于给y和x添加了两层相关性：

1. Copula本身通过U1和U2的相关性带来的x与y的关联；
2. 你手动添加的线性依赖关系beta[2]*x带来的强线性关联。

这两层相关性叠加后，最终的x和y的相关性自然会远高于Copula预设的0.5相关系数。

另外，你想要模拟的目标分布中，y是x的条件分布（y|x ~ 移位指数分布），而不是独立的边缘分布。这种场景下，根本不需要用Copula来生成——直接先生成x，再基于x生成对应的y即可，这才是符合目标分布定义的做法：

set.seed(12345)
N <- 1000
sigma_x <- 2
sigma_y <- 2
beta = c(4,1)
x <- rnorm(N, 0, sigma_x)
# 直接基于x生成条件分布的y
y <- beta[1] + beta[2]*x + sigma_y * rexp(N, rate = 1) - 1
DATA <- cbind(x, y)

# 查看相关性
cor(DATA, method = "spearman")[2]
cor(DATA, method = "pearson")[2]

如果你的真实需求是用Copula生成具有指定相关性、且边缘分布分别为正态和“移位指数混合”的变量，那你需要先推导y的边缘分布（对x积分后的边际分布），然后用该边缘分布的逆CDF去转换U2，而不是让y直接依赖于x。

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