標本相関係数の幾何学的解釈（ピアソンの積率相関係数）

p次元の標本が得られたときの標本相関係数の幾何学的解釈について解説する。

図などを交えながら、相関係数を三角比（cos）で表現できることをみていく。

また、2次元正規母集団からの標本が得られたときに、各種最尤推定量を計算し相関係数の最尤推定量である標本相関係数を計算する。

標本相関係数の幾何学的解釈

標本\((\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N) = \boldsymbol{X}\)が得られたとする。\(\boldsymbol{X}\)は次で与えられる\(p\times N\)行列である。

\begin{align}\label{eq1} \boldsymbol{X} = \begin{pmatrix}x_{11} & \cdots & x_{1N}\\ \vdots & &\vdots\\ x_{p1} & \cdots & x_{pN} \end{pmatrix} = \begin{pmatrix}\boldsymbol{u}_1^T\\ \vdots\\ \boldsymbol{u}_p^T. \tag{1}\end{pmatrix}\end{align}

ここに、\(\boldsymbol{u}_i^T\)は\(\boldsymbol{X}\)の\(i\)番目の行である。ベクトル\(\boldsymbol{u}_i\)は片方の端点が\(x_{i\alpha},\ \alpha = 1, \ldots, N\)であり、もう一方の端点は原点であることが分かる。すなわち\(\alpha\)に関する座標をもつ\(N\)次元空間上のベクトルであると考えられる。このように、標本\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N\)を\(N\)次元ユークリッド空間上の\(p\)個のベクトルで表現することができる。ユークリッド距離の定義より、\(\boldsymbol{u}_i\)の長さの2乗は

\begin{align}d(\boldsymbol{u}_i, \boldsymbol{u}_i)= \|\boldsymbol{u}_i\|^2 = \boldsymbol{u}_i^T\boldsymbol{u}_i = \sum_{\alpha=1}^N x_{i\alpha}^2\end{align}

である。

i番目とj番目の変量の\(\cos\)

\(\boldsymbol{u}_i\)と\(\boldsymbol{u}_j\)の間の角\(\theta_{ij}\)に関する\(\cos\)は

\begin{align}\label{eq2} \cos\theta_{ij} &= \cfrac{\boldsymbol{u}_i^T\boldsymbol{u}_j}{\sqrt{\boldsymbol{u}_i^T\boldsymbol{u}_i \boldsymbol{u}_j^T\boldsymbol{u}_j}} = \cfrac{\sum_{\alpha=1}^N x_{i\alpha}x_{j\alpha}}{\sqrt{\sum_{\alpha = 1}^N x_{i\alpha}^2 \sum_{\alpha=1}^N x_{j\alpha}^2}}\tag{2}\end{align}

である。

ベクトル\(d\boldsymbol{u}_j\)が\(\boldsymbol{u}_i -d \boldsymbol{u}_j\)と直交するようにスカラー\(d\)を決める。すなわち\(d\)は次のように与えられる。

\begin{align} &(d\boldsymbol{u}_j)^T(\boldsymbol{u}_i - d\boldsymbol{u}_j) = 0 \\ &\Leftrightarrow d(\boldsymbol{u}_j^T\boldsymbol{u}_i - d\boldsymbol{u}_j^T\boldsymbol{u}_j) = 0\\ &\Leftrightarrow d\boldsymbol{u}_j^T\boldsymbol{u}_j = \boldsymbol{u}_j^T\boldsymbol{u}_i \\ &\Leftrightarrow d = \cfrac{\boldsymbol{u}_j^T\boldsymbol{u}_i}{\boldsymbol{u}_j^T\boldsymbol{u}_j}\end{align}

次の図1に示すように、\(\boldsymbol{u}_i\)を\(\boldsymbol{u}_i - d\boldsymbol{u}_j\)と\(d\boldsymbol{u}_j\)に分解する。\(\boldsymbol{u}_i\)と\(\boldsymbol{u}_j\)の間の角\(\theta_{ij}\)の\(\cos\)の絶対値は\(d\boldsymbol{u}_j\)の長さを\(\boldsymbol{u}_i\)の長さで割ったものである。

図1

すなわち

\begin{align}\cos\theta_{ij} &=\sqrt{\cfrac{(d\boldsymbol{u}_j)^T(d\boldsymbol{u}_j)}{\boldsymbol{u}_i^T\boldsymbol{u}_i}}\\&= \sqrt{\cfrac{\boldsymbol{u}_i^T\boldsymbol{u}_id^2}{\boldsymbol{u}_i^T \boldsymbol{u}_i }}\\&= \sqrt{\cfrac{\boldsymbol{u}_j^T\boldsymbol{u}_j(\boldsymbol{u}_j^T\boldsymbol{u}_i / \boldsymbol{u}_j^T\boldsymbol{u}_j)^2}{\boldsymbol{u}_i^T\boldsymbol{u}_i}}\\&= \cfrac{\boldsymbol{u}_i^T\boldsymbol{u}_j}{\boldsymbol{u}_i^T\boldsymbol{u}_i\boldsymbol{u}_j^T\boldsymbol{u}_j}.\end{align}

よって、\eqref{eq2}が示された。□

i番目とj番目の変量の\(\cos\)の関係を用いて標本相関係数の幾何学的解釈を行う。\(\boldsymbol{A} = \sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T\)としたとき、\(a_{ij} / \sqrt{a_{ii}a_{jj}}\)の幾何学的解釈をするためにequiangular lineを導入する。これは原点と点\((1, 1, \ldots, 1)\)を通るベクトルである。ベクトル\(\boldsymbol{\varepsilon} = (1, 1, \ldots, 1)^T\)上での\(\boldsymbol{u}_i\)の射影は

\begin{align}\cfrac{\boldsymbol{\varepsilon}^T\boldsymbol{u}_i}{\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}}\boldsymbol{\varepsilon} &= \cfrac{\sum_{\alpha=1}^N 1\cdot x_{i\alpha}}{\sum_{\alpha=1}^N 1\cdot 1}\boldsymbol{\varepsilon}\\&= \cfrac{\sum_{\alpha=1}^Nx_{i\alpha}}{N}\boldsymbol{\varepsilon} \\&= \bar{x}_i\boldsymbol{\varepsilon}\\&= (\bar{x}_i, \bar{x}_i, \ldots, \bar{x}_i)\end{align}

となる。さらに、\(\boldsymbol{u}_i\)を\(\boldsymbol{u}_i\)のequiangular lineへの射影\(\bar{x}_i\boldsymbol{\varepsilon}\)とequiangular lineに垂直なベクトル\(\boldsymbol{u}_i - \bar{x}_i \boldsymbol{\varepsilon}\)に分解する。\(\boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon}\)の長さの2乗は

\begin{align} d(\boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon}, \boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon}) &= \|\boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon}\|^2\\&= (\boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon})^T(\boldsymbol{u}_i-\bar{x}_i\boldsymbol{\varepsilon})\\&= \sum_{\alpha=1}^N(x_{i\alpha} - \bar{x}_i)^2\end{align}

であり、これは\(a_{ii}\)である。\(\boldsymbol{u}_i - \bar{x}_i\boldsymbol{\varepsilon}\)と\(\boldsymbol{u}_j - \bar{x}_j\boldsymbol{\varepsilon}\)の端点が原点に来るように、それぞれのベクトルを移動させる。\(\boldsymbol{u}_i - \bar{x}_i\boldsymbol{\varepsilon}\)と\(\boldsymbol{u}_j - \bar{x}_j\boldsymbol{\varepsilon}\)の間の角\(\theta_{ij}\)の\(\cos\)は

\begin{align}\cos\theta_{ij} &= \cfrac{(\boldsymbol{u}_i- \bar{x}_i\boldsymbol{\varepsilon})^T(\boldsymbol{u}_j - \bar{x}_j\boldsymbol{\varepsilon})}{\sqrt{(\boldsymbol{u}_i - \bar{x}_i\boldsymbol{\varepsilon} )^T(\boldsymbol{u}_i - \bar{x}_i\boldsymbol{\varepsilon})(\boldsymbol{u}_j - \bar{x}_j\boldsymbol{\varepsilon})^T(\boldsymbol{u}_j - \bar{x}_j\boldsymbol{\varepsilon})}} \\&= \cfrac{\sum_{\alpha = 1}^N(x_{i\alpha} - \bar{x}_i)(x_{j\alpha} - \bar{x}_j)}{\sqrt{\sum_{\alpha = 1}^N(x_{i\alpha} - \bar{x}_i)^2 \sum_{\alpha = 1}^N(x_{j\alpha} - \bar{x}_j)^2}}\end{align}

であり、これは標本相関係数である。

このように標本相関係数はequiangular lineに垂直なベクトル\(\boldsymbol{u}_i - \bar{x}_i\boldsymbol{\varepsilon}\)と\(\boldsymbol{u}_j- \bar{x}_j\boldsymbol{\varepsilon}\)の間の角に関する\(\cos\)となることが示された（図2の原点に描かれている角）。

図2

標本相関係数の計算

標本相関係数を実際に計算する。今、2次元正規母集団から次の標本を得たとする。\begin{array}{ccc}\hline Sample &x_1&x_2 \\\hline 1 & 0.00002 & 0.78444 \\ 2 & 1.92039 & 0.33999\\ 3 & 3.81333 & 0.89105\\ 4 & 2.39655 & 0.71903 \\ 5 & 6.23519 & 2.77218 \\ 6 & -0.77875 & -0.36046\\ 7 & 3.11125 & 1.39780\\ 8 & -0.68311 & -0.22957\\ 9& -3.58486& -1.43994\\ 10& 5.09494 & 1.60923\end{array}

2標本正規分布の標本

このとき、標本平均ベクトル\(\hat{\boldsymbol{\mu}}\)、標本共分散行列\(\hat{\boldsymbol{\Sigma}}\)、不偏標本共分散行列\(\boldsymbol{S}\)、標本\(x_{11}, \ldots, x_{1,10}\)と\(x_{21},\ldots x_{2,10}\)の標本相関係数\(\rho_{12}\)は次となる。

\begin{align}\hat{\mu} &= \bar{x} = \begin{pmatrix}1.75249\\ 0.64837\end{pmatrix}, \\ \hat{\Sigma} &=\begin{pmatrix}8.17006 & 2.98437\\ 2.98437 & 1.23225 \end{pmatrix}, \\ S &= \begin{pmatrix}9.07785 & 3.31506 \\ 3.31506 & 1.36917 \end{pmatrix},\\ \hat{\rho}_{12} &= 0.94057.\end{align}

上の推定量の計算には次のようにR言語を用いた。

library(MASS)

N <- 10
mu <- c(2, 1)
Sigma <- rbind(c(6, 3), c(3, 2))

samples <- mvrnorm(N, mu,Sigma)
plot(samples, xlim = c(-2, 5), ylim = c(-2, 5), xlab = "x1", ylab = "x2")

x_bar <- apply(samples, 2, mean)
Sigma_hat <- (N - 1) / N * var(samples)
S <- var(samples)
rho_12 <- cor(samples[, 1], samples[, 2])

x_bar
Sigma_hat
S
rho_12