如何基于伽马分布生成均值不变的负相关双随机变量？

原代码问题分析 #

你提供的Cholesky分解方法存在两个核心问题：

• 单向调整：由于Cholesky矩阵是上三角结构，右乘后第一列变量完全等于原始X，仅第二列Y被修改，不符合双向调整的需求；
• 边缘分布偏离：线性变换后的变量不再服从伽马分布，因为伽马随机变量的线性组合不满足伽马分布的特性。

同时，原方法无法保证变换后变量的均值与原始样本一致，不符合你最小化均值偏差平方和的要求。

解决方案一：Copula方法（保证边缘为伽马分布） #

Copula是生成指定相关系数的同边缘分布变量的标准方案，能在控制变量相关性的同时，严格保留边缘分布的类型。步骤如下：

1. 生成具有目标相关系数的标准正态变量；
2. 将正态变量转换为均匀分布；
3. 通过伽马分布的分位数函数，将均匀变量转换为符合要求的伽马变量。

R代码实现 #

# 样本量
n <- 10

# 伽马分布参数（shape，默认scale=1，均值=shape*scale）
shape1 <- 2
shape2 <- 3
target_corr <- -0.1

# 1. 生成指定相关系数的标准正态变量
corr_matrix <- matrix(c(1, target_corr, target_corr, 1), nrow = 2)
norm_vars <- MASS::mvrnorm(n, mu = c(0,0), Sigma = corr_matrix)

# 2. 转换为均匀分布变量
unif_vars <- pnorm(norm_vars)

# 3. 转换为伽马分布变量（保持原始均值）
correlated_X <- qgamma(unif_vars[,1], shape = shape1)
correlated_Y <- qgamma(unif_vars[,2], shape = shape2)

# 验证结果
cat("相关系数：", cor(correlated_X, correlated_Y), "\n")
cat("原始X均值：", mean(rgamma(n, shape=shape1)), " 变换后X均值：", mean(correlated_X), "\n")
cat("原始Y均值：", mean(rgamma(n, shape=shape2)), " 变换后Y均值：", mean(correlated_Y), "\n")

解决方案二：线性变换+均值约束（基于原始样本调整） #

如果你希望基于已生成的独立伽马样本X、Y进行双向调整，同时保证均值不变、相关系数符合要求，可以通过约束优化求解线性变换参数：

1. 设定线性变换形式，加入均值保持约束；
2. 优化变换参数，使相关系数达到目标值，同时最小化样本偏差平方和。

R代码实现 #

# 样本量
n <- 10

# 伽马分布参数
shape1 <- 2
shape2 <- 3
target_corr <- -0.1

# 生成原始独立伽马样本
set.seed(123) # 固定种子保证可复现
X <- rgamma(n, shape = shape1)
Y <- rgamma(n, shape = shape2)
mu_X <- mean(X)
mu_Y <- mean(Y)

# 定义目标函数：给定变换参数b和c，计算当前相关系数与目标的偏差平方
objective <- function(params) {
  b <- params[1]
  c <- params[2]
  # 均值约束下的变换系数
  a <- 1 - b * (mu_Y / mu_X)
  d <- 1 - c * (mu_X / mu_Y)
  
  # 计算变换后的变量
  X_prime <- a * X + b * Y
  Y_prime <- c * X + d * Y
  
  # 目标：最小化相关系数偏差 + 样本偏差平方和（增强与原始样本的相似性）
  (cor(X_prime, Y_prime) - target_corr)^2 + 1e-4 * (sum((X_prime - X)^2) + sum((Y_prime - Y)^2))
}

# 初始参数猜测
initial_guess <- c(0, 0)

# 优化求解
opt_result <- optim(initial_guess, objective)
b_opt <- opt_result$par[1]
c_opt <- opt_result$par[2]

# 计算最终变换系数
a_opt <- 1 - b_opt * (mu_Y / mu_X)
d_opt <- 1 - c_opt * (mu_X / mu_Y)

# 生成变换后的变量
X_prime <- a_opt * X + b_opt * Y
Y_prime <- c_opt * X + d_opt * Y

# 验证结果
cat("相关系数：", cor(X_prime, Y_prime), "\n")
cat("原始X均值：", mu_X, " 变换后X均值：", mean(X_prime), "\n")
cat("原始Y均值：", mu_Y, " 变换后Y均值：", mean(Y_prime), "\n")
cat("是否双向调整：", !all.equal(X_prime, X) & !all.equal(Y_prime, Y), "\n")

说明 #

• Copula方法是生成指定相关伽马变量的标准方案，严格保证边缘分布为伽马，适合大多数场景；
• 线性变换方法基于原始样本调整，满足双向修改的需求，但变换后的变量不再是伽马分布，仅适合对边缘分布要求不严格的场景。

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