統計的機械学習の数理100問（05/20） - 「大人の教養・知識・気付き」を伸ばすブログ

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
1. 線形回帰
次回

前回

power-of-awareness.com

1. 線形回帰

1.5　 $\hat{\boldsymbol{\beta}}$ の分布

1.5.2　決定係数

　 $L=\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\|$ に $\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}$ を代入したものの分布を導くことにする。
　その前に $H=\boldsymbol{X}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}\in\mathbb{R}^{N\times N}$ の性質を考える。

$\begin{aligned} H^2&=\boldsymbol{X}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}\cdot \boldsymbol{X}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}=H,\\ (I-H)^2&=I-2H+H^2=I-H,\\ H\boldsymbol{X}&=\boldsymbol{X}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}\cdot \boldsymbol{X}=\boldsymbol{X} \end{aligned}$

である。更に $\hat{\boldsymbol{Y}}=\boldsymbol{X}\hat{\boldsymbol{\beta}}$ とおくとき、

$\begin{aligned} \hat{\boldsymbol{Y}}&=\boldsymbol{X}\hat{\boldsymbol{\beta}}=\boldsymbol{X}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}\boldsymbol{Y}=H\boldsymbol{Y},\\ \boldsymbol{Y}-\hat{\boldsymbol{Y}}&=(I-H)\boldsymbol{Y}=(I-H)(\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon})=(\boldsymbol{X}-H\boldsymbol{X})\boldsymbol{\beta}+(I-H)\boldsymbol{\varepsilon}=(I-H)\boldsymbol{\varepsilon},\\ RSS&=\|\boldsymbol{Y}-\hat{\boldsymbol{Y}}\|={}^{t}\{(I-H)\boldsymbol{\varepsilon}\}(I-H)\boldsymbol{\varepsilon}={}^{t}\boldsymbol{\varepsilon}(I-H)^2\boldsymbol{\varepsilon}={}^{t}\boldsymbol{\varepsilon}(I-H)\boldsymbol{\varepsilon} \end{aligned}$

である。

$H$ の性質　 $H,\ I-H$ の固有値は $0,1$ のみであって、 $H$ の固有値 $1$ と $I-H$ の固有値 $0$ の固有空間の次元は $p+1,\$ $H$ の固有値 $0$ と $I-H$ の固有値 $1$ の固有空間の次元は $N-p-1$ となる。

（ $\because$ 　 $rank(X)=p+1$ であるから、

$\begin{aligned} rank(X({}^{t}XX)^{-1}{}^{t}X)\leq \min\{rank(X({}^{t}X X)^{-1}, rank(X))\}\leq p+1 \end{aligned}$

が成り立つ。
　また $X({}^{t}XX)^{-1}{}^{t}XX=X$ についても同様に

$\begin{aligned} rank(X({}^{t}XX)^{-1}{}^{t}X)\geq rank(X({}^{t}XX)^{-1}{}^{t}XX)=rank(X)=p+1 \end{aligned}$

が成り立つ。したがって $rank(X({}^{t}XX)^{-1}{}^{t})=p+1$ が成立する。
　次に $X({}^{t}XX)^{-1}{}^{t}XX$ であるから、 $X$ の各列が $X({}^{t}XX)^{-1}{}^{t}X$ の像の基底を成し、 $X({}^{t}XX)^{-1}{}^{t}X$ の固有値 $1$ の固有ベクトルである。 $X({}^{t}XX)^{-1}{}^{t}X$ の像の次元が $p+1$ であるから、その核の次元は $N-p-1$ である。
　更に任意の $\boldsymbol{x}\in\mathbb{R}^{p+1}$ について

$\begin{aligned} (I-X({}^{t}XX)^{-1}{}^{t}X)\boldsymbol{x}=0\Leftrightarrow X({}^{t}XX)^{-1}{}^{t}X\boldsymbol{x}=\boldsymbol{x} \end{aligned}$

および

$\begin{aligned} (I-X({}^{t}XX)^{-1}{}^{t}X)\boldsymbol{x}=\boldsymbol{x}\Leftrightarrow X({}^{t}XX)^{-1}{}^{t}X\boldsymbol{x}=0 \end{aligned}$

が成り立つため、 $X({}^{t}XX)^{-1}{}^{t}X$ と $I-X({}^{t}XX)^{-1}{}^{t}X$ で固有値 $0,1$ の固有空間が入れ替わる。　 $\blacksquare$ ）

1.5.3　RSSの分布（続き）

　 $I-X({}^{t}XX)^{-1}{}^{t}X$ について

$\begin{aligned} {}^{t}(I-X({}^{t}XX)^{-1}{}^{t}X)(I-X({}^{t}XX)^{-1}{}^{t}X)&=(I-X\{{}^{t}X X)^{-1} {}^{t}X\}^2\\ &=I-2X({}^{t}XX)^{-1} {}^{t}X +\{X\{{}^{t}X X)^{-1}{}^{t}X\}^2\\ &=I-X({}^{t}X X)^{-1}{}^{t}X \end{aligned}$

と実対称行列である。したがってある直交行列 $P$ が存在しそれを用いて対角行列

$\begin{aligned} P(I-X({}^{t}XX)^{-1}{}^{t}X) \end{aligned}$

を構築できる。その固有値のうち $N-p-1$ 個が $1$ でそれ以外の $p+1$ 個が $0$ であるから、対角成分の最初の $N-p-1$ 個を $1$ にすることができる。
　また $\boldsymbol{v}=P\boldsymbol{\varepsilon}\in\mathbb{R}^{N}$ とすると、 $\boldsymbol{\varepsilon}={}^{t}P\boldsymbol{v}$ であるから、 $\boldsymbol{v}=(v_1,\cdots,v_N)$ とおけば

$\begin{aligned} L&={}^{t}\boldsymbol{\varepsilon}\{I-X({}^{t}XX)^{-1}{}^{t}X\}\boldsymbol{\varepsilon}\\ &={}^{t}({}^{t}P\boldsymbol{v})(I-X({}^{t}X X)^{-1}{}^{t}X){}^{t}P\boldsymbol{v}\\ &=\displaystyle{\sum_{i=1}^{N-p-1}{v_i}^2} \end{aligned}$

が成り立つ。また $\boldsymbol{v}$ の最初の $N-p-1$ 個の成分のみからなるベクトルを $\boldsymbol{w}$ とおけば、

$\begin{aligned} E[\boldsymbol{w}]=\boldsymbol{0}\ (\because E[{}^{t}P\boldsymbol{v}]=\boldsymbol{0}) \end{aligned}$

である。また $\tilde{I}$ を対角成分の最初の $N-p-1$ 個が $1$ 、残り $p+1$ 個が $0$ であるような対角行列として

$\begin{aligned} E[\boldsymbol{v}{}^{t}\boldsymbol{v}]&=E[P\boldsymbol{\varepsilon}{}^{t}(P\boldsymbol{\varepsilon})]\\ &=E[P\boldsymbol{\varepsilon}{}^{t}\boldsymbol{\varepsilon}{}^{t}P]\\ &=P \sigma^2 \tilde{I}{}^{t}P=\sigma^2\tilde{I} \end{aligned}$

であるから、 $E[\boldsymbol{w}{}^{t}\boldsymbol{w}]=\sigma^2I$ である。
　したがって

$\begin{aligned} \displaystyle{\frac{\|\boldsymbol{y}-\hat{\boldsymbol{y}}\|^2}{\sigma^2}}\sim \chi_{N-p-1}^2 \end{aligned}$

が成り立つ。

1.5.4　カイ二乗分布の確率密度

1.6　 $\hat{\beta}_j=0$ の仮説検定

　 $\hat{\beta}_j=0,\ j=1,2,\cdots,p$ の仮説検定、すなわち

$\begin{aligned} &\ \ \ \ \ \ H_0^j:\hat{\beta}_j=0\\ v.s.&\ \ \ \ \ \ H_1^j:\hat{\beta}_j\neq0 \end{aligned}$

を考える。
　 $\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\in\mathbb{R}^{p}$ および $\boldsymbol{\beta}\in\mathbb{R}^{p+1}$ が固定されているとして

$\begin{aligned} \boldsymbol{y}_1=\beta_0+\boldsymbol{x}_1{}^{t}[\beta_1,\cdots,\beta_p]+\varepsilon_1,\cdots,\boldsymbol{y}_n=\beta_0+\boldsymbol{x}_n{}^{t}[\beta_1,\cdots,\beta_p]+\varepsilon_n, \end{aligned}$

が得られる。

####################
### 推定のテスト ###
####################

library("ggplot2")

# テスト回数の設定
num_test <- 1000

# テスト結果を格納
df_results <- data.frame(matrix(NA, nrow = num_test, ncol = 2))
colnames(df_results) <- c("beta_0","beta_1")

# テストの試行
for(i in 1:num_test){
  # 乱数シードの設定
  set.seed(i)
  
  # 標本数の設定
  num_samples <- 1000
  
  # 観測値の生成
  vc_x_obs <- rnorm(num_samples, mean = 2, sd = 1)
  vc_res <- rnorm(num_samples, mean = 0, sd = 1)
  
  # 母数の真の値
  num_beta_0 <- 1
  num_beta_1 <- 1
  
  # y = 1 + x + e というモデルを考える
  vc_y <- num_beta_0 + num_beta_0 * vc_x_obs + vc_res
  
  # モデルの指定
  frml_model <- formula(vc_y~vc_x_obs)
  
  # 推定
  lst_lm <- lm(formula = frml_model,data = data.frame(vc_x_obs,vc_y))
  
  df_results[i,] <- lst_lm$coefficients
}

# 
g <- ggplot(data = df_results,aes(x = beta_0, y = beta_1)) + geom_point()
g <- g + theme_bw(base_family = "")  + labs(x = expression(beta[0]),y = expression(beta[1]),color = "")
g <- g + theme(plot.title = element_text(hjust = 0.5),legend.position = "none",
               legend.title=element_text(size = 7),
               legend.text=element_text(size = 6))
g <- g + xlim(num_beta_0*c(0.7,1.3)) + ylim(num_beta_1*c(0.7,1.3))

g

$\beta_0,\ \beta_1$ のプロット(*上記の実行結果)

1.6.1　 $t$ 分布の導出

　有意水準 $\alpha$ を設定したとして、統計量

$\begin{aligned} T=\displaystyle{\frac{U}{\sqrt{\displaystyle{\frac{V}{m}}}}} \end{aligned}$

を考える。ここで $U\sim N(0,1),\ V\sim \chi_m^2$ とする。
　もし $T$ が両端の上位下位 $\alpha/2$ の確率に含まれるのであれば、帰無仮説 $H_0:\beta_j=0$ を棄却する。 $\beta_j=0$ ならば理論的に $T\sim t_{N-p-1}$ 、すなわち自由度 $N-p-1$ の $t$ 分布に従うはずである。
　回帰モデル

$\begin{aligned} Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon \end{aligned}$

における誤差項 $\varepsilon$ の標準偏差を

$\begin{aligned} \hat{\sigma}=\sqrt{\displaystyle{\frac{RSS}{N-p-1}}} \end{aligned}$

で推定し、 $\hat{\beta}_j$ の標準偏差を

$\begin{aligned} SE(\hat{\beta}_j)=\hat{\sigma}\sqrt{B_j} \end{aligned}$

とする。ここで $B_j$ は $\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}$ の第 $j$ 対角成分である。
　このような状況設定で統計量

$\begin{aligned} t=\displaystyle{\frac{\hat{\beta}_j-\beta}{SE(\hat{\beta}_j)}}\sim t_{N-p-1} \end{aligned}$

を示そう。まず

$\begin{aligned} t=\displaystyle{\frac{\hat{\beta}_j-\beta}{SE(\hat{\beta}_j)}}&=\displaystyle{\frac{\hat{\beta}_j-\beta}{\hat{\sigma}\sqrt{B_j}}}\\ &=\displaystyle{\frac{\hat{\beta}_j-\beta}{\sqrt{B_j}}}\sqrt{\displaystyle{\frac{N-p-1}{RSS}}}\\ &=\displaystyle{\frac{\hat{\beta}_j-\beta}{\sigma\sqrt{B_j}}}\sqrt{\displaystyle{\frac{N-p-1}{\displaystyle{\frac{RSS}{\sigma^2}}}}}\\ \end{aligned}$

であるから、

$\begin{aligned} U&:=\displaystyle{\frac{\hat{\beta}_j-\beta}{\sigma\sqrt{B_j}}}\sim N(0,1)\\ V&:=\displaystyle{\frac{RSS}{\sigma^2}}\sim \chi_{N-p-1}^2 \end{aligned}$

が成立する。
　あとは $U,V$ が独立であることを示す。そのために $RSS$ が $\boldsymbol{Y}-\hat{\boldsymbol{Y}}$ の関数であることから、 $\boldsymbol{Y}-\hat{\boldsymbol{Y}}$ および $\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}$ の独立性を示せば充分である。

$\begin{aligned} (\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}(\boldsymbol{Y}-\hat{\boldsymbol{Y}})=({}^{t}\boldsymbol{X}\boldsymbol{X})^{-1}{}^{t}\boldsymbol{X}\boldsymbol{\varepsilon}{}^{t}\boldsymbol{\varepsilon}(I-H) \end{aligned}$

に注意すれば、 $\mathbb{E}[\boldsymbol{\varepsilon}{}^{t}\boldsymbol{\varepsilon}]=\sigma^2 I$ および $H\boldsymbol{X}=\boldsymbol{X}$ より

$\begin{aligned} \mathbb{E}[(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}){}^{t}(y-\hat{y})]=\boldsymbol{0} \end{aligned}$

が成立する。 $\boldsymbol{Y}-\hat{\boldsymbol{Y}}=(I-H)\boldsymbol{\varepsilon},\ \hat{\boldsymbol{\beta}}-\boldsymbol{\beta}$ ともに正規分布に従うから、両者の共分散が $\boldsymbol{0}$ であることから、これらは独立である。