如何在要重复 5000 次的 for 循环中找到相关系数？并保存统计信息答案

【问题标题】：how to find correlation coefficient in a for loop that is to be repeated 5000 times? and save the statistic如何在要重复 5000 次的 for 循环中找到相关系数？并保存统计信息
【发布时间】：2019-06-24 19:17:14
【问题描述】：

对于 2 个独立的正态分布变量 x 和 y，它们是使用 x = rnorm(50) 和 y = rnorm(50) 找到的。计算相关性 5000 次，每次保存结果。计算绝对值大于 0.3 的相关性的可能性有多大？（默认 set.seed(42) 并绘制系数分布的直方图）

这是我迄今为止尝试过的......

set.seed(42)

n <- 50 #length of random sequence
x_norm <- rnorm(n)
y_norm <- rnorm(n)

nrun <- 5000  
corr <- numeric(nrun) 

for (i in 1:nrun) {
corrxy <-  cor(x_norm,y_norm)
corr[i] <- sum(abs(corrxy > 0.3)) / n  #save statistic in the vector
    }

hist(corr)

预计我会在 [i] 中保存 5000 个不同的系数编号，并且当使用 hist(0) 绘制时，这些系数应遵循近似正态分布。但我不明白 for 循环是如何工作的，以及如何合并大于 0.3 的系数值。

【问题讨论】：

corrxy > 0.3 为您提供布尔输出，TRUE 或 FALSE，对于 TRUE 的计算结果为 1，对于 FALSE 的计算结果为 0。
很好的问题。两个rnorm() 不应该在循环内吗？否则cor() 将始终相同。

标签： r for-loop correlation normal-distribution

【解决方案1】：

我想你快到了。您只需要在 for 循环内外移动一些代码。

您希望每次循环运行都有新数据（否则您会获得相同的相关性 5000 次），并且每次循环运行时都需要保存相关性。这会产生一个包含 5000 个相关性的向量，您可以使用它来查看在 for 循环之外高于 0.3 的相关性比例（除以运行次数，而不是观察次数）。

编辑：在绝对函数的包围中需要进行最后一次校正。您想找到绝对相关性 > .3 而不是 corrxy > .3 的绝对值。

set.seed(42)

n <- 50 #length of random sequence


nrun <- 5000  
corrxy <- numeric(nrun) # The correlation is the statistic you want to save

for (i in 1:nrun) {
x_norm <- rnorm(n) # Compute a new dataset for each run (otherwise you get the same correlation)
y_norm <- rnorm(n)
     corrxy[i] <-  cor(x_norm,y_norm) # Calculate the correlation

}
hist(corrxy)
sum(abs(corrxy) > 0.3) / nrun # look at the proportion of runs that have cor > .3

下面是 5000 个相关性的结果直方图。高于|.3|的相关比例在这种情况下是 0.034。

【讨论】：

哇，非常感谢。我终于明白为什么总和必须在循环之外
不客气，如果你喜欢这个答案，你可以点击投票计数器下方的复选标记接受它。

【解决方案2】：

这是在不显式调用循环的情况下进行此类模拟的另一种方法：

首先定义你的模拟：

my_sim <- function(n) { # n is the norm distribution size
  x <- rnorm(n)
  y <- rnorm(n)

  corrxy <- cor(x, y)
  corrxy # return the correlation (single value)
}

现在我们可以用replicate()多次调用这个函数：

set.seed(123)
nrun <- 10
my_results <- replicate(nrun, my_sim(n=50))
#my_results
# [1] -0.0358698314 -0.0077403045 -0.0512509071 -0.0998484901  0.1230261286  0.1001124010 -0.0002023124
# [8]  0.2017120443  0.0644662387  0.0567232640

现在在my_results 中，您可以获得每个模拟的所有相关性（例如，只有 10 个）。

你可以计算你的统计数据：

sum(abs(my_results)> 0.3) / nrun  # nrun is 10

或情节：

hist(my_results)

【讨论】：