내삽과 외삽

회귀는 많은 사람들에게 생소하게 느껴질 단어인 내삽(interpolation)과 외삽(extrapolation)으로 구분해 볼 수 있다. 손수회귀의 관점에서 볼 때, 내삽이란 주위에 데이터가 많을 때, 결과값을 예측하는 것이고, 외삽이란 대부분의 데이터와 동떨어진 점에서 결과값을 예측하는 것이라고 생각할 수 있다.

다음의 그림에서 \(x=0\) 일 때, \(\mathbb{E}[y|x=0]\) 을 구하는 것을 내삽이라고 한다면, \(x=4\) 또는 \(x=6\) 일 때, \(\mathbb{E}[y|x]\) 를 구하는 것을 외삽이라고 볼 수 있다.

n = 10000
x <- runif(n, -3, 3)
e <- rnorm(n)
y <- sin(1.4*x)+sqrt(x+4) + log(x+5) + e
dat <- data.frame(x=x, y=y)

plot(y~x, xlim=c(-3, 6))

적합 모형들

다음의 몇 가지 모형에 대해, 위의 데이터로 \(\mathbb{E}[y|x=0]\) , \(\mathbb{E}[y|x=3]\) , \(\mathbb{E}[y|x=4]\) , \(\mathbb{E}[y|x=6]\) 를 추정해 보면 다음과 같다.

• 선형 모형 \(y = b_0 + b_1 x\)
• 3차 다항 모형 \(y = b_0 + b_1 x + b_2 x^2 + b_3 x^3\)
• 랜덤 포레스트
• 그래디언트 부스팅 머신
• 인공신경망

fitLM <- lm(y~x, data=dat)
#fitRF <- randomForest::randomForest(y~x, data=dat)

fitPOLY <- lm(y ~ x + I(x^2) + I(x^3), data=dat)

library(ranger)
fitRF <- ranger(y~x, data=dat, mtry=1, max.depth=3)
# max.depth is manully selected

library(xgboost)
fitXGB <- xgboost(as.matrix(dat[,'x']), label=dat[,"y"], 
                  objective = "reg:linear", verbose=0,
                  nrounds=10)
# nrounds is manully selected

i <- order(x)
library(scales)

op <- par(mar = rep(0, 4)) #https://stackoverflow.com/questions/5663888/trying-to-remove-all-margins-so-that-plot-region-comprises-the-entire-graphic
par(mfcol=c(2,2))
plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(fitLM)[i], col='red', lwd=2)
text(2,0,'lm')

plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(fitPOLY)[i], col='orange', lwd=2)
text(2,0,'polynomial')

plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(fitRF, data=dat[i,])$predictions, col='blue', lwd=2)
text(2,0,'ranger')

plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(fitXGB, newdata=as.matrix(dat[i,'x'])), col='purple', lwd=2)
text(2,0,'xgboost')

#par(op)

par(mfcol=c(1,1))
library(keras)
model <- keras_model_sequential()
model %>%
  layer_dense(units = 32, input_shape = 1) %>%
  layer_activation_leaky_relu() %>%
  layer_dense(units=32) %>%
  layer_activation_leaky_relu() %>%
  layer_dense(units = 1, activation = 'linear')

model %>% compile(
  optimizer = "adam",
  #lr=0.01,
  loss = "mse",
  metrics = "mae"
)

history <-
  model %>% fit(x,
                y,
                batch_size = 100,
                epochs = 20)

par(op)
plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(model, x=x)[i,], col='green', lwd=2)
text(2,0,'keras')

예측값을 모두 모아보면, 선형회귀를 제외하곤, 크게 차이가 없음을 확인할 수 있다.

#par(mar = rep(0, 4))
par(op)
par(mfcol=c(1,1))
plot(y~x, xlim=c(-3, 6), col=alpha('black', 0.1))
lines(x[i], predict(fitLM)[i], col='red', lwd=2)
lines(x[i], predict(fitPOLY)[i], col='orange', lwd=2)
lines(x[i], predict(fitRF, data=dat[i,])$predictions, col='blue', lwd=2)
lines(x[i], predict(fitXGB, newdata=as.matrix(dat[i,'x'])), col='purple', lwd=2)
lines(x[i], predict(model, x=x)[i,], col='green', lwd=2)

이제 외삽을 해보자. 데이터가 없는 \(x=-3\) 에서 \(x=6\) 까지 예측을 해본다.

par(op)
par(mfcol=c(1,1))
plot(y~x, xlim=c(-6, 6), col=alpha('black', 0.1))

newdat = data.frame(x=seq(-6,6,0.1))
lines(newdat$x, predict(fitLM, newdat), col='red', lwd=2)
lines(newdat$x, predict(fitPOLY, newdat), col='orange', lwd=2)
lines(newdat$x, predict(fitRF, data=newdat)$predictions, col='blue', lwd=2)
lines(newdat$x, predict(fitXGB, newdata=as.matrix(newdat[,'x'])), col='purple', lwd=2)
lines(newdat$x, predict(model, x=newdat$x), col='green', lwd=2)

이 엄청난 차이는 무엇 때문일까?

내삽과 외삽의 차이

적당히 복잡한 모형의 경우에는 내삽에서 큰 차이가 나지 않는다. 왜냐하면 주위에 데이터가 많기 때문에 모형은 이들 데이터를 잘 설명해야 할 의무가 있다.

하지만 외삽의 경우는 다르다. 주위에 데이터가 없기 때문에, 멀리에 떨어져 있는 데이터에서 일반화(generalization)을 해야 한다. 그리고 일반화는 모형이 가지고 있는 가정에 의해 좌우된다.

그렇다면 모형이 가지고 있는 가정 또는 태생적 제약은 무엇일까?

• 선형 회귀와 ReLU 활성화 함수 : 선형 회귀에서 \(dy/dx\) (기울기)는 상수로 유지된다. 데이터가 없는 공간에서도 마찬가지이다. (마치 진공의 우주에서 등속운동을 하는 물체와 같다.) ReLU 활성화 함수를 쓴 모형은 piecewise 선형이다. 데이터가 없는 진공 상태에서 \(dy/dx\) 는 상수로 유지되기 마련이다.

• 3차 다항 회귀 : 3차 다항함수는 무엇이 상수인가? 2차 다항 함수는 기울기의 변화율이 상수이다. 3차 다항 함수는 기울기의 변화율의 변화율 \(d^3y/dx^3\) 이 상수이다. 그렇다면 왜 고차 다항함수가 과적합되고, 특히 일반화에 젬병인지 이해할 수 있을 것이다. 세상에 변화율의 변화율의 변화율의 변화율의 …의 변화율이 상수로 유지되는 경우가 거의 없다!

• 랜덤 포레스트와 그래디언트 부스팅 머신 : 수많은 의사 결정 나무로 구성된 이들 모형에서 데이터가 전혀 없는 진공 속에서 변화율 \(dy/dx\) 는 \(0\) 으로 유지된다. 다시 말해 \(\mathbb{E}[y|x]\) 는 상수가 되기 싶다.

마지막으로, 그렇다면 \(\tanh\) 를 활성화 함수로 한다면 어떤 가정을 통해 외삽이 이루어지는가?

par(mfcol=c(1,1))
library(keras)
model <- keras_model_sequential()
model %>%
  layer_dense(units = 32, input_shape = 1) %>%
  layer_activation(activation=activation_tanh) %>%
  layer_dense(units=32) %>%
  layer_activation(activation=activation_tanh) %>%
  layer_dense(units = 1, activation = 'linear')

model %>% compile(
  optimizer = "adam",
  #lr=0.01,
  loss = "mse",
  metrics = "mae"
)

history <-
  model %>% fit(x,
                y,
                batch_size = 100,
                epochs = 20)


par(op)
par(mfcol=c(1,1))
plot(y~x, xlim=c(-6, 6), col=alpha('black', 0.1))

newdat = data.frame(x=seq(-6,6,0.1))
#lines(newdat$x, predict(fitLM, newdat), col='red', lwd=2)
#lines(newdat$x, predict(fitPOLY, newdat), col='orange', lwd=2)
#lines(newdat$x, predict(fitRF, data=newdat)$predictions, col='blue', lwd=2)
#lines(newdat$x, predict(fitXGB, newdata=as.matrix(newdat[,'x'])), col='purple', lwd=2)
lines(newdat$x, predict(model, x=newdat$x), col='green', lwd=2)

text(2,0,'dnn with tanh activation model extrapolation')

작은 데이터

데이터가 작을 때에 모형이 어떻게 적합되는지를 보면, 모형이 각자 나름의 방법으로 추측하고 있음을 알 수 있다.

n = 100
x <- runif(n, -3, 3)
e <- rnorm(n)
y <- sin(1.4*x)+sqrt(x+4) + log(x+5) + e
dat <- data.frame(x=x, y=y)

fitLM <- lm(y~x, data=dat)
#fitRF <- randomForest::randomForest(y~x, data=dat)

fitPOLY <- lm(y ~ x + I(x^2) + I(x^3), data=dat)

fitRF <- ranger(y~x, data=dat, mtry=1, max.depth=3)
# max.depth is manully selected

fitXGB <- xgboost(as.matrix(dat[,'x']), label=dat[,"y"], 
                  objective = "reg:linear", verbose=0,
                  nrounds=10)

model <- keras_model_sequential()
model %>%
  layer_dense(units = 16, input_shape = 1) %>%
  layer_activation_leaky_relu() %>%
  #layer_dense(units=32) %>%
  #layer_activation_leaky_relu() %>%
  layer_dense(units = 1, activation = 'linear')

model %>% compile(
  optimizer = "adam",
  #lr=0.01,
  loss = "mse",
  metrics = "mae"
)

model %>% fit(x,
              y,
              batch_size = 100,
              epochs = 400)

par(mfcol=c(2,1))
plot(y~x, xlim=c(-6, 6))
plot(y~x, xlim=c(-6, 6))
newdat = data.frame(x=seq(-6,6,0.1))
lines(newdat$x, predict(fitLM, newdat), col='red', lwd=2)
lines(newdat$x, predict(fitPOLY, newdat), col='orange', lwd=2)
lines(newdat$x, predict(fitRF, data=newdat)$predictions, col='blue', lwd=2)
lines(newdat$x, predict(fitXGB, newdata=as.matrix(newdat[,'x'])), col='purple', lwd=2)
lines(newdat$x, predict(model, x=newdat$x), col='green', lwd=2)