Bias-Variance decomposition

MikuHatsune2017-01-14

機械学習などで予測モデルを立てたときに絶対に突っ込まれるのが「それ他のデータセットでも言えんの?」ということで、モデルの性能を評価しないといけない。一般的には手持ちのデータからモデルを作って、それとは別のデータに対してそのモデルが有用かという汎化性能をみないといけない。しかしながら、モデルだけ作って丸投げとか、モデルを作ったデータに再度モデルを適応して性能を評価したとか、そういうのがほとんど散見される。
 
Bias-Variance decomposition (日本語)によって、あるモデルに寄る誤差は
\begin{align}E[(Y-\hat{f}(x))^2]&=\sigma^2 + (E_T[\hat{f}(x)]-f(x))^2 + E_T[(\hat{f}(x)-E_T[\hat{f}(x)])^2]\\&=\textrm{irreducible}~\textrm{error} + \textrm{Bias}^2 + \textrm{Variance}\end{align}
とかける。このとき、bias というのはモデルと訓練データから計算できる、モデルの説明具合で、variance というのはも出るとテストデータから計算できる、期待値まわりの分散である。
パラメータ1個のような単純なモデルだと、そもそも手持ちの訓練データを表現することは(たいていの場合)不可能であるが、新しいデータが入ってきても「まあそんなもんだよね」という程度の性能で新しいデータを表現できる。すなわち、bias が大きくてvarianve が小さい、という状況で、underfitting というらしい。
一方で、パラメータが超多いようなモデルでは、手持ちの訓練データをほぼ完璧に表現できる(n個の未知変数はn個の連立方程式が必要というアレ)一方で、手持ちの訓練データに(のみ)ガチガチに当てはまっているから、新しいデータが入力されたときにものすごい外れた値が出力される。すなわち、bias が小さくてvariance が大きい、という状況で、これが巷でよくきくoverfitting である。
 
Bias-Variance decomposition や、それを図示するvalidation curve というのは、機械学習に強い(とよく言われる)Python では実装されているようで、よく見かける。しかし、R では素人でも一発でできるような実装はざっと調べたかぎりではないようで、R-blogger でも2016年の記事があった程度なので、写経してやっておく。
 
いま、多項式で適当に発生させたデータがあるとする。真のモデルは3次式(切片含めて自由度4)で、ランダムエラー\epsilon を入れて
y=X\beta+\epsilon
でデータが発生している。ここで、X は1次元から30次元までの多項式で、x は[0,100] で適当に発生させたものである。\beta多項式の係数。


The Elements of Statistical Learning の図7.1 では、bias とvariance のシュミレーションがあるので、これを再現してみる。原著はPDF があるのでやる気のある人は参照してほしい。

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

統計的学習の基礎 ―データマイニング・推論・予測―

統計的学習の基礎 ―データマイニング・推論・予測―

  • 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman,杉山将,井手剛,神嶌敏弘,栗田多喜夫,前田英作,井尻善久,岩田具治,金森敬文,兼村厚範,烏山昌幸,河原吉伸,木村昭悟,小西嘉典,酒井智弥,鈴木大慈,竹内一郎,玉木徹,出口大輔,冨岡亮太,波部斉,前田新一,持橋大地,山田誠
  • 出版社/メーカー: 共立出版
  • 発売日: 2014/06/25
  • メディア: 単行本
  • この商品を含むブログ (6件) を見る
 
シミュレーションでは50個のデータ点(訓練データ)でモデルを作り、10000個の新規データ(テストデータ)での性能を100回検証している。
灰色は訓練データへのbias である。横軸は多項式の次数の多さで、みぎに行くほど複雑なモデルになる。そのため、訓練データをほとんど完璧に記述できる(エラーが小さくなる)ようになる。これだけみると、モデルの次元は大きく複雑なほうがいいように思える。
赤色はテストデータへの当てはまりである。単純なモデルから次数を増やしてくと、あるところで最小値をとり、その後エラーが増大するようになる。オレンジの点はテストエラーが最小になるときの次数であり、このとき5だったようである。
 
このデータがそもそもどのように出来たかは、我々は(Rコンソール内では)神なので、4次元だと知っている。4次元には近いがちょっと違う。ただし、それは我々が神だから知っているだけの話なので、現実的には最小となればいいだろうと思って採用する。

 
R-blogger のほうではクロスバリデーションを使ったhyperparameter の選択とか書いてある。

# Generate the training and test samples
seed <- 1809
set.seed(seed)
gen_data <- function(n, beta, sigma_eps) {
    eps <- rnorm(n, 0, sigma_eps)
    x <- sort(runif(n, 0, 100))
    X <- cbind(1, poly(x, degree = (length(beta) - 1), raw = TRUE))
    y <- as.numeric(X %*% beta + eps)
    return(data.frame(x = x, y = y))
}

# Fit the models
require(splines)
n_rep <- 100
n_df <- 30
df <- 1:n_df
beta <- c(5, -0.1, 0.004, -3e-05)
n_train <- 50
n_test <- 10000
sigma_eps <- 0.5

xy <- res <- list()
xy_test <- gen_data(n_test, beta, sigma_eps)
for (i in 1:n_rep) {
    xy[[i]] <- gen_data(n_train, beta, sigma_eps)
    x <- xy[[i]][, "x"]
    y <- xy[[i]][, "y"]
    res[[i]] <- apply(t(df), 2, function(degf) lm(y ~ ns(x, df = degf)))
}

# Plot the data
x <- xy[[1]]$x
X <- cbind(1, poly(x, degree = (length(beta) - 1), raw = TRUE))
y <- xy[[1]]$y
plot(y ~ x, col = "gray", lwd = 2, pch=16)
lines(x, X %*% beta, lwd = 3, col = "black")
lines(x, fitted(res[[1]][[1]]), lwd = 3, col = "palegreen3")
lines(x, fitted(res[[1]][[4]]), lwd = 3, col = "darkorange")
lines(x, fitted(res[[1]][[25]]), lwd = 3, col = "steelblue")
legend(x = "topleft", legend = c("True function", "Linear fit (df = 1)", "Best model (df = 4)", 
    "Overfitted model (df = 25)"), lwd = rep(3, 4), col = c("black", "palegreen3", 
    "darkorange", "steelblue"), text.width = 32, cex = 0.85)

# Compute the training and test errors for each model
pred <- list()
mse <- te <- matrix(NA, nrow = n_df, ncol = n_rep)
for (i in 1:n_rep) {
    mse[, i] <- sapply(res[[i]], function(obj) deviance(obj)/nobs(obj))
    pred[[i]] <- mapply(function(obj, degf) predict(obj, data.frame(x = xy_test$x)), res[[i]], df)
    te[, i] <- sapply(as.list(data.frame(pred[[i]])), function(y_hat) mean((xy_test$y-y_hat)^2))
}

# Compute the average training and test errors
av_mse <- rowMeans(mse)
av_te <- rowMeans(te)

# Plot the errors
plot(df, av_mse, type = "l", lwd = 2, col = gray(0.4), ylab = "Prediction error", 
    xlab = "Flexibilty (spline's degrees of freedom [log scaled])", ylim = c(0, 1), log = "x")
abline(h = sigma_eps, lty = 2, lwd = 0.5)
for (i in 1:n_rep) {
  lines(df, te[, i], col = "lightpink")
}
for (i in 1:n_rep) {
  lines(df, mse[, i], col = gray(0.8))
}
lines(df, av_mse, lwd = 2, col = gray(0.4))
lines(df, av_te, lwd = 2, col = "darkred")
points(df[1], av_mse[1], col = "palegreen3", pch = 17, cex = 1.5)
points(df[1], av_te[1], col = "palegreen3", pch = 17, cex = 1.5)
points(df[which.min(av_te)], av_mse[which.min(av_te)], col = "darkorange", pch = 16, cex = 1.5)
points(df[which.min(av_te)], av_te[which.min(av_te)], col = "darkorange", pch = 16, cex = 1.5)
points(df[25], av_mse[25], col = "steelblue", pch = 15, cex = 1.5)
points(df[25], av_te[25], col = "steelblue", pch = 15, cex = 1.5)
legend(x = "top", legend = c("Training error", "Test error"), lwd = rep(2, 2), 
    col = c(gray(0.4), "darkred"), text.width = 0.3, cex = 0.85)

ここで、deviance 関数というのがあるが、これはbias を計算していて、fit という回帰後のオブジェクトがあったとすれば

sum(fit$residuals^2)

もしくは

sum((y - fit$fitted.values)^2)

で求められる。
 
Python ではおなじみだが、入力数が増えると判別性能がどう改善されていくかというlearning curve というのがR ではおなじみではなさそう。caret パッケージでできるようだ。

library(caret)
library(kernlab)
data(spam)
spam0 <- spam
lda_data <- learing_curve_dat(dat = spam, verbose=FALSE,
                              outcome = "type",
                              test_prop = 1/5, 
                              ## `train` arguments:
                              method = "lda", 
                              metric = "ROC",
                              trControl = trainControl(classProbs = TRUE, 
                                                       summaryFunction = twoClassSummary))

v <- split(lda_data, lda_data$Data)
yl <- c(0.9, 1)
plot(0, type="n", xlim=c(1, max(lda_data$Training_Size)), ylim=yl, xlab="Training size", ylab="ROC", las=1)
lines(v$Training$Training_Size, v$Training$ROC, type="o", pch=16, col=2, lwd=3)
lines(v$Testing$Training_Size, v$Testing$ROC, type="o", pch=16, col=4, lwd=3)
legend("bottomright", legend=c("Training", "Testing"), col=c(2, 4), lty=1, pch=16, lwd=3)