LKJ分布に従う相関行列
相関係数の3次元散布図
相関行列の相関係数は、相関行列の行列式が正値であるという制約に従って分布している。いま、変数が３個、３次の相関行列について考える。
３次の相関行列を次式のように置く。
相関行列は対称である。すなわち、
が成り立っている。したがって、
に含まれる相関係数は次の3個
、
、
であると考えることができる。
これらの相関係数は、次式
を満たす範囲に存在する。
いま、行列式を、、の関数として見て、行列式が正の場合は赤の点で、非正であれば青の点で行列式の値の３D散布図を描くと図１のようになる。
図１
図１はリスト１のプログラムで描かれている。プログラムの実行により表示されているグラフは、グラフのドラッグにより回転することができる。なお、プログラムのソースファイルは圧縮ファイルlkj_d_scatter_file.zipにまとめられている。
3個、、の値が３変数の相関行列の相関係数であるためには、図１において赤い点で表されている領域内の値でなければならない。
この制約を満たす相関行列の分布にLKJ相関行列確率分布があるが、これはパラメータ１つ
をもち、次式で与えられる。
の値に対して相関係数、、の分布を描くプログラムをリスト２のように作成した。プログラムファイルは圧縮ファイルlkj_d_scatter_file.zipにまとめてある。
LKJ分布に従う相関行列は、次のStanスクリプトによりサンプルを得ている。
data {
real<lower = 0.0> eta;
}
parameters {
corr_matrix[3] corr;
}
model {
corr
~ lkj_corr(eta);
}
リスト２のプログラムScatter_LKJ.pyを実行するとの値の入力が求められる。
(py39) PS D:\lkj_d_scatter_file>
python Scatter_LKJ.py
eta = 3.0
「eta = 3.0」と3.0を入力してEnterキーを押すと計算が始まる。計算が終了すると図２のグラフが表示される。
図２
グラフは、相関係数の値を0.2幅で区切ってStanによるサンプリングの集計を行い、各ブロック内のサンプル数が最大数の1／2以上の場合は赤い点で、1／2以下で1／4以上の場合は緑の点で、1／4以下であれば薄青の点で表している。グラフ上をドラッグすればグラフを回転して、図３のようなグラフも得ることができる。
図３
図２，図３から原点に近いところの密度が高いことが分かる。これは、式（１）の密度関数は行列式の値が大きくなると式（１）の値も大きくなることによる。
「eta = 1」と入力すると図４のグラフを得る。
図４
のときは、式（１）の値は定数となり、一様分布を表す。図４では、緑の点と薄青の点も分布の周辺に見られるが、これは相関係数を幅0.2のブロックで区切って処理しているためであり、図４は一様分布を反映していると見ることができる。
「eta = 0.3」と入力すると、図５の散布図を得る。
図５
のときは、式（１）の値は行列式の値が０に近付くと大きくなる。行列式の値が正という制約条件の下で相関係数の値が大きくなると式（１）の値も大きなる。図５では４つの角に赤い点が見られる。
なお、PythonおよびPyStanの解説を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
PythonとPyStanのインストールの説明ウェブサイトも用意している。
リスト１ 行列式の符号の散布図描画Pythonスクリプト（ファイル名Scatter_Det.py）
import numpy as np
import matplotlib.pyplot as plt
XYZ_pos = []
XYZ_npos = []
for i12 in np.arange(-0.9, 1.0, 0.2):
for i13 in np.arange(-0.9,
1.0, 0.2):
for
i23 in np.arange(-0.9, 1.0, 0.2):
A = np.array([[1, i12, i13],
[i12, 1, i23],
[i13, i23, 1 ]])
v = np.linalg.det(A)
if v > 1.0e-9:
XYZ_pos.append([i12, i13, i23])
else:
XYZ_npos.append([i12, i13, i23])
XYZ_pos = np.array(XYZ_pos)
XYZ_npos = np.array(XYZ_npos)
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(projection =
'3d')
ax.scatter((XYZ_pos.T)[0],
(XYZ_pos.T)[1], (XYZ_pos.T)[2], color = 'r',
label = 'Positive')
ax.scatter((XYZ_npos.T)[0],
(XYZ_npos.T)[1], (XYZ_npos.T)[2], color = 'b',
label = 'NotPositive')
ax.set_xlabel(r'$\rho_{12}$', fontsize
= 18)
ax.set_ylabel(r'$\rho_{13}$', fontsize
= 18)
ax.set_zlabel(r'$\rho_{23}$', fontsize
= 18)
plt.legend()
plt.show()
リスト２ 3次の相関行列のLKJ分布（ファイル名Scatter_LKJ.py）
import numpy as np
import pystan
import matplotlib.pyplot as plt
StanCode = """
data {
real<lower = 0.0> eta;
}
parameters {
corr_matrix[3] corr;
}
model {
corr
~ lkj_corr(eta);
}
"""
eta = float(input('eta = '))
sm = pystan.StanModel(model_code =
StanCode)
fit = sm.sampling(data = {'eta': eta},
iter = 100_000, n_jobs = 1)
print(fit)
def idx_r(r):
return int((r + 1) *
5)
r_smpls = fit['corr']
cnt_r = np.zeros((10, 10, 10))
for v in r_smpls:
cnt_r[idx_r(v[0,1])][idx_r(v[0,2])][idx_r(v[1,2])] += 1
cnt_max = cnt_r.max()
print(cnt_max)
def r_from_idx(i):
return 2 * (i / 10) - 0.9
Low_cnt = []
Med_cnt = []
High_cnt = []
for i12 in range(10):
for i13 in range(10):
for
i23 in range(10):
if cnt_r[i12][i13][i23] > cnt_max/2:
High_cnt.append([r_from_idx(i12),
r_from_idx(i13),
r_from_idx(i23)])
elif cnt_r[i12][i13][i23] > cnt_max/4:
Med_cnt.append([r_from_idx(i12),
r_from_idx(i13),
r_from_idx(i23)])
elif cnt_r[i12][i13][i23] > 0:
Low_cnt.append([r_from_idx(i12),
r_from_idx(i13),
r_from_idx(i23)])
High_cnt = np.array(High_cnt)
Med_cnt = np.array(Med_cnt)
Low_cnt = np.array(Low_cnt)
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(projection =
'3d')
if len(High_cnt) > 0:
ax.scatter(High_cnt.T[0],
High_cnt.T[1], High_cnt.T[2],
color = 'r', label = '>max/2')
if len(Med_cnt) > 0:
ax.scatter(Med_cnt.T[0],
Med_cnt.T[1], Med_cnt.T[2],
color = 'g', label = 'max/4 ~ max/2')
if len(Low_cnt) > 0:
ax.scatter(Low_cnt.T[0],
Low_cnt.T[1], Low_cnt.T[2],
color = 'b', alpha = 0.3, label = '<max/4')
ax.set_xlabel(r'$\rho_{12}$', fontsize
= 14)
ax.set_ylabel(r'$\rho_{13}$', fontsize
= 14)
ax.set_zlabel(r'$\rho_{23}$', fontsize
= 14)
plt.title(r'$\eta = ${}'.format(eta),
fontsize = 20)
plt.legend()
plt.show()