ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【データ解析】相関係数を式で考えてみる

相関係数の話の続きです.
chaos-kiyono.hatenablog.com

 今回は確率変数を使って考えてみます.

Y=aX+bWの係数abと,XYの散布図の関係.

相関って線形な関係

 平均0,分散\sigma^2の確率変数Xと,これと独立な平均0,分散\sigma^2の確率変数Wを仮定します.XWの分布については,分散が有限ということだけで形は指定しません.

 Xと線形な依存性がある確率変数Yを,定数abを使って次の形で表します.

 \begin{eqnarray} Y = a X + b W\end{eqnarray}

この場合,a \neq 0のとき,YにはXと線形な関係があり,Wの係数の|b|が大きいほどその関係性が薄くなります.つまり,XYの相関は,aの絶対値が,bの絶対値よりも大きいほど強いと言えそうです.今回は実際に相関係数を計算してこのことを足しためてみます.

以下の計算で使う基本事項として,Yの期待値と分散を計算しておきます.Yの期待値は,

 \begin{eqnarray} \mathbb{E} [Y] &=& \mathbb{E} [a X + b W ] \\
 &=& a \mathbb{E} [X] + b \mathbb{E} [W ] \\ 
&=& a \cdot 0 + b  \cdot 0 \\ 
&=& 0 \end{eqnarray}

Yの分散は,

 \begin{eqnarray} {\rm Var} [Y] &=& {\rm Var} [a X + b W ] \\
&=& {\rm Var} [a X]  + {\rm Var} [b W ] \\
&=& a^2 {\rm Var} [X]  + b^2 {\rm Var} [ W ] \\
 &=& a^2 \sigma^2 + b^2 \sigma^2 \\ 
&=& \left( a^2 + b^2 \right) \sigma^2 \end{eqnarray}

です.

 上の例では,a=0ならXYは独立 (WXとは独立なので)です.また,b=0ならXYは単に比例しています.

 相関係数rの定義は,

 r = \frac{\displaystyle \mathbb{E} \big[ \left(X- \mathbb{E} \left[ X \right] \right)\left(Y- \mathbb{E} \left[ Y \right] \right)\big]}{\displaystyle \sqrt{\mathbb{E}  \big[ \left(X- \mathbb{E} \left[ X \right] \right)^2\big]} \sqrt{\mathbb{E}  \big[ \left(Y- \mathbb{E} \left[ Y \right] \right)^2\big]}}

です.

 上で定義した,XY相関係数を計算すると,

 \begin{eqnarray} r &=& \frac{ \mathbb{E} \big[ X \left(a X + bW \right) \big]}{\sqrt{\mathbb{E}  \big[ X^2 \big]} \sqrt{\mathbb{E}  \big[ \left(a X + b W \right)^2\big]}} \\
&=& \frac{ \mathbb{E} \big[ a X^2 + b X W \big]}{\sqrt{{\rm Var} \big[ X \big]} \sqrt{{\rm Var} \big[ a X + b W \big]}} \\
&=& \frac{ a\, \sigma^2}{\sigma \sqrt{\left( a^2 + b^2 \right) \sigma^2}} \\
&=& \frac{a}{ \sqrt{a^2+b^2}} 
\end{eqnarray}

になります.最後の式は,直角三角形をイメージしてもらって,直角をはさむ2辺の長さをabとしたときに,斜辺の長さ\sqrt{a^2+b^2}aの比になっています.

直角三角形の辺の長さ.

 上の結果を使うと,a=0のとき,つまり,YXに相関がないケースでは,r=0になります.

 また,b=0のとき,つまり,YXの定数倍のケースでは,r=1か,r=-1になります.

 そして,相関係数の範囲は,

 a > 0 のとき,0 < r \le 1

 a < 0 のとき,-1 \le r < 0

になることが分かります.

 ということで,相関係数がとる値の範囲は,-1 \le r \le 1です.

 また,「相関係数の推定量

 \hat{r} = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)\left(y_i-\bar{y} \right)}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)^2} \sqrt{\frac{1}{N}\sum_{i=1}^N \left(y_i-\bar{y} \right)^2}}

についても,-1 \le \hat{r}  \le 1になります.上の計算はこの証明にはなっていないので,この範囲については自分で考えるか,インターネットで検索して調べてみてください.当然,いろんな教科書に載っています.

相関する確率変数の生成

 上の計算を簡単にすることを考えます.簡単化のため,互いに独立な確率変数XWは,平均0,分散1とします.そして,Yの分散も1にします.このとき,

  \begin{eqnarray} \mathbb{E} [Y^2] &=& \mathbb{E} [\left(a X + b W \right)^2 ] \\
 &=& a^2 + b^2 = 1 \end{eqnarray}

が成り立つ必要があるので,b^2=1-a^2となります.

 以上の条件でもう一度,相関係数rを計算してみると,

 \begin{eqnarray} r &=& \frac{ a}{\sqrt{a^2 +b^2}} \\
&=& \frac{a}{\sqrt{a^2 + \left( 1 - a^2 \right)}} \\
&=& a
\end{eqnarray}

になります.結局, a = rになりました.

 ということで,相関係数rの確率変数XYを作りたいとき,まず,平均0,分散1の互いに独立な確率変数XWを用意して,

 Y = r X + \sqrt{1-r^2} W

とすれば良いと言うことです.

Rで互いに相関する正規乱数の生成

 真の相関係数rになる正規乱数をRで生成してみます.使ったスクリプトはこれです.

# データの長さ
N <- 1000
# 真の相関係数
r <- 0.5
# 2変量乱数の生成
x <- rnorm(N) # 正規乱数
w <- rnorm(N) # 正規乱数
y <- r*x+sqrt(1-r^2)*w
# 相関係数の推定 
r.est <- cor(x,y)
###
# 散布図を描画
y.max <- max(abs(c(x,y)))
par(pty="s",cex.main=1.5)
plot(x,y,xlim=c(-y.max*1.3,y.max*1.3),ylim=c(-y.max*1.3,y.max*1.3),pch=16,col="blue",main=sprintf("r (est) = %.2f",r.est))
if(r >= 0){
 text(y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=2,cex=1.4,col=2)
 text(-y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=4,cex=1.4,col=1)
}else{
 text(-y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=4,cex=1.4,col=2)
 text(y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=2,cex=1.4,col=1)
}
Rスクリプトの実行結果