ANN: 평균과 표준편차 출력 모형
keras와 인공신경망(ANN; Artificial Neural Network)
인공신경망이 대단해 보일지라도, 그 핵심을 따지고 들어간다면, 다음의 함수를 추정하는 것이라고 볼 수 있다.
\[f(x) = \mathbb{E}[y|x]\]
이는 손수회귀(Manual Regression)에서도 설명한 것처럼, 직관적으로 생각해도 그리 어렵지 않은 개념이다. 하지만 고차원의 설명변수 \(\vec{x}\) 와 샘플의 수가 충분하지 않을 때에는 앞에서 설명한 손으로(?) 또는 순진한(naive?) 방법으로는 한계가 부딪히게 된다.
여러 블랙 박스 모형들은 저마다의 가정(assumption)을 가지고 고차원성을 해결하려고 한다고 생각하면 될 것이다. 그 중에 최근에 각광을 받고 있는 방법은 심층 신경망이다. keras을 사용하여 심층 신경망 모형을 본격적으로 사용할 수 있다.
첫 번째 모형 : \(f(x) = \mathbb{E}[y|x]\)
데이터는 다음과 같다.
## Data
N = 1000
x = runif(N, -2, 2)
#y = sin(x*pi) + abs(x)/2 * rnorm(N)
y = sin(x*pi) + 1/2*(abs(cos(x))+1/8) * rnorm(N)
#y = sin(x*pi) + abs(cos(x)) * rnorm(N)
plot(y ~ x)
그리고, 실제 평균과 \(\pm\) 1 표준편차의 그래프는 다음과 같다.
curve(sin(x*pi), xlim=c(-2,2), ylim=c(-2,2))
curve(sin(x*pi) + 1/2*(abs(cos(x))+1/8), xlim=c(-2,2), add=TRUE, lty='dotted')
curve(sin(x*pi) - 1/2*(abs(cos(x))+1/8), xlim=c(-2,2), add=TRUE, lty='dotted')
이를 심층 인공신경망으로 추정해 보자.
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 = 100)
plot(history)
이제 모형의 예측과 실제 평균을 함께 그려보자. 회색은 실제 평균과 \(\pm 1\) 표준편차를 보여주고, 굵은 검정색 점선은 심층신경망으로 추정한 평균(회귀선)이다.
# PREDICTION
test_x <- seq(-2,2,0.05)
pred_y <- predict(model, test_x)[,1]
datPred <- data.frame(x=test_x, y=pred_y)
# TRUE
true_y <- sin(test_x*pi)
true_yUpper <- sin(test_x*pi) + 1/2*(abs(cos(test_x))+1/8)
true_yLower <- sin(test_x*pi) - 1/2*(abs(cos(test_x))+1/8)
datTrue <- data.frame(x=test_x, y=true_y, ymax=true_yUpper, ymin=true_yLower)
require(ggplot2)
ggplot() +
geom_line(data=datPred, aes(x=x, y=y), linetype='dotted') +
geom_line(data=datTrue, aes(x=x, y=y), col='grey20') +
geom_ribbon(data=datTrue, aes(x=x, ymin=ymin, ymax=ymax), alpha=0.15)
자료의 중심부에서는 상당히 정확하지만, 양끝에서는 그리 정확하지 못했다. 하지만 여기서 한 가지 더 주목할 점은 실제 데이터는 \(\mathbb{V}\textrm{ar}[y|x]\) 역시 변하고 있다는 점이다. 우리가 다룬 위의 심층 모형은 조건부 분산은 다루지 못한다.
두 번째 모형 : \(f_1(x) = \mathbb{E}[y|x], f_2(x) = \mathbb{V}\textrm{ar}[y|x]\)
심층신경망의 출력을 두 개 이상으로 만들 수 있다. 흔히 다중 출력(multiple output, multi-head)라고 얘기한다.
keras를 사용하면 쉽다. 위의 모형에서 마지막 output만 바꿔주면 된다.
model %>%
layer_dense(units = 32, input_shape = 1) %>%
layer_activation_leaky_relu() %>%
layer_dense(units=32) %>%
layer_activation_leaky_relu() %>%
layer_dense(units = 2, activation = 'linear')
평균-표준편차 출력 모형의 손실함수
문제는 loss이다.
가장 손쉬운 방법은 \(-\log \mathbb{P}(y|x)\) 로 설정하는 것이다. 이때 \(\mathbb{P}(y|x)\) 에 \(\log\) 를 취하는 것은 우리가 이 손실합의 최소화하는 모수(신경망 weights)를 구하기 때문이다. 다음 식에서 확인할 수 있듯이, 이는
\[\log \mathbb{P}(y_1|x_1) + \log \mathbb{P}(y_2|x_2) = \log \mathbb{P}(y_1|x_1)\mathbb{P}(y_2|x_2)\]
결국 모든 관찰 \((x_1, y_1), (x_2, y_2), \cdots\) 이 서로 독립일 때, 모든 관찰에 대해 확률을 최대로 하는 최대우도추정값을 구하는 것과 같다.
여기서 \(y\) 의 조건부 확률 분포 \(\mathbb{P}(y|x)\) 는 다음과 같이 가장 만만한 정규분포를 가정하자.
\[\mathbb{P}(y | x) \sim \mathcal{N}(\mu(x), \sigma^2(x))\]
그리고 신공신경망의 두 출력 \(y_1, y_2\) 를 \(\mu(x), \sigma(x)\) 로 사용할 것이다.
\[\mathbb{P}(y | x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\]
\[\log \mathbb{P}(y | x) = -\log{\sigma}-\frac{1}{2}\log{2\pi}-\frac{(x-\mu)^2}{2\sigma^2}\]
\[-\log \mathbb{P}(y | x) = \log{\sigma}+\frac{1}{2}\log{2\pi}+\frac{(x-\mu)^2}{2\sigma^2}\]
여기서 \(\frac{1}{2}\log{2\pi}\) 는 상수이므로 \(-\log \mathbb{P}(y | x)\) 의 합을 최소화하는 \(\mu(x)\) , \(\sigma(x)\) 를 찾는 것은 다음을 최소화하는 \(\mu(x)\) , \(\sigma(x)\) 를 찾는 문제와 동일하다.
\[ \log{\sigma} + \frac{(x-\mu)^2}{2\sigma^2}\]
손실함수 구현
위의 손실함수는 keras에서 다음과 같이 구현할 수 있다.
gauss_loss <- function(t, pred) {
#k_mean(k_square(pred[1,] - t[1,]) + pred[2,])
#k_mean(k_square((pred[,1] - t[,1])/pred[,2])/2 + pred[,2])
k_mean(k_square((pred[,1] - t[,1])/k_exp(pred[,2]))/2 + pred[,2] )
}
model %>% compile(
optimizer = "adam",
loss = function(y_true, y_pred)
# tilted_loss(quantile, y_true, y_pred),
gauss_loss(y_true, y_pred),
metrics = "mae"
)
조금만 설명해 보자면, loss=function(y_true, y_pred) {} 에서 y_true와 y_pred는 모두 길이 2의 벡터이다. y_pred는 인공신경망이 출력하는 두 값을 담고 있는 길이 2의 벡터가 되고, y_true는 동일하게 길이 2의 벡터로 상정했다. 하지만 우리가 가지고 있는 데이터는 \((x,y)\) 형태로 \(y\) 뿐이기 때문에, y_true의 두 번째 원소는 0으로 채워넣었다.
그리고 gauss_loss는 y_true와 y_pred를 받아, 위에서 소개한 손실함수를 구하게 된다. 아, 아니다. 인공신경망의 출력은 음수도 가능하기 때문에 \(\exp\) 를 적용한 값이 표준편차가 되도록 조정을 하였다.
train_Y <- matrix(c(y, rep(0, length(y))), length(y), 2)
history <-
model %>% fit(x,
train_Y,
batch_size = 100,
epochs = 1000)
plot(history)
이제 다시 한번 실제 모형과 추정된 모형을 비교해보자.
# PREDICTION
test_x <- seq(-2,2,0.05)
pred_y <- predict(model, test_x)[,1]
pred_sd <- predict(model, test_x)[,2]
datPred <- data.frame(x=test_x, y=pred_y,
y1=pred_y+exp(pred_sd),
y2=pred_y-exp(pred_sd))
# TRUE
true_y <- sin(test_x*pi)
true_yUpper <- sin(test_x*pi) + 1/2*(abs(cos(test_x))+1/8)
true_yLower <- sin(test_x*pi) - 1/2*(abs(cos(test_x))+1/8)
datTrue <- data.frame(x=test_x, y=true_y, ymax=true_yUpper, ymin=true_yLower)
require(ggplot2)
ggplot() +
geom_line(data=datPred, aes(x=x, y=y), linetype='dotted') +
geom_line(data=datPred, aes(x=x, y=y1), linetype='dotted') +
geom_line(data=datPred, aes(x=x, y=y2), linetype='dotted') +
geom_line(data=datTrue, aes(x=x, y=y), col='grey20') +
geom_ribbon(data=datTrue, aes(x=x, ymin=ymin, ymax=ymax), alpha=0.15)
어째 output이 하나였던 경우보다 더 정확한 듯 보인다. multi-head인 경우에 regularization 효과가 있다고 알려져 있다.
빅데이터 시뮬레이션
## Data
N = 10000
x = runif(N, -2, 2)
#y = sin(x*pi) + abs(x)/2 * rnorm(N)
y = sin(x*pi) + 1/2*(abs(cos(x))+1/8) * rnorm(N)
#y = sin(x*pi) + abs(cos(x)) * rnorm(N)
plot(y ~ x)
## To initialize the weights
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 = 2, activation = 'linear')
gauss_loss <- function(t, pred) {
#k_mean(k_square(pred[1,] - t[1,]) + pred[2,])
#k_mean(k_square((pred[,1] - t[,1])/pred[,2])/2 + pred[,2])
k_mean(k_square((pred[,1] - t[,1])/k_exp(pred[,2]))/2 + pred[,2] )
}
model %>% compile(
optimizer = "adam",
loss = function(y_true, y_pred)
# tilted_loss(quantile, y_true, y_pred),
gauss_loss(y_true, y_pred),
metrics = "mae"
)
train_Y <- matrix(c(y, rep(0, length(y))), length(y), 2)
history <-
model %>% fit(x,
train_Y,
batch_size = 100,
epochs = 100)
plot(history)
# PREDICTION
test_x <- seq(-2,2,0.05)
pred_y <- predict(model, test_x)[,1]
pred_sd <- predict(model, test_x)[,2]
datPred <- data.frame(x=test_x, y=pred_y,
y1=pred_y+exp(pred_sd),
y2=pred_y-exp(pred_sd))
# TRUE
true_y <- sin(test_x*pi)
true_yUpper <- sin(test_x*pi) + 1/2*(abs(cos(test_x))+1/8)
true_yLower <- sin(test_x*pi) - 1/2*(abs(cos(test_x))+1/8)
datTrue <- data.frame(x=test_x, y=true_y, ymax=true_yUpper, ymin=true_yLower)
require(ggplot2)
ggplot() +
geom_line(data=datPred, aes(x=x, y=y), linetype='dotted') +
geom_line(data=datPred, aes(x=x, y=y1), linetype='dotted') +
geom_line(data=datPred, aes(x=x, y=y2), linetype='dotted') +
geom_line(data=datTrue, aes(x=x, y=y), col='grey20') +
geom_ribbon(data=datTrue, aes(x=x, ymin=ymin, ymax=ymax), alpha=0.15)
거의 완벽하다!
생각할 문제
- \(y\) 의 조건부 분포 \(\mathbb{P}[y|x]\) 를 정규분포가 아니라, laplace 분포, t-분포 등으로 바꿔 보자. 어떤 잇점이 있는가?
- 설명변수가 1차원인 경우임에도 정확한 추정을 위해 10000개의 데이터포인트가 필요했다. 과연 고차원의 자료에서 이 방법을 유용하게 사용할 수 있을까?
- 만약 \(\mathbb{V}\textrm{ar}[y|x]\) 는 크게 변함이 없다는 사전 지식(혹은 사전 분포)가 있다면 어떻게 손실 함수를 바꿀 수 있는가?
개인적인 의문
- 솔직히 통계학을 안다면 누구나 이 방법을 생각해 낼 수 있을 것이다. 하지만 간단히 해 본 인터넷 검색에서는 자료를 찾을 수 없었다. (이 방법은 평균의 불확실성을 표현할 수 있는 베이지안 인공신경망과 다름을 유의하자.) 과연 이 모형은 어떤 이름으로 해야 검색할 수 있을까?
- 생각해보니 이 모형의 이름은 따로 없지만 y|x의 분포(parameteric)를 예측한다는 점에서 MDN(Mixture Density Network)와 동일하다고 볼 수 있을 것 같습니다.
- 손실함수만 조금 바꿔서 할 수 있는 quantile regression에 대해서는 다음에서 언급되고 있다.(사실 다음의 사이트를 기초로 위의 예시를 만들었다.)
Leave a comment