（著）山たー
分散共分散行列を相関行列や偏相関行列に変換する関数がPythonに無かったので実装してみました。こういうところでMATLABとの差があるのでしょうか。
分散共分散行列
先に分散共分散行列(variance–covariance matrix)の定義とPythonでの計算方法について説明しておきます。
定義
変数の列ベクトル$X$を $$\boldsymbol{X} ={\begin{bmatrix}X_{1}\\ \vdots \\X_{n}\end{bmatrix}}$$ とします。このとき、分散共分散行列$\Sigma(=\textrm{Cov}(\boldsymbol{X})=\textrm{Var}(\boldsymbol{X}))$は$X_{i}$の平均を
$\mu_{i}=\textrm{E}(X_{i})$ として、 $$ \Sigma ={\begin{bmatrix}\textrm {E} [(X_{1}-\mu _{1})(X_{1}-\mu _{1})]&\textrm {E} [(X_{1}-\mu _{1})(X_{2}-\mu _{2})]&\cdots &\textrm {E} [(X_{1}-\mu
_{1})(X_{n}-\mu _{n})]\\\\\textrm {E} [(X_{2}-\mu _{2})(X_{1}-\mu _{1})]&\textrm {E} [(X_{2}-\mu _{2})(X_{2}-\mu _{2})]&\cdots &\textrm {E} [(X_{2}-\mu _{2})(X_{n}-\mu _{n})]\\\\\vdots
&\vdots &\ddots &\vdots \\\\\textrm {E} [(X_{n}-\mu _{n})(X_{1}-\mu _{1})]&\textrm {E} [(X_{n}-\mu _{n})(X_{2}-\mu _{2})]&\cdots &\textrm {E} [(X_{n}-\mu _{n})(X_{n}-\mu
_{n})]\end{bmatrix}} $$ と表せます。これは $$ \Sigma =\textrm {E} \left[\left(\boldsymbol{X} -\textrm {E} [\boldsymbol{X} ]\right)\left(\boldsymbol{X} -\textrm {E} [\boldsymbol{X} ]\right)^{\rm {T}}\right] $$
とも表せます。
Pythonでの求め方
データから分散共分散行列を求めるには、np.cov()を使えば面倒な行列計算をせずに済みます。以下ではnp.cov()などを用いて分散共分散行列を取得し、それを相関行列や偏相関行列に変換したい場合を考えます。
変換する関数の実装
分散共分散行列を相関行列に変換
相関行列(Correlation matrix)を$\textrm{Corr}$とすると、 $$ \textrm{Corr}(\boldsymbol{X} )=\left(\textrm{diag}(\Sigma)\right)^{-{\frac {1}{2}}} \Sigma \left(\textrm{diag} (\Sigma )\right)^{-{\frac {1}{2}}} $$
となります。ただし、$\textrm{diag}(\Sigma)$は$\Sigma$の対角成分であり、$\left(\textrm{diag}(\Sigma)\right)^{-{\frac {1}{2}}}$は$\Sigma$の対角成分の各要素を$-\frac{1}{2}$乗した行列を意味します。
import numpy as np
#分散共分散行列を相関行列に変換
def cov2corr(cov):
D=np.diag(np.power(np.diag(cov),-0.5))
corr=np.dot(np.dot(D,cov),D)
return corr
なお、データから直接、相関行列を求めるにはnp.corrcoef()を使います。
分散共分散行列を偏相関行列に変換
偏相関行列(Partial correlation matrix)を$P=[p_{ij}]$とします。$\Sigma$の逆行列を$\Omega=[\omega_{ij}]$とすると、 $$ p_{ij}=-\frac{\omega_{ij}}{\sqrt{\omega_{ii}\omega_{jj}}} $$
が成り立ちます。これは偏相関係数の、分散共分散行列の逆行列を用いた求め方です。これを行列化すると、 $$ P(\boldsymbol{X} )=-\left(\textrm{diag}(\Omega)\right)^{-{\frac {1}{2}}} \Omega \left(\textrm{diag} (\Omega )\right)^{-{\frac {1}{2}}} $$
となります。ただし、このままでは対角成分が$-1$になるので、実装の上は対角成分がすべて$2$の行列を加えることで、対角成分を$1$にしています。
import numpy as np
#分散共分散行列を偏相関行列に変換
def cov2partialcorr(cov):
#逆行列を求める
omega=np.linalg.inv(cov)
D=np.diag(np.power(np.diag(omega),-0.5))
partialcorr=-np.dot(np.dot(D,omega),D)
#対角成分を-1から1にする
partialcorr+=2*np.eye(cov.shape[0])
return partialcorr
使用例
検証という意味も込めて使用例を書いておきます。偏相関行列の確認のため、MATLABのpartialcorrと比較します。
import numpy as np
#分散共分散行列を相関行列に変換
def cov2corr(cov):
D=np.diag(np.power(np.diag(cov),-0.5))
corr=np.dot(np.dot(D,cov),D)
return corr
#分散共分散行列を偏相関行列に変換
def cov2partialcorr(cov):
#逆行列を求める
omega=np.linalg.inv(cov)
D=np.diag(np.power(np.diag(omega),-0.5))
partialcorr=-np.dot(np.dot(D,omega),D)
#対角成分を-1から1にする
partialcorr+=2*np.eye(cov.shape[0])
return partialcorr
#分散共分散行列
cov=np.array([[ 0.2516, 0.3014, 0.0507, 12.5960],
[ 0.3014, 52.0622, 0.2069, 17.5152],
[ 0.0507, 0.2069, 0.2267, 2.7273],
[12.5960, 17.5152, 2.7273, 706.0404]])
#相関行列
corr=cov2corr(cov)
#偏相関行列
partialcorr=cov2partialcorr(cov)
#結果表示
print("相関行列")
print(corr)
print("偏相関行列")
print(partialcorr)
結果
相関行列 [[ 1. 0.08327731 0.2122887 0.94506691] [ 0.08327731 1. 0.06022454 0.09135642] [ 0.2122887 0.06022454 1. 0.21557201] [ 0.94506691 0.09135642 0.21557201 1. ]]
偏相関行列 [[ 1. -0.0105331 0.02723055 0.94213313] [-0.0105331 1. 0.04194991 0.03687756] [ 0.02723055 0.04194991 1. 0.04516695] [ 0.94213313 0.03687756 0.04516695 1. ]]
コメントをお書きください
辛い (月曜日, 17 6月 2019 22:03)
初めまして。こちらのサイトで勉強させてもらっています。
お聞きしたいことがあるのですが、相互相関係数を求める場合は分散共分散行列の対角成分に一つでもマイナス項がある場合、偏相関係数を求める場合は分散共分散行列の逆行列のすべての対角成分の内マイナスが奇数個の場合はさせている処理的にエラーが出ると思うのですが、その場合は各相関係数は求められないということでいいのでしょうか・・・
例えば
cov=np.array.([[1,2,5],
[1,-1,1],
[0,1,2]])
の場合等・・・
通りすがり (金曜日, 21 6月 2019 13:48)
>辛いさん
分散共分散行列の各要素は共分散(特に対角成分は分散そのもの)なのでマイナスにならないのでは？
通りすがり (金曜日, 21 6月 2019 13:52)
あ、共分散は負の値になりますね。
ただ、対角の分散は正ですね。