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

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

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

作者:讓, 鈴木
共立出版

Amazon

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

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
1. 線形回帰
- 1.7.　決定係数と共線性の検出
次回

前回

power-of-awareness.com

1. 線形回帰

1.7.　決定係数と共線性の検出

　以降、

$\begin{aligned} W&=\begin{bmatrix} \displaystyle{\frac{1}{n}}&\displaystyle{\frac{1}{n}}&\cdots&\displaystyle{\frac{1}{n}}\\ \vdots&\ddots& &\vdots\\ \vdots& &\ddots&\vdots\\ \displaystyle{\frac{1}{n}}&\cdots&\displaystyle{\frac{1}{n}}&\displaystyle{\frac{1}{n}} \end{bmatrix}\in\mathbb{R}^{n\times n}\\ W\boldsymbol{Y}&=\begin{bmatrix} \bar{Y}&\cdots&\cdots&\bar{Y}\\ \vdots&\ddots& &\vdots\\ \vdots& &\ddots&\vdots\\ \bar{Y}&\cdots&\bar{Y}&\bar{Y} \end{bmatrix},\bar{Y}=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}Y_i} \end{aligned}$

とする。
　観測ベクトル $\boldsymbol{X}_1,\cdots,\boldsymbol{X}_n\in\mathbb{R}^n,Y_1,\cdots,Y_n\in\mathbb{R}$ に対して

$\begin{aligned} 残差変動:\ RSS&=\|\hat{Y}-\boldsymbol{Y}\|=\|(I-H)\boldsymbol{\varepsilon}\|^2=\|(I-H)\boldsymbol{Y}\|^2,\\ 回帰変動:\ ESS&=\|\hat{Y}-\bar{Y}\|^2=\|\hat{Y}-W\boldsymbol{Y}\|^2=\|(H-W)\boldsymbol{Y}\|^2,\\ 全変動:\ TSS&=\|\boldsymbol{Y}-\bar{Y}\|^2=\|(I-W)\boldsymbol{Y}\|^2 \end{aligned}$

を定義する。 $TSS$ に比して $RSS$ が小さければそのモデルのデータへの適合度が高いと言える。
　これらの間には

$\begin{aligned} TSS=RSS+ESS \end{aligned}$

が成り立つ。ここで $H\boldsymbol{X}=\boldsymbol{X}$ であり、 $\boldsymbol{X}$ の最も左の列はすべて $1$ である。またすべての成分が $1$ であるようなベクトルおよびそれの定数倍のベクトルも固有値 $1$ の固有ベクトルになり、 $HW=W$ が成立する。したがって

$\begin{aligned} (I-H)(H-W)=O \end{aligned}$

も成り立つ。
　次にベクトル $(I-W)\boldsymbol{Y}=(I-H)\boldsymbol{Y}+(H-W)\boldsymbol{Y}$ の両辺で大きさの二乗を取ることで

$\begin{aligned} \{(I-W)\boldsymbol{Y}\}^2=&\{(I-H)\boldsymbol{Y}+(H-W)\boldsymbol{Y}\}^2\\ =&\{(I-H)\boldsymbol{Y}\}^2+{}^t\boldsymbol{Y}{}^t(I-H)(H-W)\boldsymbol{Y}\\ &+{}^{t}\boldsymbol{Y}{}^{t}(H-W)(I-H)\boldsymbol{Y}+\{(H-W)\boldsymbol{Y}\}^2\\ =&\{(I-H)\boldsymbol{Y}\}^2+{}^t\boldsymbol{Y}(I-H)(H-W)\boldsymbol{Y}\\ &+{}^{t}\boldsymbol{Y}(H-W)(I-H)\boldsymbol{Y}+\{(H-W)\boldsymbol{Y}\}^2\\ =&\{(I-H)\boldsymbol{Y}\}^2+\{(H-W)\boldsymbol{Y}\}^2 \end{aligned}$

が成り立つ。したがって $\|(I-W)\boldsymbol{Y}\|=\|(I-H)\boldsymbol{Y}\|+\|(H-W)\boldsymbol{Y}\|$ である。
　さらに $RSS$ と $ESS$ が独立であることも示すことができる。実際

$\begin{aligned} \mathbb{V}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{Y}]&=\mathbb{V}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{X}\boldsymbol{\beta}+(H-W)\boldsymbol{\varepsilon}]\\ &=\mathrm{Cov}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{X}\boldsymbol{\beta}]+ \mathrm{Cov}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{\varepsilon}]\\ &=\mathrm{Cov}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{\varepsilon}] \end{aligned}$

である。そして $(I-H)(H-W)=O$ に注意すれば

$\begin{aligned} \mathrm{Cov}[(I-H)\boldsymbol{\varepsilon},(H-W)\boldsymbol{\varepsilon}]=\mathbb{E}[(I-H)\boldsymbol{\varepsilon}{}^{t}\boldsymbol{\varepsilon}(H-W)]=O \end{aligned}$

である。 $RSS,ESS$ ともに正規分布に従うから、したがって両者は独立である。
　指標

$\begin{aligned} R^2=\displaystyle{\frac{RSS}{TSS}}=1-\displaystyle{\frac{ESS}{TSS}} \end{aligned}$

を決定係数という。単回帰( $p=1$ )のときには決定係数は

$\begin{aligned} \hat{\rho}=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}}{\sqrt{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}\displaystyle{\sum_{i=1}^{n}(Y_i-\bar{Y})^2}}}} \end{aligned}$

の二乗に等しい。この意味で決定係数は説明変数と目的変数の相関の大きさを表していると考えられる。実際、 $\hat{Y}=\hat{\beta}_0+\hat{\beta}_1 X$ および $\hat{\beta}_0=\bar{Y}-\hat{\beta}_1\bar{X}$ より

$\begin{aligned} \hat{Y}&=(\bar{Y}-\hat{\beta}_1\bar{X})+\hat{\beta}_1 X\\ \Leftrightarrow\ \hat{Y}-\bar{Y}&=\hat{\beta}_1(X-\bar{X}) \end{aligned}$

が成り立つ。したがって

$\begin{aligned} \hat{\beta}_1=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}}{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}}} \end{aligned}$

および

$\begin{aligned} \|\boldsymbol{X}-\bar{\boldsymbol{X}}\|^2&=\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}\\ \|\boldsymbol{Y}-\bar{\boldsymbol{Y}}\|^2&=\displaystyle{\sum_{i=1}^{n}(Y_i-\bar{Y})^2} \end{aligned}$

より

$\begin{aligned} \displaystyle{\frac{ESS}{TSS}}&=\displaystyle{\frac{{\hat{\beta}_1}^2\|\boldsymbol{X}-\bar{\boldsymbol{X}}\|^2}{\|\boldsymbol{Y}-\bar{\boldsymbol{Y}}\|^2}}\\ &=\left\{\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}}{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}}}\right\}^2\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}}{\displaystyle{\sum_{i=1}^{n}(Y_i-\bar{Y})^2}}}\\ &=\displaystyle{\frac{\left\{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}\right\}^2}{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}\displaystyle{\sum_{i=1}^{n}(Y_i-\bar{Y})^2}}}=\hat{\rho}^2 \end{aligned}$

を得る。
　また調整済み決定係数という、分子の $RSS$ と分母の $ESS$ をそれぞれ自由度で割ったものに置き換えた

$\begin{aligned} 1-\displaystyle{\frac{\displaystyle{\frac{RSS}{N-p-1}}}{\displaystyle{\frac{TSS}{N-1}}}} \end{aligned}$

が用いられることもある。これは $p$ が大きくなればなるほど普通の決定係数が大きくなることを踏まえ、変数が大きくなると決定係数が小さくなるような罰則を課した決定係数に定義し直したものである。
　決定係数は、説明変数によって目的変数がどれだけ説明できているかを表す。他方で説明変数同士に冗長なものが無いかを測る尺度として $\mathrm{VIF}$ がある。 $j$ 番目の説明変数の $\mathrm{VIF}_j$ は

$\begin{aligned} \mathrm{VIF}_j=\displaystyle{\frac{1}{1-R_{X_j|X_{-j}}^2}} \end{aligned}$

である。ここで $R_{X_j|X_{-j}}^2$ は目的変数として $\boldsymbol{X}\in\mathbb{R}^{n\times p}$ の第 $j$ 列を、説明変数としてそれ以外の変数を用いたときの重回帰における決定係数である。 $\mathrm{VIF}$ はその説明変数の他の変数に対する共線性を検出するものである。

1.7.1　決定係数のシミュレーション結果(R)

##########################
### 決定係数を定義する ###
##########################

# 決定係数を与える関数
R2 <- function(vc_x,vc_y){
  vc_y_hat <- lm(formula(vc_y~vc_x))$fitted.values # y_hat
  vc_y_bar <- mean(vc_y) # y_bar
  
  # TSS, RSSを定義
  nm_RSS <- sum((vc_y - vc_y_hat)^2)
  nm_TSS <- sum((vc_y - vc_y_bar)^2)
  
  return(1-nm_RSS/nm_TSS)
}

# 調整済み決定係数
R2_adj <- function(vc_x,vc_y){
  vc_y_hat <- lm(formula(vc_y~vc_x))$fitted.values # y_hat
  vc_y_bar <- mean(vc_y) # y_bar
  
  # TSS, RSSを定義
  nm_RSS <- sum((vc_y - vc_y_hat)^2)
  nm_TSS <- sum((vc_y - vc_y_bar)^2)
  
  # 自由度を調整する
  nm_TSS <- nm_TSS/(nrow(vc_x)-1)
  nm_RSS <- nm_RSS/(nrow(vc_x)-ncol(vc_x)-1)

  # 返り値
  return(1-nm_RSS/nm_TSS)
}

# テストデータを作成
nm_sim <- 10000 # シミュレーション数
nm_n <- 500 # 標本数
nm_m <- 1 # 説明変数の数

# テスト結果を格納
df_result_1 <- data.frame(matrix(nrow = 3, ncol = nm_sim))

rownames(df_result_1) <- c("R^2","rho^2","AE")
colnames(df_result_1) <- 1:nm_sim

for(i in 1:nm_sim){
  set.seed(i)
  
  vc_x <- matrix(rnorm(nm_n * nm_m), ncol = nm_m)
  vc_y <- rnorm(nm_n)
  
  # 決定係数を計算
  df_result_1[,i] <- c(R2(vc_x,vc_y),
                     as.numeric(cor(vc_x,vc_y)^2),
                     R2(vc_x,vc_y) - as.numeric(cor(vc_x,vc_y)^2))
}

ret <- c(apply(df_result_1,1,mean),apply(df_result_1,1,sd))


# 説明変数を増やすとどうなるか調べてみる
nm_max_var <- 500 # 説明変数の数の上限

vc_result <- c()
vc_result_adj <- c()

for(nm_m in 1:nm_max_var){
  vc_x <- matrix(rnorm(nm_n * nm_m), ncol = nm_m)
  vc_y <- rnorm(nm_n)
  
  # 決定係数を計算
  vc_result <- c(vc_result,R2(vc_x,vc_y))
  vc_result_adj <- c(vc_result_adj,R2_adj(vc_x,vc_y))
}

テスト結果

1.7.2　決定係数のシミュレーション結果(Python)

import numpy as np

# 決定係数R^2を定義する
def R2(x,y):
    n = x.shape[0]
    xx = np.insert(x, 0, 1, axis = 1)
    beta = np.linalg.inv(xx.T@xx)@xx.T@y
    
    y_hat = xx@beta
    y_bar = np.mean(y)
    
    RSS = np.linalg.norm(y - y_hat)**2
    TSS = np.linalg.norm(y - y_bar)**2
    
    return 1 - RSS/TSS

# シミュレーションを実施
n = 300
m = 1

x = np.random.randn(n,m)
y = np.random.randn(n)

# 決定係数
print(R2(x,y))

# 相関係数の二乗
xx = x.reshape(n)
print(np.corrcoef(xx,y)[0,1]**2)

# 絶対誤差
print(R2(x,y) - np.corrcoef(xx,y)[0,1]**2)

1.7.3　VIFのシミュレーション(R)

#####################
### VIFを計算する ###
#####################

R2 <- function(vc_x,vc_y){
  vc_y_hat <- lm(formula(vc_y~vc_x))$fitted.values # y_hat
  vc_y_bar <- mean(vc_y) # y_bar
  
  # RSS, TSSを定義
  nm_RSS <- sum((vc_y - vc_y_hat)^2)
  nm_TSS <- sum((vc_y - vc_y_bar)^2)
  
  return(1-nm_RSS/nm_TSS)
}

VIF <- function(vc_x){
  p <- ncol(vc_x)
  vc_values <- array(dim = p)
  
  for(j in 1:p){
    vc_values[j] <- 1/(1 - R2(vc_x[,-j],vc_x[,j]))
  }
  
  return(vc_values)
}

library("MASS")

x <- as.matrix(Boston)
VIF(x)

1.7.4　VIFのシミュレーション(Python)

import numpy as np
from sklearn.datasets import load_boston

# Bostonを読み込む
boston = load_boston()
x = boston.data

print(x.shape) # 行数(標本数)と列数(説明変数の数)を出す

# 決定係数R^2を定義する
def R2(x,y):
    n = x.shape[0]
    xx = np.insert(x, 0, 1, axis = 1)
    beta = np.linalg.inv(xx.T@xx)@xx.T@y
    
    y_hat = xx@beta
    y_bar = np.mean(y)
    
    RSS = np.linalg.norm(y - y_hat)**2
    TSS = np.linalg.norm(y - y_bar)**2
    
    return 1 - RSS/TSS

# VIFを定義する
def VIF(x):
    p = x.shape[1] # 列数を取得
    values = list() # リストを定義

    for j in range(p):
        S = list(set(range(p)) - {j}) # 説明変数のセットのうちjのみ除く
        values.append(1/(1 - R2(x[:,S],x[:,j])))
    return values

VIF(x)

次回

power-of-awareness.com

前回

1. 線形回帰

1.7. 決定係数と共線性の検出

1.7.1 決定係数のシミュレーション結果(R)

1.7.2 決定係数のシミュレーション結果(Python)

1.7.3 VIFのシミュレーション(R)

1.7.4 VIFのシミュレーション(Python)

次回

1.7.　決定係数と共線性の検出

1.7.1　決定係数のシミュレーション結果(R)

1.7.2　決定係数のシミュレーション結果(Python)

1.7.3　VIFのシミュレーション(R)

1.7.4　VIFのシミュレーション(Python)