「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。

MENU

統計的機械学習の数理100問(05/20)

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

を用いることにする。

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固有値1I-H固有値0の固有空間の次元はp+1,\ H固有値0I-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}XI-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-1t分布に従うはずである。
 回帰モデル


\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}であることから、これらは独立である。

プライバシーポリシー お問い合わせ