Situation
「自由度」の悩み
統計学を勉強していると、自由度(degree of freedom)という概念に苦戦している方が多いのではないでしょうか。
例えば、
「そもそもの自由度の定義は?」
「カイ二乗検定で自由度を求める公式を覚えたけど、理論的な根拠がわからない...」
「回帰分析の残差平方和の自由度がn-pになる理由は?」
こうした疑問に対し、多くの教材では「とにかく公式を覚えましょう」という暗記頼みのアプローチや、数式ではなく直感的な説明で覚えることが取られがちです。
しかし、これでは根本的な理解には至らず、統計検定などの試験対策としても不十分と考えます。
本記事で解決したい課題
本記事では、暗記に頼らない理論的な理解を提供します。
特に回帰分析の決定係数を具体例として、RSSとTSSの自由度がそれぞれ$n-p$,$n-1$となる理由を、べき等行列の性質等を用いて厳密に導出します。自由度関連の問題に対して、確固たる理論的基盤を持ってアプローチできるようになることを目指します。
Complication
自由度の定義
統計量に対しての自由度の定義について、明確に数式として定義されているものは見つけられませんでした(もしご存知の方がいれば教えてください)。
ただ、自由度は「カイ二乗分布の自由度」や「t分布の自由度」など、特定の確率分布に関連して使用されることが多く、
統計量に対して「自由度 xxx の統計量」という言葉で用いられている場合は、暗黙に特定の分布が想定されていると考えられます。
自由度の計算の具体例
決定定数に登場する、各統計量の自由度とその証明
最小二乗法における回帰分析の文脈のモデルを考えます(ここでは真の分布が正規分布であると仮定します)。
\begin{align}
y = X\beta + \epsilon \\
y \in \mathbb{R}^n, & X \in \mathbb{R}^{n \times p}, \beta \in \mathbb{R}^p, \epsilon \sim N(0, \sigma^2 I_n) \\
X = [\mathbf{1}\ X_1\ X_2\ \ldots\ X_{p-1}] & (\text{ただし、} \mathbf{1} \text{は } n \text{ 次元の全ての要素が } 1 \text{ のベクトル}) \\
\end{align}
ここで、真の分布が正規分布として仮定できる場合において、モデルパラメータ数 $p$ を決める方法として、AICやBICなどの情報量規準もありますが、決定係数ならびに自由度調整済み決定係数 (adjusted R-squared) が利用されることがあるそうです。
統計検定準一級の試験範囲にもある内容であり、決定係数は以下のように定義されます。
\begin{align}
R^2 = 1 - \frac{RSS}{TSS}
\end{align}
このとき、残差平方和 (RSS: Residual Sum of Squares) は以下のように定義されます。
\begin{align}
RSS &= \| y - X\hat{\beta} \|^2 \\
\hat{\beta} &= (X^\top X)^{-1}X^\top y
\end{align}
また、全平方和 (TSS: Total Sum of Squares) は以下のように定義されます。
\begin{align}
TSS = \| y - \bar{y} \|^2 \\
\bar{y} = \frac{1}{n} \sum_{i=1}^n y_i
\end{align}
ここで、$RSS$ と $TSS$ の従う分布を考えます。
RSS が従う分布
\begin{align}
RSS &= \| y - X\hat{\beta} \|^2 \\
&= (y - X(X^\top X)^{-1}X^\top y)^\top (y - X(X^\top X)^{-1}X^\top y) \\
&= (y - Hy)^\top (y - Hy) & (\because H := X(X^\top X)^{-1}X^\top ) \\
&= y^\top (I - H)^\top (I - H)y \\
\end{align}
ここで $I-H$が べき等行列であることを利用できる。というのも下記の式が成立するため、$I-H$ はべき等行列となります:
\begin{align}
(I - H)(I - H) &= I - 2H + H^2 \\
&= I - 2H + X(X^\top X)^{-1}X^\top X(X^\top X)^{-1}X^\top \\
&= I - 2H + X I (X^\top X)^{-1}X^\top \\
&= I - 2H + H \\
&= I - H
\end{align}
よって、
\begin{align}
RSS &= y^\top (I - H)y
\end{align}
ここで、べき等行列の固有値は0か1であるため、べき等行列 $I-H$ のランクは、$I-H$ のトレースに等しいことを利用します。
\begin{align}
\text{rank}(I - H) &= \text{trace}(I - H) \\
&= \text{trace}(I) - \text{trace}(H) \\
&= n - \text{trace}(X(X^\top X)^{-1}X^\top ) \\
&= n - \text{trace}((X^\top X)^{-1}X^\top X) \\
&= n - \text{trace}(I) \\
&= n - p
\end{align}
ここで、$I-H$に対して固有値分解を行うと、
\begin{align}
I - H = P \Lambda P^\top
\end{align}
ここで、$I-H$は対称行列のため、ある直交行列$P$を用いて対角化可能で、その結果を$\Lambda$としました。$\Lambda$ の対角成分は $I-H$ の固有値です。すなわち、$\Lambda$ の対角成分は0か1であり、1の個数は $\text{rank}(I-H) = n - p$ であることから、下記が成立します。
\begin{align}
RSS &= y^\top P \Lambda P^\top y \\
&= (P^\top y)^\top \Lambda (P^\top y) \\
&= \sum_{i=1}^{n} \lambda_i (P^\top y)_i^2 \\
&= \sum_{i=1}^{n-p} (P^\top y)_i^2 & (\because \lambda_i = 1 \text{ for } i=1,\ldots,n-p \text{ and } \lambda_i = 0 \text{ for } i=n-p+1,\ldots,n) \\
\end{align}
したがって、$P^\top y$ は $y \sim N(X \beta, \sigma^2 I_n)$ に対して直交変換を受けたものであるため、$P^\top y \sim N(P^\top X \beta, \sigma^2 I_n)$ が成立します。
PTXβ の平均の計算
ここで、平均を考えます。
\begin{align}
P^\top X \beta &= P^\top I X \beta \\
&= P^\top (H + (I - H)) X \beta \\
&= P^\top H X \beta + P^\top (I - H) X \beta
\end{align}
まず、右辺の第二項目を考えます。
\begin{align}
(I-H)X\beta
&= X\beta - H X\beta \\
&= X\beta - X(X^\top X)^{-1} X^\top X\beta \\
&= X\beta - X\beta \\
&= 0
\end{align}
よって、 $P^\top (I - H) X \beta = 0$ となります。
次に、右辺の第一項目を考えます。
$(I-H)$ の固有値分解は下記のように表される。
\begin{align}
I-H = P \Lambda P^\top , \quad \Lambda = \begin{bmatrix} I_{n-p} & 0 \\ 0 & 0_p \end{bmatrix}, \quad P^\top P = I
\end{align}
ここで $P = [P_1\ P_2]$ と分割し、$P_1 \in \mathbb{R}^{n \times (n-p)}, P_2 \in \mathbb{R}^{n \times p}$ とする。すると、
- $P_1$ は固有値 1 の固有ベクトルを列に持つ
- $P_2$ は固有値 0 の固有ベクトルを列に持つ
固有値 0 の定義より:
\begin{align}
(I-H) P_2 = 0 \quad\Rightarrow\quad P_2^\top (I-H) P_2 = 0
\end{align}
また $H = I - (I-H)$ なので:
\begin{align}
P_1^\top H = P_1^\top (I - (I-H)) = P_1^\top - P_1^\top (I-H) = P_1^\top - P_1^\top I = 0
\end{align}
よって
\begin{align}
P_1^\top X \beta = 0
\end{align}
となります。
以上のことから、$P^\top y$の平均は$0$ となります。
よって、$RSS$ は以下のような分布に従います。
\begin{align}
\frac{RSS}{\sigma^2} = \sum_{i=1}^{n-p} \left(\frac{(P^\top y)_i}{\sigma}\right)^2 \sim \chi^2(n - p)
\end{align}
TSS の従う分布
\begin{align}
TSS &= \| y - \bar{y} \|^2 \\
&= (y - \bar{y}\mathbf{1})^\top (y - \bar{y}\mathbf{1}) \\
&= y^\top (I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top ) y
\end{align}
ここで、$I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top$ がべき等行列であるため、同様にrankを計算します。
\begin{align}
\text{rank}(I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top ) &= \text{trace}(I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top ) \\
&= \text{trace}(I) - \text{trace}(\frac{1}{n}\mathbf{1}\mathbf{1}^\top ) \\
&= n - \frac{1}{n} \text{trace}(\mathbf{1}\mathbf{1}^\top ) \\
&= n - \frac{1}{n} \cdot n \\
&= n - 1
\end{align}
ここで、$I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top$に対して固有値分解を行うと、
\begin{align}
I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top = Q \Gamma Q^\top
\end{align}
ここで、$I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top$は対称行列のため、ある直交行列$Q$を用いて対角化可能で、その結果を$\Gamma$としました。$\Gamma$ の対角成分は $I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top$ の固有値です。すなわち、$\Gamma$ の対角成分は0か1であり、1の個数は $\text{rank}(I - \frac{1}{n}\mathbf{1}\mathbf{1}^\top) = n - 1$ であることから、下記が成立します。
\begin{align}
TSS &= y^\top Q \Gamma Q^\top y \\
&= (Q^\top y)^\top \Gamma (Q^\top y) \\
&= \sum_{i=1}^{n} \gamma_i (Q^\top y)_i^2 \\
&= \sum_{i=1}^{n-1} (Q^\top y)_i^2 & (\because \gamma_i = 1 \text{ for } i=1,\ldots,n-1 \text{ and } \gamma_i = 0 \text{ for } i=n) \\
\end{align}
よって、$Q^\top y \sim N(Q^\top X \beta, \sigma^2 I_n)$ です。
ここで、$Q^\top X \beta$ の平均は一般的に$0$ではないことに注意します。なぜなら、$X$ の第一列は $\mathbf{1}$ であるため、$Q^\top \mathbf{1} \neq 0$ となるためです。そのため、$TSS$ は非心カイ二乗分布に従います。
\begin{align}
\frac{TSS}{\sigma^2} = \sum_{i=1}^{n-1} \left(\frac{(Q^\top y)_i}{\sigma}\right)^2 \sim \chi^2(n - 1, \lambda),\quad \lambda = \sum_{i=1}^{n-1} \left(\frac{(Q^\top X \beta)_i}{\sigma}\right)^2
\end{align}
自由度調整済み決定係数
上記の結果を用いて、自由度調整済み決定係数 (adjusted R-squared) は以下のように定義されます。
\begin{align}
\bar{R}^2 = 1 - \frac{RSS/(n - p)}{TSS/(n - 1)}
\end{align}
ここで、$RSS/(n - p)$ は $\sigma^2$ の不偏推定量であり、$TSS/(n - 1)$ は $y$ の分散の不偏推定量であることがわかります。したがって、$\bar{R}^2$ は $R^2$ よりも過学習に対して頑健な指標となります。
Summary
統計検定準1級における自由度の理解において、最も重要なのは自由度 = 行列のランクという関係性だと考えます。
本資料では、回帰分析の決定係数を例に、以下の具体的な計算結果を示しました:
-
RSS(残差平方和): $\frac{RSS}{\sigma^2} \sim \chi^2(n-p)$
- 自由度 = $\text{rank}(I-H) = n-p$(サンプル数 - パラメータ数)
-
TSS(全平方和): $\frac{TSS}{\sigma^2} \sim \chi^2(n-1, \lambda)$
- 自由度 = $\text{rank}(I-\frac{1}{n}\mathbf{1}\mathbf{1}^T) = n-1$(サンプル数 - 1)
これらの自由度は、べき等行列の性質($\text{rank} = \text{trace}$)を用いて導出できます。
実践的な意義: 自由度調整済み決定係数 $\bar{R}^2 = 1 - \frac{RSS/(n-p)}{TSS/(n-1)}$ では、各統計量を適切な自由度で割ることで、不偏推定量を得て過学習を防いでいます。
このような行列計算による自由度の導出過程を理解することで、完全な暗記に頼らず、理論的な裏付けを持って自由度の概念を把握でき、納得を持ちながら統計学の勉強を進められると考えます。
余談
よく、ネットなどを見ると自由度の説明について「自由に動けるデータの数」などと説明されることがありますが、これは誤解を招く表現だと思います。というのも、この説明は具体的にどのような状況を指しているのかが曖昧であり、数学的に厳密ではないからです。
加えて、自由に動けるという表現自体が誤解を招く可能性があります。ひとまず統計量の値を固定して、その統計量を構成する確率変数のうち固定する必要でない数と解釈するとします。これについても問題があると考えます。例えば、今回の$TSS$などにおいて,$y$のうちの$n-1$個のデータが限りなく大きい値を保持するように動かしたとします。このとき、残りの一つのデータは虚数などの非現実的な値を取らない限りは統計量の値が変わってしまいます。つまり、$n-1$個のデータが自由に動けるからといって、必ずしも統計量の値が自由に動けるわけではありません(もし、虚数などの非現実的な値を許容するならば、統計量の値は自由に動けることになりますが、いずれにせよ、曖昧な表現であることには変わりありません)。
どんな分野でも共通なことがいえると思うですが、曖昧な表現を避け、できるだけ具体的かつ数学的に厳密な表現を用いることが、理解を深める上で重要だと考えます。用語を独自に使う場合は、説明の前にその用語の定義を明示することも重要です。