いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として
を用いることにする。
3. リサンプリング
- 真の統計モデルよりも複雑なモデルを選択してしまう過学習による影響を受けずに学習の性能を評価する方法であるクロスバリデーションを学ぶ。
- 学習結果のバラつきを評価する汎用的な方法であるブートストラップを学ぶ。
3.1 クロスバリデーション
組のデータがあるときにそのうちの一部を推定結果を評価するのに用いることにする方法は、推定のための標本が少なくなり推定精度が劣化するというジレンマが発生する。これに対抗するための方法が(
-fold)クロスバリデーション(
)である。
を
を満たすように取った上で、評価値を以下のような手順で推定を行う：
| (1) | 全データを |
|
| (2) | ||
| (3) | 番号 |
|
| (4) | 番号 |
|
| (5) | ||
| (6) | 各 |
3.1.1 クロスバリデーションの実施例
重回帰モデルにおいてとし説明変数
のうち
および切片に依存すると仮定して、乱数で生成したテストデータを用いて
-
クロスバリデーションを行なった。
############################ ### クロスバリデーション ### ############################ 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)
線形回帰の母数推定における分割数別
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) | ロジスティック回帰の変数選択で |
|
| (2) | 線形判別および二次判別とで |
|
| (3) |
3.2.1 適用例
に
を変えつつ
近傍法を適用して誤り率を評価する。
################################### ### 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)
