如何使用for循环更新协方差矩阵中的对角线？

Question

我已经使用for规范为两个变量创建了模拟数据，并希望在循环中将这些变量0、.5、.7和.9关联起来。但是，每次我运行for循环时，我只能关联.9上的值，而不能关联其他任何关联条件下的值。

library(MASS) #library I needed to create simulated data with mvrnorms

num_iter <- 75
N <- 30                       # setting my sample size
mu <- c(50.5, 10.5)           # setting the std
R <- c(0,.5,.7,.9)            # this vector defines the different correlation conditions I will add

# saving files
dir.create("simulated1data") # This creates a directory to store files

# performing 75 iterations and so there should be 75 data files in the folder I made
for(i in 1:num_iter){
  for(j in 1:4){
    cov <- matrix(c(1,R[j],R[j],1),2,2)
    x <- mvrnorm(N,mu,cov)
    write.table(x, file=paste("simulated1data/simdata_",i,"_",j,".txt",sep="")) # writing to separate txt file
  }
}

根据我的理解，我的(对于1:4中的j)没有适当地遍历我的R向量中的所有jth值，这就是为什么X中的变量总是在.9上相关的原因。有人知道怎么解决这个问题吗？谢谢您抽时间见我!

Answer 1

要分配R的值，请预先创建cov矩阵并使用逻辑索引矩阵imat。

第一个代码块如问题中所示。

library(MASS) #library I needed to create simulated data with mvrnorms

num_iter <- 75
N <- 30                       # setting my sample size
mu <- c(50.5, 10.5)           # setting the std
R <- c(0, 0.5, 0.7, 0.9)      # this vector defines the different correlation conditions I will add

这是为了在我的系统上测试。

# saving files
dirsimdata <- "~/tmp/simulated1data"
dir.create(dirsimdata) # This creates a directory to store files

现在是cov和imat矩阵。

# index matrix used to assign values from R
imat <- matrix(c(FALSE, TRUE, TRUE, FALSE), nrow = 2)
# start with all 1's
cov <- matrix(1, nrow = 2, ncol = 2)

最后是双for循环。

# performing 75 iterations and so there should be 75 data files in the folder I made
for(i in 1:num_iter){
  for(j in 1:4){
    cov[imat] <- R[j]
    x <- mvrnorm(N, mu, cov)
    flname <- paste0("simdata_", i, "_", j, ".txt")
    flname <- file.path(dirsimdata, flname)
    write.table(x, file = flname) # writing to separate txt file
  }
}

Answer 2

我在你的代码中没有发现任何错误。您错误地将mu标识为标准差，但它是每个变量的平均值，而R是协方差而不是相关性。您可以在协方差矩阵中设置1处每个变量的标准差。如果我在进入循环之前设置num_iter <- 2并使用set.seed(42)，考虑到样本大小只有30，我得到了合理的相关性：

cor(read.table("simulated1data/simdata_1_1.txt"))
#          V1       V2
# V1 1.000000 0.204011
# V2 0.204011 1.000000
cor(read.table("simulated1data/simdata_1_2.txt"))
#           V1        V2
# V1 1.0000000 0.2706851
# V2 0.2706851 1.0000000
cor(read.table("simulated1data/simdata_1_3.txt"))
#           V1        V2
# V1 1.0000000 0.6727047
# V2 0.6727047 1.0000000
cor(read.table("simulated1data/simdata_1_4.txt"))
#           V1        V2
# V1 1.0000000 0.9306898
# V2 0.9306898 1.0000000
cor(read.table("simulated1data/simdata_2_1.txt"))
#            V1         V2
# V1 1.00000000 0.06184222
# V2 0.06184222 1.00000000
cor(read.table("simulated1data/simdata_2_2.txt"))
#           V1        V2
# V1 1.0000000 0.3686962
# V2 0.3686962 1.0000000
cor(read.table("simulated1data/simdata_2_3.txt"))
#           V1        V2
# V1 1.0000000 0.7660853
# V2 0.7660853 1.0000000
cor(read.table("simulated1data/simdata_2_4.txt"))
#           V1        V2
# V1 1.0000000 0.8589621
# V2 0.8589621 1.0000000