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

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

MENU

統計的機械学習の数理100問(11/X)

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

を用いることにする。

3. リサンプリング

  • 真の統計モデルよりも複雑なモデルを選択してしまう過学習による影響を受けずに学習の性能を評価する方法であるクロスバリデーションを学ぶ。
  • 学習結果のバラつきを評価する汎用的な方法であるブートストラップを学ぶ。

3.1 クロスバリデーション

 N組のデータがあるときにそのうちの一部を推定結果を評価するのに用いることにする方法は、推定のための標本が少なくなり推定精度が劣化するというジレンマが発生する。これに対抗するための方法が(k-fold)クロスバリデーション(\mathrm{CV})である。
 k\in\mathbb{N}N\equiv0(\mathrm{mod}\ k)を満たすように取った上で、評価値を以下のような手順で推定を行う:

   (1) 全データをk分割し、各データセットに番号i=1,2,\cdots,kを振る。
   (2) i=1とする。
   (3) 番号i以外のデータセットを用いて母数\thetaの推定値\hat{\theta}_{(i)}を計算する。
   (4) 番号iのデータセットおよび推定値\hat{\theta}_{(i)}を用いて結果を評価する。
   (5) i=kならば(6)に進み、そうでなければi=i+1として(3)に戻る。
   (6) i=1,2,\cdots,kについて(4)で得た結果の算術平均を以て最終的な評価値とする。
3.1.1 クロスバリデーションの実施例

 重回帰モデルにおいてp=5とし説明変数X_1,\cdots,X_5のうちX_3,X_4,X_5および切片に依存すると仮定して、乱数で生成したテストデータを用いてk-\mathrm{fold}クロスバリデーションを行なった。

############################
### クロスバリデーション ###
############################

library("ggplot2")

# クロスバリデーションによる回帰係数の推定誤差の評価
cv_linear <- function(mt_x, vc_y, k){
  if(length(vc_y)%%k==0){
    n <- length(vc_y)
    m <- n/k
    S <- 0
    
    for(j in 1:k){
      vc_test <- ((j - 1) * m + 1):(j * m) # テスト対象の行番号列
      
      vc_beta <- solve(t(mt_x[-vc_test,]) %*% mt_x[-vc_test,]) %*% t(mt_x[-vc_test,]) %*% vc_y[-vc_test]
      vc_residuals <- vc_y[vc_test] - mt_x[vc_test,] %*% vc_beta
      S <- S + drop(t(vc_residuals) %*% vc_residuals)
    }
    
    return(S/n)
  }else{
    return(NA)
  }
}

##################################
### テストデータを人工的に作成 ###
##################################

n <- 100
p <- 5

mt_x <- matrix(rnorm(n * p), ncol = p) 
mt_x <- cbind(1,mt_x) # 回帰切片に相当する部分を追加
vc_beta <- rnorm(p + 1) # 推定対象となる回帰係数
vc_beta[c(2,3)] <- 0 # beta_2,beta_3=0とする
vc_eps <- rnorm(n) # 誤差項

vc_y <- mt_x %*% vc_beta + vc_eps

# 評価を実施

cv_linear(mt_x[,c(1,4:6)], vc_y, 10)
cv_linear(mt_x[,c(1:4)], vc_y, 10)
cv_linear(mt_x, vc_y, 10)

###
n <- 100
p <- 5
mt_x <- matrix(rnorm(n * p), ncol = p)
mt_x <- cbind(1,mt_x)
vc_beta <- rnorm(p+1)
vc_beta[c(2,3)] <- 0
vc_U <- NULL
vc_V <- NULL

for(j in 1:n){
  vc_eps <- rnorm(n)
  vc_y <- mt_x %*% vc_beta +vc_eps
  vc_U <- c(vc_U, cv_linear(mt_x[,c(1,4:6)], vc_y, 10))
  vc_V <- c(vc_V, cv_linear(mt_x,vc_y,10))
}

df_cv <- data.frame("var_456" = vc_U,"var_total" = vc_V)

g <- ggplot(df_cv,aes(x = var_456,y = var_total)) 
g <- g + geom_point() + geom_abline(slope = 1,color = "red",
                                    intercept = 0)
g <- g + theme_light() + theme(plot.title = element_text(hjust = 0.5),legend.position = "bottom",
                               legend.title=element_text(size = 10),
                               legend.text=element_text(size = 10))
g <- g + ggtitle("全変数を選択した方が予測誤差が大きい")
g <- g + xlab("変数4,5,6を選択した際の二乗誤差") + ylab("全変数を選択した際の二乗誤差")
g <- g + xlim(c(0.6,1.6)) + ylim(c(0.6,1.6))
plot(g)

# CV
nm_test <- 30

df_test_results <- data.frame()

for(i in 1:nm_test){
  set.seed(i)
  vc_test_no <- c()
  vc_results <- c()
  
  mt_x <- matrix(rnorm(n * p), ncol = p)
  mt_x <- cbind(1,mt_x)
  vc_beta <- rnorm(p+1)
  vc_beta[c(2,3)] <- 0
  
  vc_eps <- rnorm(n)
  
  vc_y <- mt_x %*% vc_beta + vc_eps
  
  for(j in 2:n){
    if(n%%j==0){
      vc_test_no <- c(vc_test_no,j)
      vc_results <- c(vc_results,cv_linear(mt_x, vc_y, j))
    }
  }
  
  df_test_results <- rbind(df_test_results,data.frame("test_no" =rep(i,length(vc_test_no)), "valid_no" = vc_test_no, "results" = vc_results))
}

df_test_results[,1] <- factor(df_test_results[,1])

g <- ggplot(df_test_results,aes(x=valid_no,y=results,color=test_no)) + geom_line() + geom_point()
g <- g + theme_light() + theme(plot.title = element_text(hjust = 0.5),legend.position = "none",
                                        legend.title=element_text(size = 10),
                                        legend.text=element_text(size = 10))
g <- g + ggtitle(paste0("(テスト回数:",nm_test ,"回)"))
g <- g + xlab("k") + ylab("CV")
g <- g + xlim(c(0,100))
plot(g)


線形回帰の母数推定における分割数別\mathrm{CV}予測誤差の推移

import numpy as np
from numpy.random import randn

def cv_linear(X,y,K):
    n = len(y)
    m = int(n/K)
    S = 0
    for j in range(K):
        test = list(range(j * m,(j+1) * m)) # テストデータの添字を取る
        train = list(set(range(n))-set(test)) # 訓練データの添字を取る
        beta = np.linalg.inv(X[train,].T@X[train,])@X[train,].T@y[train]
        e = y[test]-X[test,]@beta
        S = S + np.linalg.norm(e)**2
    return S/n

n = 100
p = 5
X = randn(n, p)
X = np.insert(X,0,1,axis=1)
beta = randn(p+1)
beta[[1,2]] = 0
y = X@beta + randn(n)
cv_linear(X[:,[0,3,4,5]],y,10)

3.2 クロスバリデーションの適用範囲

 クロスバリデーションは線形回帰における変数選択のみならず、さまざまなものに応用できる:

   (1) ロジスティック回帰の変数選択で\mathrm{CV}による誤り率の評価値を比較して最適な変数の組み合わせを求める。
   (2) 線形判別および二次判別とで\mathrm{CV}による誤り率の評価値を比較してより良い方を選択する。
   (3) K近傍法において\mathrm{CV}で各Kにおける誤り率を比較して最適なKを決定する。
3.2.1 適用例

 \mathrm{iris}Kを変えつつK近傍法を適用して誤り率を評価する。

###################################
### CVによるK近傍法の誤り率評価 ###
###################################

# K=1,...,10に対して誤り率をグラフ化

library(ggplot2)

knn_1 <- function(x,y,z,k){
  x <- as.matrix(x)
  n <- nrow(x); p <- ncol(x)
  dis <- array(dim = n)
  
  for(i in 1:n)dis[i] <- norm(z-x[i,],"2")
  
  S <-order(dis)[1:k]
  u <- sort(table(y[S]),decreasing = T)
  
  v <- names(u)[1]
  
  if(length(u)>1) u <- names(u)[2]
  while(length(u)>1 && v==w){
    k <- k-1
    u <- u[k]
    v <- names(u)[1]
    w <- names(u)[2]
  }
  return(v)
}

knn <- function(x,y,z,k){
  n <- nrow(z)
  w <- array(dim = n)
  for(i in 1:n)w[i] <- knn_1(x,y,z[i,],k)
  return(w)
}

########################
### 以下、テスト実施 ###
########################
data(iris)

#
df <- iris[sample(1:150,size = 150,replace = F),]

n <- nrow(df)
U <- NULL; V <- NULL
K <- 10

for(k in 1:K){
  top_seq <- 1 + seq(0, 135, 10)
  S <- 0
  
  for(top in top_seq){
    idx <- top:(top+14)
    knn_ans <- knn(x = df[-idx,1:4],y = df[-idx,5],z = df[idx,1:4],k = k)
    ans <- df[idx,5]
    
    S <- S + sum(knn_ans!=ans)
  }
  
  S <- S/n
  U <- c(U,k)
  V <- c(V,S)
}

df_error_ratio <- data.frame("category" = rep("CV",length(U)),"K"=factor(U),"err"=V * 100)

# 図示
g <- ggplot(df_error_ratio,aes(x = K, y = err, group = category)) + geom_point() + geom_line()
g <- g + ggtitle("CVによる誤り率の評価")
g <- g + xlab("K") + ylab("謝り率(%)")
g <- g + ylim(0,8)
g <- g + theme_classic() + theme(plot.title = element_text(hjust = 0.5),legend.position = "none",
                                 legend.title=element_text(size = 7),
                                 legend.text=element_text(size = 7))
g <- g + geom_hline(yintercept = 2, linetype =3)
g <- g + geom_hline(yintercept = 4, linetype =3)
g <- g + geom_hline(yintercept = 6, linetype =3)
g <- g + geom_hline(yintercept = 8, linetype =3)
plot(g)
プライバシーポリシー お問い合わせ