二つの確率変数$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)との間にずれがあります．