分析の実験をするときに、相関を持った2つの変数の乱数を発生させたいときがよくある。その関数をRで作った。ここで、ある相関を持った乱数を発生させるには以下の性質を用いる。すなわち、ある確率分布Aに従う乱数\(x\)に対して相関\(r\)を持つ乱数\(y\)は以下のように定義される。なお\(z\)は確率分布Aに従う乱数である。
\[ y = r \times x + \sqrt{(1 - r^2)} \times z \]
ただし、確率分布Aに従う完全な乱数を発生させるのは不可能である。そのことが原因で\(x\)と\(y\)の相関係数は\(r\)と若干異なってしまう。そこで、何回も\(y\)を生成して最も\(r\)に近くなった時の\(y\)を用いることにした。
以下は、そのコードである。
cor_get <- function(r, n, family = 'normal', rep = 1000, seed = 1234,
mean = 0, sd = 1, lambda = 1, size = 10, prob = 0.5,
min = 0, max = 1){
set.seed(seed)
if(family == 'normal'){
x <- rnorm(n, mean, sd)
y <- list()
cor <- NULL
for(i in 1:rep){
z <- rnorm(n, mean = mean, sd = sd)
y[[i]] <- r * x + sqrt(1 - r^2)*z
cor[i] <- cor(x, y[[i]])}
}
else if(family == 'binomial'){
x <- rbinom(n, size, prob)
y <- list()
cor <- NULL
for(i in 1:rep){
z <- rbinom(n, size, prob)
y[[i]] <- r * x + sqrt(1 - r^2)*z
cor[i] <- cor(x, y[[i]])}
}
else if(family == 'poisson'){
x <- rpois(n, lambda)
y <- list()
cor <- NULL
for(i in 1:rep){
z <- rpois(n, lambda)
y[[i]] <- r * x + sqrt(1 - r^2)*z
cor[i] <- cor(x, y[[i]])}
}
else if(family == 'uniform'){
x <- runif(n, min, max)
y <- list()
cor <- NULL
for(i in 1:rep){
z <- runif(n, min, max)
y[[i]] <- r * x + sqrt(1 - r^2)*z
cor[i] <- cor(x, y[[i]])}
}
else{warning('Warning: No such distribution', call. = F)}
diff <- NULL
for(i in 1:rep){
diff[i] <- abs(cor[i] - r)
}
library(tidyverse)
number <-
diff %>%
tibble() %>%
mutate(id = 1:rep) %>%
filter(. == min(diff)) %>%
select(id) %>%
as.numeric()
df <- tibble(x, y[[number]])
colnames(df) <- c("x", "y")
return(list(df, cor(df[1], df[2])))
}第1引数のrは目標の相関係数、第2引数のnはサンプル数、第3引数のfamilyは確率分布A、第4引数のrepは\(y\)の生成回数（これを多くすれば精度がよくなる可能性が高まる）、第5引数のseedは乱数のseedである。それ以降の引数は、各確率分布に必要なパラメータである。ここでは、その確率分布として正規分布、二項分布、ポワソン分布、一様分布を設定することができる。 この関数はリスト形式でそのデータフレームと相関関係を出力する。なお、この関数を使うにはtidyverseをインストールしておく必要があることには留意されたい。