相関のある2つの変数の乱数の発生方法

Rで相関関係のある2つの変数の乱数を発生させる時、いつも以下の方法に頼っていました。

n<-10000
b<-0.9 
set.seed(1234)
x<-rnorm(n)
e<-rnorm(n)
y<-x*b+e
cor(x,y)

この時のbをいじる事で相関の強さを調整していました。なんと効率の悪い方法でしょう。

しかし、 以下の式を用いる事で、目標となる相関係数rの関係にある2つの変数（xとy）の乱数を発生させる事ができるようです。

 

なお、ここでxとzは独立な乱数です。なぜそうなるかは置いておいて（勉強します）、非常に便利に相関関係のある2つの変数の乱数を発生させる事ができるようです。

早速、Rで試します。なお、ここでは正規分布に従う乱数を用いています。

n<-10000
r<-0.7

set.seed(1234)
x<-rnorm(n)
z<-rnorm(n)
y<-x*r+sqrt(1-r^2)*z

cor(x,y)

この結果、相関係数は0.6973729となりました。

おおよそ、目的の相関係数0.7に近づいたと言えるでしょう。

ただし、 若干の違い（今回は（0.7 - 0.6973729）の分だけ）があるので、どれくらいの違いがどれくらいの頻度で生まれるのか検証してみましょう。

目標の相関係数は0.7、繰り返し数は1000回です。

おおよそですが、目標である相関係数0.7を中心とした正規分布になっています。

ちなみに、0.69以下となる確率は0.028%、0.71以上となるの確率は0.019%でした。目標の相関係数から0.01以上離れることはほとんどないと言えるようです。

なお、コードは以下の通りです。

n<-10000
r<-0
rep<-1000 #繰り返し数
rec<-numeric(rep)

for(i in 1:rep){
	x<-rnorm(n)
    z<-rnorm(n)
    y<-x*r+sqrt(1-r^2)*z
    rec[i]<-cor(x,y)
}

library(ggplot2)
g<-ggplot()+geom_histogram(mapping=aes(x=rec),binwidth=0.001)
print(g)

p_under<-ifelse(rec<=(r-0.1),1,0)
p_over<-ifelse((r+0.1)<=rec,1,0)

p_under<-sum(p_under)/length(rec)
p_over<-sum(p_over)/length(rec)

print(p_under)
print(p_over)

それでは、目標となる相関係数の値を変化させたらどうなるでしょうか。

目標となる相関係数を-1〜1まで0.25ずつ増加させてみました。（実際は-1〜0まで0.25ずつ増加させたもののヒストグラムを出力し、その後、0〜1まで0.25増加させたものを出力した）

目標の相関係数の絶対値が1に近づくほど、その誤差は少なくなっています。目標の相関係数が1もしくは-1の時は、必ず目標の相関係数を得られるようです。

なおコードは以下の通りです。

n<-10000
r<-seq(0,1,0.25)
rep<-1000 #繰り返し数
rec<-numeric(rep*length(r))
a<-0

for(r in r){
	
for(i in 1+a:rep){
	x<-rnorm(n)
    z<-rnorm(n)
    y<-x*r+sqrt(1-r^2)*z
    rec[i]<-cor(x,y)
    }
    
    rep<-rep+1000
    a<-a+1000
}

library(ggplot2)
g<-ggplot()
g<-g+geom_histogram(mapping=aes(x=rec[1:1000]),binwidth=0.001,position="identity")+geom_histogram(mapping=aes(x=rec[1001:2000]),binwidth=0.001,position="identity")+geom_histogram(mapping=aes(x=rec[2001:3000]),binwidth=0.001,position="identity")+geom_histogram(mapping=aes(x=rec[3001:4000]),binwidth=0.001,position="identity")+geom_histogram(mapping=aes(x=rec[4001:5000]),binwidth=0.001,position="identity")
print(g)

目標である相関係数が0である場合が最も目標とする相関係数を得られにくいようです。目標である相関係数が0である場合、それを1000回繰り返したところ、-0.01以下となる確率は16%、0.01以上となる場合は17.9%でした。

以上から、上記をRで試した時、目標とする相関係数の絶対値が大きければ大きいほど、その目標値に近い値を得られると言えるでしょう。

ただし、最も目標値が得られにくいと考えられる相関係数0の場合でも、それが0.1以上離れる事はありませんでしたから、これくらいの誤差は許容できる！という場合は、とても良い方法であると言えると思います。