生成具有受控相关性的 x 和 y

Question

我想生成具有受控相关性的 x 和 y 变量。目前，我正在使用以下代码来生成 x 和 y（如下）。但是，我想创建一个函数，在模拟 x 和 y 时可以指定 x 和 y 之间的相关性。有谁知道如何调整此代码以实现 x 和 y 之间的特定相关性？误差项 erorr_x 和 erorr_y 均值为零并且是独立的。

n = 100
sigma_x = 3
sigma_y = 3
beta = c(5,2)

SIM_DATA <- function(n=100, sigma_x, sigma_y, beta, rho, MRep=50){
  df <- list(NULL)
  for(i in 1:MRep){
    sim_data <- matrix(NA, nrow = n, ncol = 2)
    erorr_x <- rnorm(n, mean = 0, sd = 1)
    erorr_y <- rexp(n, rate = 1) - 1
    x <- sigma_x * erorr_x 
    y <- beta[1] + beta[2] * x + sigma_y * erorr_y 
    sim_data[, 1] <- x
    sim_data[, 2] <- y
    df[[length(df) + 1]] <- sim_data
  return(df)
}

Answer 1

下面，

fcor

是一个接受任何函数

f

的函数，它对二元分布进行采样，并使用随机求根来返回与

f

具有相同边际分布的样本，但具有指定的相关性

rho

。

library(mvtnorm)

fcor <- function(n, m, f, rho, k = 1e3, reps = 1e3) {
  # n = number of observations per sample
  # m = number of samples
  # f = function that generates random observations from a bivariate distribution
  # rho = target correlation
  # k, reps: parameters for the stochastic root finding process
  r <- rho
  mu <- c(0, 0)
  num <- den <- 0

  for (i in seq_len(reps)) {
    z <- rmvnorm(k, mu, r - diag(r - 1, 2))
    xy <- f(k)
    num <- num + r^2
    den <- den + r*cor(sort(xy[,1])[rank(z[,1])], sort(xy[,2])[rank(z[,2])])
    r <- rho*num/den
  }
  
  sigma <- r - diag(r - 1, 2)
  replicate(m, {
    z <- rmvnorm(n, mu, sigma)
    xy <- f(n)
    cbind(sort(xy[,1])[rank(z[,1])], sort(xy[,2])[rank(z[,2])])
  }, FALSE)
}

观察到

xy

来自

f

，并且每列中的值只是简单地重新排序以获得所需的相关性。因此，边际分布得以保留。

接下来，定义采样函数：

sigma_x <- sigma_y <- 3
beta = c(5, 3)

fgen <- function(n) {
  cbind(
    x <- rnorm(n, 0, sigma_x),
    beta[1] + beta[2]*x + sigma_y*rexp(n, rate = 1) - 1
  )
}

测试：

# generate 50 sets of 100 observations with an underlying correlation
# coefficient of 0.1
sapply(fcor(100, 50, fgen, 0.1), cor)[2,]
#>  [1]  0.094406279  0.284145751  0.192933797  0.158482675  0.057744828
#>  [6]  0.083726118  0.311058839  0.163806143  0.063969503 -0.051479752
#> [11]  0.063455103  0.248141890  0.225047262 -0.169708868  0.134908044
#> [16]  0.043553690  0.105127515  0.135446660  0.034585970  0.195204328
#> [21]  0.122711017  0.136277440  0.204716668  0.259685417  0.221983749
#> [26] -0.023070684  0.183586680  0.166254882  0.030147390  0.229761794
#> [31]  0.014598950  0.040468710  0.193315940  0.132643005  0.105137032
#> [36]  0.147129518  0.264126411  0.113230597  0.009142460  0.051503916
#> [41]  0.034976333  0.105221879  0.190592943  0.103760516  0.159057939
#> [46] -0.019006156 -0.008168006 -0.002263670 -0.005456012  0.201508298

只有 100 个观测值，相关系数变化很大。生成具有 10M 个观测值的样本：

cor(fcor(1e7, 1, fgen, 0.1)[[1]])[2]
#> [1] 0.09926747

也可以使用负相关系数：

cor(fcor(1e6, 1, fgen, -0.1)[[1]])[2]
#> [1] -0.1004658