お元気ですか？最近読んだ『違国日記』という漫画がよかったのでお勧めです。しばらくしたらどっかで感想を書こうと思います。
https://www.shodensha.co.jp/ikokunikki/
https://www.amazon.co.jp/dp/B077GQL19W/ref=dp-kindle-redirect?_encoding=UTF8&btkr=1
偏相関
さて、最近、偏相関というものについて定義と計算方法をうにうにと考えていました。偏相関は相関関係から疑似相関を引っこ抜こうとがんばる作業です。ネット上の情報は分かるような分らんような感じだったので私の理解をメモ代わりに書いておきます。本を買ってもよかったのですがそこまでするようなことか...と思って自分でやってみると時間がかかってしまった。
相関係数
まず相関係数（相関行列）からおさらいをする。以下ではすべて実数の範囲で考える。
分散共分散
実数に値を持つ分散が有限な N 個の確率変数 (X^1 \cdots X^N) に対して、分散共分散行列 \Sigma \in \mathrm{Mat}(N, N) が
\Sigma^{ij} = \mathrm{Cov}(X^i, X^j)\ \ \ \ (1\leq i, j \leq N)
で定義できる。ただし \mathrm{Cov}(X^i, X^j) = E[(X^i - \mu^i)(X^j - \mu^j)] である。いま、各確率変数X^iを\mu^iだけ平行移動したものを新たにX^i（つまり\mu^i = 0）とすれば
\begin{align*}
\Sigma
&= (E[X^i X^j])_{ij} \\
&= E[X^t X]
\end{align*}
となる。
相関行列
相関行列\rho \in \mathrm{Mat}(N, N; \mathbf{R})が
\begin{align*}
\rho
&= \left( \dfrac{\mathrm{Cov}(X^i, X^j)}{\sigma(X^i) \sigma(X^j)} \right)_{ij} \\
&= \left( \dfrac{E[X^i X^j]}{ \sqrt{E[(X^i)^2]} \sqrt{E[(X^j)^2]}} \right)_{ij}
\end{align*}
で定義される。これは共分散を標準偏差で正規化したものである。実際、| \rho^{ij} | \leq 1が簡単な計算でわかる。
さて、確率変数X^iを独立にM回測定したときの測定値をx^i = (x^i_1 \cdots x^i_M)^tとすれば、相関係数は
\begin{align*}
\rho^{ij}
&= \dfrac{\dfrac{1}{M}\sum_{m = 1}^M x^i_m x^j_m }{\sqrt{\dfrac{1}{M} \sum_{m = 1}^M (x^i_m)^2}\sqrt{\dfrac{1}{M} \sum_{m = 1}^M (x^i_m)^2}} \\
&= \dfrac{\langle x^i, x^j \rangle}{\|x^i\| \|x^j\|}
\end{align*}
と書ける。明らかに|\rho^{ij}| \leq 1である。また、相関行列を分散共分散行列 \Sigma を用いて書けば
\begin{align*}
\rho
&= \left( \dfrac{\langle x^i, x^j \rangle}{\|x^i\| \|x^j\|} \right)_{ij} \\
&= \mathrm{diag}\left(\dfrac1{\|x^i\|}\right)_i \left(\langle x^j, x^k \rangle \right)_{jk} \mathrm{diag} \left( \dfrac1{\|x^l\|} \right)_l \\
&= ( \mathrm{diag} \Sigma )^{-1/2} \Sigma ( \mathrm{diag} \Sigma )^{-1/2}.
\end{align*}
ただし、\Sigma は、x^i を横に並べた行列 x = (x^1\ \cdots \ x^N) を用いて
\begin{align*}
\Sigma
&= \dfrac1M \left( \langle x^i, x^j \rangle \right)_{ij} \\
&= \dfrac1{M} \left( (x^i)^t x^j \right)_{ij} \\
&= \dfrac1{M} x^t x
\end{align*}
である。
偏相関係数
相関係数はその名の通り確率変数どうしの相関を表す。しかし、それらは「真の相関」であろうか？二つの確率変数 Y^1 と Y^2 が相関している(つまり、|\rho^{Y^1 Y^2}| が 1 に近い)として、それらは直接に相関しているのだろうか？一般にはNOである。例えば小学生の集団を考える。小学生一般では、身長(Y^1)が高いほど語彙(Y^2)が多いが、これらは明らかに第三の変数である学年(X としよう)に依存している。このとき、身長 Y^1 と語彙数 Y^2 は疑似相関の関係にあると言う。
可能な限り Y^1 と Y^2 の間に発生するような疑似相関を排除したい。そのためには X を固定してしまえば良いのではないのだろうか、という考えがうまれる。 学年 X を固定してしまえば身長 Y^1 と語彙数 Y^2 の間の相関は小さくなるはずである（学年が同じでも誕生月による変化はもちろんあるだろうから、完全に 0 とまではならないと思われるが）。このような相関を Y^1 と Y^2 の偏相関と呼ぶ。
偏相関の定義
記号の定義
確率変数の列 D = (Y^1\ Y^2\ X) = (Y^1\ Y^2\ X^1 \cdots X^N) を考え、確率変数 Y^1, Y^2 の間の偏相関（X = (X^1 \cdots X^N) の観測値を固定した時のY^1 と Y^2 の相関)を r_{Y|X} と書くことにする。また、確率変数 D を M 回測定したとする。第 m 回目の測定値を d_m = (y_m^1\ y_m^2\ x_m^1 \cdots x_m^N) \in \mathbf{R}^{N+2} とする（アルファベットの小文字は測定値を表す。上つき添字と下付き添字はそれぞれ確率変数の番号と測定された順番を示す）。D を M 回測定した結果を縦に並べて行列を作る:
\begin{align*}
d
&=
\begin{pmatrix}
d_1 \\
\vdots \\
d_M
\end{pmatrix} \\
&=
\begin{pmatrix}
y_1^1 & y_1^2 & x_1^1 & \cdots & x_1^N \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
y_M^1 & y_M^2 & x_M^1 & \cdots & x_M^N
\end{pmatrix} \\
&=
\begin{pmatrix}
y^1 & y^2 & x^1 & \cdots & x^N
\end{pmatrix} \\
&=
\begin{pmatrix}
y & x
\end{pmatrix} \\
&=
\begin{pmatrix}
d^1 & d^2 & d^3 & \cdots & d^{N + 2}
\end{pmatrix}
\end{align*}
ここで y^a = (y_1^a \cdots y_M^a)^t, x^n = (x_1^n \cdots x_M^n)^t \in \mathbf{R}^M, y = (y^1\ y^2) \in \mathrm{Mat}(M,2), x = (x^1 \cdots x^N) \in \mathrm{Mat}(M, N), d^i = y^i\ (i = 1,2), d^i = x^{i - 2}\ (i = 3, \ldots, N + 2) とした。
以下では、添え字に a, b が現れれば、それは Y^1, Y^2 等の添え字を表し、n が現れればそれは X^1 \cdots X^N 等の添え字であることを約束しておく。同様に、i や j は d の列の列番号を表すものとする。
また、d による分散共分散行列 \Sigma は
\Sigma = \dfrac1{M} d^t d
であり、相関行列 \rho は
\begin{align*}
\rho
&= (\mathrm{diag} \Sigma)^{-1/2} \Sigma (\mathrm{diag} \Sigma)^{-1/2} \\
&= \left( \dfrac{\langle d^i, d^j \rangle}{\|d^i\| \|d^j\|} \right)_{ij}
\end{align*}
となる。
偏相関の定義
\beta_1, \beta_2, \alpha_1, \ldots, \alpha_N を任意の実数として
\begin{align*}
d
\begin{pmatrix}
\beta_1 \\
\beta_2 \\
\alpha_1 \\
\vdots \\
\alpha_N
\end{pmatrix}
&=
\begin{pmatrix}
y & x
\end{pmatrix}
\begin{pmatrix}
\beta_1 \\
\beta_2 \\
\alpha_1 \\
\vdots \\
\alpha_N
\end{pmatrix} \\
&=
\begin{pmatrix}
y^1 & y^2 & x^1 & \cdots & x^N
\end{pmatrix}
\begin{pmatrix}
\beta_1 \\
\beta_2 \\
\alpha_1 \\
\vdots \\
\alpha_N
\end{pmatrix} \\
&= \beta_1 y^1 + \beta_2 y^2 + \alpha_1 x^1 + \cdots + \alpha_N x^N
\end{align*}
であるから、行列 d の像 \mathrm{Im}{d} \subset \mathbf{R}^M について \mathrm{Im}{d} = \mathrm{Im}{y} \cup \mathrm{Im}{x} が成立し、また一般に \mathrm{Im}{y} \not\perp \mathrm{Im}{x} である。そして、ベクトル y^a は行列y の像 \mathrm{Im}{y} に、 x^n は x の像 \mathrm{Im}{x} に含まれる。
偏相関 r_{Y|X} は X = (X^1 \cdots X^N) を固定した場合の Y = (Y^1, Y^2) の相関であったから、それは確率変数 Y について得られた測定値 y^a から、確率変数 X について得られた測定値 x = (x^n) の影響を排除したデータの相関である。それはどのようなものであろうか？
幾何的な解釈が示唆を与えてくれる。次の図のような N = 1 の場合 (d = (y\ x) = (y^1\ y^2\ x^1)) の例が想像の助けとなるだろう。
この例では M > N+2 = 3 とし、 d はフルランク（つまりランク3）であるとしている。したがって \mathrm{Im}{y} \cap \mathrm{Im}{x} = \{0\}, \mathrm{dim} (\mathrm{Im}{y}) = 2, \mathrm{dim} (\mathrm{Im}{x}) = 1 となっている。また、図では (\mathrm{Im}{x})^\perp の次元が 2 であるように見えているが、実際には M - 2 次元空間であることに注意されたい。さて、緑色のベクトルを y^a \in \mathrm{Im}{y} とすると、y^a は \mathrm{Im}{x} へ直交射影(P とする)したオレンジのベクトル P y^a と、それに垂直な (\mathrm{Im}{x})^\perp へ直交射影(Q とする。Q = 1-P である)した青色のベクトル Q y^a の和である。図からも明らかなように、Q y^a こそが y^a から x の影響を排除したものとなる。
Q によって「影響が排除されて」いることを確認してみよう。同じ記号を用いて Y^1 と Y^2 の間の共分散を計算してみると、
\begin{align*}
\Sigma^{Y^1 Y^2}
&= \dfrac1{M} \langle y^1, y^2 \rangle \\
&= \dfrac1{M} \langle P y^1 + Q y^1, P y^2 + Q y^2 \rangle \\
&= \dfrac1{M} (\langle P y^1, P y^2 \rangle + \langle P y^1, Q y^2 \rangle + \langle Q y^1, P y^2 \rangle + \langle Q y^1, Q y^2 \rangle) \\
&= \dfrac1{M} (\langle P y^1, P y^2 \rangle + \langle Q y^1, Q y^2 \rangle)
\end{align*}
となる。この式から、y^1 と y^2 の間の共分散は \mathrm{Im}{x} に射影された部分どうしの共分散 \langle P y^1, P y^2 \rangle / M とその直交補空間 \mathrm{Im}{x}^\perp に射影された部分どうしの共分散 \langle Q y^1, Q y^2 \rangle / M との和であることがわかる。確かに Q y^a は共分散の意味において「y^a から x の影響を排除して」いる。
以上から、偏相関係数 r_{Y|X} は分散共分散の (\mathrm{Im}{x})^\perp への射影に関する部分 \langle Q y^a, Q y^b \rangle / M を（相関係数の時とおなじように）正規化することで定義できるであろう:
\begin{align*}
r_{Y|X}^{ab}
&= \dfrac{\langle Q y^a, Q y^b \rangle}{\|Q y^a\| \|Q y^b\|} \\
&= \dfrac{\langle y^a, Q y^b \rangle}{\|Q y^a\| \|Q y^b\|}.
\end{align*}
ここで直交射影の性質 Q^t = Q = Q^2 を用いた。
相関行列が正則な時の偏相関係数の計算
以下では相関行列が正則とする。このためには M > N + 2 が必要である。相関行列が正則な場合、偏相関係数行列 R = (r^{ab}) と相関行列 \rho = (\rho^{ij})_{1 \leq i, j \leq N + 2} の間には次のような関係がある; 相関行列の逆行列を \rho^{-1} = (\rho_{ij}), Kroneckerのデルタを \delta_{ab} と書けば
\begin{align*}
r_{Y|X}^{ab} = (-1)^{\delta_{ab} + 1} \dfrac{\rho_{ab}}{\sqrt{\rho_{aa}}\sqrt{\rho_{bb}}}.
\end{align*}
式の意味
証明の前にこの式の意味を考える。相関行列の逆行列は、分散共分散行列 \Sigma を用いて
\begin{align*}
\rho^{-1}
&= (\mathrm{diag} \Sigma)^{1/2} \Sigma^{-1} (\mathrm{diag} \Sigma)^{1/2} \\
&= \left( \|d^i\| \|d^j\| \Sigma_{ij} \right)_{ij}
\end{align*}
と書ける。ここで \Sigma^{-1} = (\Sigma_{ij}) とした。分散共分散行列 \Sigma の逆行列 \Sigma^{-1} を精度行列と呼ぶ。この表現を用いれば偏相関係数と精度行列の関係は
\begin{align*}
r^{ab}
&= (-1)^{\delta_{ab} + 1} \dfrac{\|d^a\| \|d^b\| \Sigma_{ab}}{\sqrt{\|d^a\|^2 \Sigma_{aa}}\sqrt{\|d^b\|^2 \Sigma_{bb}}} \\
&= (-1)^{\delta_{ab} + 1} \dfrac{\Sigma_{ab}}{\sqrt{\Sigma_{aa}} \sqrt{\Sigma_{bb}}}
\end{align*}
となる。この式は精度行列の成分 \Sigma_{ab} が確率変数 Y^a と Y^b の間の直接的な関係を表している事を示している。この解釈は多変量正規分布の場合には直接理解できる。
式の証明
ブロック行列
G =
\begin{pmatrix}
A & B \\
C & D
\end{pmatrix}
を考える。A, D は正方行列、C, D は G が正方行列となるような矩形行列とする。いま、A, D が正則とすれば、 G の逆行列は
G^{-1} =
\begin{pmatrix}
(A - B D^{-1} C)^{-1} & \ast \\
\ast & \ast
\end{pmatrix}
の形をしている。ところで\rho = d^t d / M は
\begin{align*}
\rho
&= \dfrac1M
\begin{pmatrix}
y^t \\
x^t
\end{pmatrix}
\begin{pmatrix}
y & x
\end{pmatrix} \\
&= \dfrac1M
\begin{pmatrix}
y^t y & y^t x \\
x^t y & x^t x
\end{pmatrix}
\end{align*}
である。いま、\rho は正則であるとしているから、A = y^t y, D = x^t x は明らかに正則となり、\rho^{-1} の左上のブロックは
\begin{align*}
\dfrac1M
\begin{pmatrix}
\rho_{11} & \rho_{12} \\
\rho_{21} & \rho_{22}
\end{pmatrix}
&= (y^t y - y^t x (x^t x)^{-1} x^t y)^{-1})^{-1} \\
&= (y^t (1 - x (x^t x)^{-1} x^t) y)^{-1} \\
&= (y^t (1 - P) y)^{-1} \\
&= (y^t Q y)^{-1} \\
&=
\begin{pmatrix}
\langle y^1, Q y^1 \rangle & \langle y^1, Q y^2 \rangle \\
\langle y^2, Q y^1 \rangle & \langle y^2, Q y^2 \rangle
\end{pmatrix}
^{-1} \\
&=
\dfrac1{\langle y^1, Q y^1 \rangle \langle y^2, Q y^2 \rangle - \langle y^1, Q y^2 \rangle \langle y^2, Q y^1 \rangle}
\begin{pmatrix}
\langle y^2, Q y^2 \rangle & - \langle y^1, Q y^2 \rangle \\
- \langle y^2, Q y^1 \rangle & \langle y^1, Q y^1 \rangle
\end{pmatrix}
\end{align*}
となる。x (x^t x)^{-1} x^t = P（\mathrm{Im}{x} への直交射影）であることは x の特異値分解から分かる。退屈な計算の結果、
\begin{align*}
r_{Y|X}^{ab} = (-1)^{\delta_{ab} + 1} \dfrac{\rho_{ab}}{\sqrt{\rho_{aa}}\sqrt{\rho_{bb}}}
\end{align*}
を得る。
Discussion