二つの確率変数$X, Y$の相関を調べたい場合,相関係数$\rho$を推定するのが一般的です.$n$個のサンプル$(X_1, Y_1), \cdots, (X_n, Y_n)$による相関係数の推定値の一つとして,
$$
\hat{\rho} = \dfrac{\sum_i (X_i - \bar{X}) (Y_i - \bar{Y}) }{\sqrt{\sum_i (X_i - \bar{X})^2 \sum_i (Y_i - \bar{Y})^2 }}
$$
があります(Pearsonの積率相関係数).ただし,算術平均
\bar{X} = \dfrac{1}{n} \sum_i X_i, \quad \bar{Y} = \dfrac{1}{n} \sum_i Y_i
を定義しています.
相関係数の推定値$\hat{\rho}$について調べていると,意外と難しい(そして興味深い)ことに気づいたので,共有します。また説明もなく使われることが多い以下の$t$検定についても,前提条件や導出を示します.
帰無仮説$H_0$:相関係数$\rho = 0$,対立仮説$H_1$:相関係数$\rho \neq 0$という検定を考える.$H_0$のもとで
$$
t = \hat{\rho} \sqrt{ \dfrac{n-2}{1-\hat{\rho}^2} }
$$
は自由度$n-2$のt分布に従うため,$|t| > t_{n-2, \alpha/2}$であるとき帰無仮説$H_0$は棄却される.ただし$t_{n-2,\alpha/2}$は自由度$n-2$のt分布の$\alpha$分位点である.
相関係数の定義とその最尤推定量
説明に入る前に,「パラメータ」とその「推定量」を明確に区別する必要があります.本記事では,確率分布を規定する$\rho$のような量を「パラメータ」と呼び,得られたサンプル$(X_1, Y_1), \cdots, (X_n, Y_n)$による「パラメータ」の推定値を「推定量」と呼びます.「推定量」に関して,本記事では$\hat{\rho}$のようにハット付きで表記します.
設定と問題定義
確率変数$X, Y$について,それぞれの期待値を$(\mu_{X}, \mu_{Y})$,分散を$(\sigma_X^2, \sigma_Y^2)$とします。共分散は
$$
\sigma_{XY} = \mathbb{E}[(X-\mu_X)(Y-\mu_Y)]
$$
として定義されます.このとき,確率変数$X, Y$の相関係数(≠推定値)$\rho$を
$$
\rho = \dfrac{\sigma_{XY}}{\sigma_X \sigma_Y}
$$
として定義します.
ここでの問題は,$\rho$の推定量$\hat{\rho}$を適切に選び,その性質を明らかにすることです.
最尤推定量(MLE)
あるパラメータ$\boldsymbol{\theta}$に基づく確率変数$Z$の分布$f(z \mid \boldsymbol{\theta})$が与えられているものとします.サンプル$(z_1, \cdots, z_n)$が(各々独立に)分布$f(z \mid \boldsymbol{\theta})$から得られたとき,
$$
L(\boldsymbol{\theta}) = \sum_i \log f(z_i \mid \boldsymbol{\theta})
$$
のことを(対数)尤度関数と呼びます.このとき,$\boldsymbol{\theta}$の最尤推定量(Maximum Likelihood Estimator: MLE)は,尤度関数$L(\boldsymbol{\theta})$を最大化する$\hat{\boldsymbol{\theta}}$として定義されます:
$$
\hat{\boldsymbol{\theta}} = {\rm argmax}_{\boldsymbol{\theta}} L(\boldsymbol{\theta})
$$
例えば,以下はそれぞれ$\sigma^2_X, \sigma_Y^2, \sigma_{XY}^2$のMLEとなります(少し計算すればわかる).
\hat{\sigma}^2_X = \dfrac{1}{n} \sum_i (x_i - \bar{x})^{2} \\
\hat{\sigma}^2_Y = \dfrac{1}{n} \sum_i (y_i - \bar{y})^2 \\
\hat{\sigma}_{XY} = \dfrac{1}{n} \sum_i (x_i - \bar{x}) (y_i - \bar{y})
$n$で割ることに関する注意
分散を計算する場合に$n-1$で割ることが標準的であると知っている人にとっては,$n$で割ることに違和感があるかもしれません.分散の推定の際に$n-1$で割っていたのは,推定量の不偏性(推定値の期待値がパラメータに一致するという性質)を保証するためです.一般に,MLEは不偏推定量ではありません.実際, $$ \mathbb{E}[\hat{\sigma}^2_X] = \dfrac{n-1}{n} \sigma^2_X \neq \sigma^2_X $$ であるため,パラメータ$\sigma^2_X$からのずれが生じています.ただし,$n\to \infty$の極限では,MLEは都合がよい性質を持つため,この記事ではMLEを用います.MLEは,(それほど強くない仮定の下で)以下のうれしい性質を持ちます.
$\boldsymbol{\theta}$のMLEを$\hat{\boldsymbol{\theta}}$とすると,
$$
\sqrt{n} (\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}) \to_{d} \mathcal{N}({\bf 0}, \boldsymbol{I}(\boldsymbol{\theta})^{-1})
$$
ただし$\mathcal{N}({\bf 0}, \boldsymbol{I}(\boldsymbol{\theta})^{-1})$は平均${\bf 0}$,共分散行列$\boldsymbol{I}(\boldsymbol{\theta})^{-1}$の多変量正規分布であり,$\boldsymbol{I}(\boldsymbol{\theta})$はFisher情報行列である:
$$
\boldsymbol{I}(\boldsymbol{\theta}) = \mathbb{E}\left[ \dfrac{\partial}{\partial \boldsymbol{\theta}} \log f(Z \mid \boldsymbol{\theta}) \dfrac{\partial}{\partial \boldsymbol{\theta}^t} \log f(Z \mid \boldsymbol{\theta}) \right]
$$
雑に言うと,MLEはサンプル数が十分に多い極限で真の値$\boldsymbol{\theta}$に収束することが保証されます.
さらに,MLEは不変性(≠不偏性)というありがたい性質も持ちます.
パラメータ$\boldsymbol{\tau} = g(\boldsymbol{\theta})$の推定を考えるとき,$\hat{\boldsymbol{\theta}}$が$\boldsymbol{\theta}$のMLEであるならば,$g(\hat{\boldsymbol{\theta}})$が$\boldsymbol{\tau}$のMLEとなる.
これを積率相関係数
$$
\hat{\rho} = \dfrac{\sum_i (X_i - \bar{X}) (Y_i - \bar{Y}) }{\sqrt{\sum_i (X_i - \bar{X})^2 \sum_i (Y_i - \bar{Y})^2 }} = \dfrac{\hat{\sigma}_{XY}}{\hat{\sigma}_X \hat{\sigma}_Y}
$$
に当てはめると,$\hat{\rho}$が相関係数$\rho$のMLEとなることがただちに導かれます.
以上から,積率相関係数$\hat{\rho}$は相関係数$\rho$の適切な推定量(の一つ)であることがわかります.
MLEが従う分布と相関係数に対する検定
積率相関係数$\hat{\rho}$は$X, Y$の分布に依らずMLEとなります.以下では具体的な分布を仮定することで,MLE $\hat{\rho}$の従う確率分布や,仮説検定を導出します.
相関係数のMLEが従う分布
確率変数$X, Y$が相関係数$\rho$の正規分布に従う場合,MLE $\hat{\rho}$の従う分布を計算することができます.共分散行列(をスケール倍したもの)がWishart分布に従うことを用いつつ,適切に変数変換して,ガンマ関数$\Gamma(z)$の性質を用いることで,
f(\hat{\rho} = r \mid \rho) = \dfrac{ (1-\rho^2)^{(n-1)/2} (1-r^2)^{(n-4)/2} }{ \sqrt{\pi} \Gamma \left(\dfrac{n-1}{2} \right)\Gamma\left(\dfrac{n-2}{2} \right)}\sum_{k=0}^{\infty} \dfrac{(2\rho r)^k}{k!} \Gamma \left(\dfrac{n+k-1}{2} \right)^2
が導かれます.この式は$\rho = 0$の場合に限ると,以下のような簡単な形となります.
f(\hat{\rho} = r \mid \rho) \propto (1-r^2)^{(n-4)/2}
いま
$$
t = r \sqrt{\dfrac{n-2}{1-r^2}}
$$
によって変数変換すると,$t$の確率分布として
$$
f(t \mid \rho) \propto \left( 1 + \dfrac{t^2}{n-2} \right)^{-(n-1)/2}
$$
が得られます.これはまさに自由度$n-2$のStudentの$t$分布となっています.
相関係数に関する検定で用いられていた仮定
先ほどの導出に依れば,以下の仮説検定には$X, Y$が正規分布に従うという仮定がなされていたことがわかります.逆に言えば,正規分布にさえ従えば,サンプル数$n$がそれほど多くなくても厳密にStudentの$t$分布によって検定ができることがいえます.
帰無仮説$H_0$:相関係数$\rho = 0$,対立仮説$H_1$:相関係数$\rho \neq 0$という検定を考える.$H_0$のもとで
$$
t = \hat{\rho} \sqrt{ \dfrac{n-2}{1-\hat{\rho}^2} }
$$
は自由度$n-2$のt分布に従うため,$|t| > t_{n-2, \alpha/2}$であるとき帰無仮説$H_0$は棄却される.ただし$t_{n-2,\alpha/2}$は自由度$n-2$のt分布の$\alpha$分位点である.
数値シミュレーション
Wikipediaによると,MLE $\hat{\rho}$は不偏推定量ではないため,(絶対値で見ると)小さめに見積もりがちとのことでした.少し気になったので,簡単なコードで試してみました.
以下は,$\rho = 0.0, 0.5, 0.8$のそれぞれについて,サンプル数$n = 20, 100$で計算したMLEの分布を示したものです.試行回数(発生させた$n$個のサンプルに基づいてMLEを計算する回数)は$N = 4000$回としています.
相関係数は$-1$から$1$の値しかとらない都合上,確かに相関係数の推定値(の絶対値)が上振れることはあまりなさそうです.
$n = 10$から$n = 500$について,$\rho = 0.8$に固定しつつ$\hat{\rho}$の分布を示したものが下図です.青色の線はMLEの算術平均,水色がかかった範囲はMLEの標準偏差をあらわしています.確かに,MLE $\hat{\rho}$は実際のパラメータ$\rho$より下振れる傾向がありそうです.
まとめとその他の相関係数
積率相関係数$\hat{\rho}$が相関係数$\rho$のMLEになるのは,MLEの持つ強い性質(不変性)に因ります.
たかが相関係数ですが,分布はかなり複雑な形状となり,$\rho = 0$の場合に限って単純化されます.
何も考えずに$t$検定することは簡単ですが,その前提($X, Y$が正規分布に従うこと)を理解しておくのも重要です.
またMLEが若干下振れることがあるということも,意外と知られていないことかなと思いました.
相関係数はPearsonの積率相関だけでなく,Spearmanの順位相関など,多くの種類があります.順位相関は確率変数の順位のみに着目するため,外れ値に対しロバストになるという良い性質を持ちます.順位相関の持つ数理的な性質にも興味はありますが,積率相関よりもさらに難しそうな気がしています.
参考文献
久保川達也:「現代数理統計学の基礎」とその補足資料
https://sites.google.com/site/ktatsuya77/xian-dai
※MathStat_hosoku.pdfの式(6)の$(-\rho^2)$は$(1-\rho^2)$の誤植と考えられます.
また久保川さんのpdfではサンプル数ではなく自由度を$n$としているため,式(6)との間にずれがあります.