gradient_descent
This is an old revision of the document!
Table of Contents
Gradient Descent
task: y hat = a + b*x 의 regression 에서 a 와 b를 구하기
- regressin line으로 이루어지는 에측치로 결과되는 오차의 제곱의 합은 (sum of square residuals, ssr 혹은 sse) 최소값이 되어야 한다.
- 방법 1. (버릴 방법)
- 우선 b를 안다고 하고 정해두면 (가령 -2.4와 같이 – 비현실적이지만)
- 예측할 수 있는 a값의 range를 두고
- 하나씩 대입하여 그 때의 sum of square residual 을 (ssr 혹은 sse) 구해서 기록한다.
- 이때 우리는 sse값이 ()^2의 형태를 띄기 때문에 알파벳 U 형태를 띄는 것을 안다.
- 따라서 이 sse 그래프의 최소값이 될 때의 a값이 무엇인가를 구해보면 된다.
- 아래 r script는 이런 방법으로 구해 본것인데 아래처럼 두가지 단점이 있다.
- 우선 정확한 b값을 알 수는 없다.
- a 값을 장님 문고리 잡듯이 -40에서 40 사이의 점수를 0.1씩 증가시켜 보는 것이 너무 비효율적이다.
- 더불어 이렇게 구한 a값은 1/10 단위의 숫자를 구하게 되지 정확한 a값을 구하는 것은 아니다 (예를 들면 1.4로 a값을 구했는데 정확한 값은 1.38976 일 수도 있다.
rs01
library(ggplot2)
library(ggpmisc)
library(tidyverse)
library(data.table)
library(MASS)
# settle down
rm(list=ls())
rnorm2 <- function(n,mean,sd) {
mean + sd * scale(rnorm(n))
}
ss <- function(x) {
sum((x-mean(x))^2)
}
sp <- function(x, y) {
sum((x-mean(x))*(y-mean(y)))
}
# 1. mean of three variables
# iq, allowance, gpa
# x, x2, y
means <- c(108, 42, 3.2)
# allowance, iq, gpa
sds <- c(11, 8, 0.4) # standard deviations
# correlation matrix. note that correlation
# between v1 and v2 is near none (0.05)
corr_matrix <- matrix(c(1, 0.05, 0.4,
0.05, 1, 0.2,
0.4, 0.2, 1),
nrow = 3) # Target correlation matrix
# 2. covariance matrix
# diagonal matrix of SDs
sd_diag <- diag(sds)
sd_diag
# Calculate covariance matrix
cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
# 3. population data
# set.seed(32)
set.seed(432)
# n_pop <- 100000
# pop <- mvrnorm(n = n_pop, mu = means, Sigma = cov_matrix)
# pop <- as.data.frame(pop)
# colnames(pop) <- c("x", "x2", "y")
# cor(pop)
# n_sam <- 100
# sam <- pop |> slice_sample(n = n_sam, replace = T)
# head(sam)
# cor(sam)
n_sam <- 100
sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
sam <- as.data.frame(sam)
colnames(sam) <- c("x", "x2", "y")
cor(sam)
y <- sam$y
x <- sam$x
x2 <- sam$x2
# check with regression
mo <- lm(y ~ x, data = sam)
summary(mo)
# graph
ggplot(data = sam, aes(x = x, y = y)) +
geom_point() +
stat_poly_line() +
stat_poly_eq(use_label(c("eq", "R2"))) +
theme_classic()
a <- summary(mo)$coefficient[1]
b <- summary(mo)$coefficient[2]
cat(a, b)
# We know what b and a are from the proof
bc <- sp(y, x)/ss(x)
ac <- mean(y) - bc * mean(x)
cat(ac, bc)
# but,
# finding a with b being fixed
# set.seed(191)
# Initialize random betas
# 우선 b를 고정하고 a만
# 변화시켜서 이해
b <- summary(mo)$coefficients[2] # 원래 b값
a <- 0
b.init <- b
a.init <- a
b.init
a.init
# Predict function:
predict <- function(x, a, b){
return (a + b * x)
}
# we use sum of square of error which oftentimes become big
sseloss<- function(predictions, y) {
residuals <- (y - predictions)
return(sum(residuals^2))
}
sses <- c() # for sum of square residuals
as <- c() # for as (intercepts)
summary(x)
for (i in seq(from = -40, to = 40, by = 0.01)) {
pred <- predict(x, i, b)
sse <- sseloss(pred, y)
sses <- append(sses, sse)
as <- append(as, i)
}
length(sses)
length(as)
min(sses)
max(sses)
min.pos.sses <- which(sses == min(sses))
min.pos.sses
print(as[min.pos.sses])
summary(mo)
plot(as, sses, type='l', lwd=2)
abline(v=as[min.pos.sses])
text(x = as[min.pos.sses], y = median(sses), col='red',
labels = paste(" sse = ", round(min(sses),4),
"\n is minimum value
when a =", as[min.pos.sses]),
pos=4)
ro01
> library(ggplot2)
> library(ggpmisc)
> library(tidyverse)
> library(data.table)
> library(MASS)
>
> # settle down
> rm(list=ls())
>
> rnorm2 <- function(n,mean,sd) {
+ mean + sd * scale(rnorm(n))
+ }
> ss <- function(x) {
+ sum((x-mean(x))^2)
+ }
> sp <- function(x, y) {
+ sum((x-mean(x))*(y-mean(y)))
+ }
>
> # 1. mean of three variables
> # iq, allowance, gpa
> # x, x2, y
> means <- c(108, 42, 3.2)
>
> # allowance, iq, gpa
> sds <- c(11, 8, 0.4) # standard deviations
> # correlation matrix. note that correlation
> # between v1 and v2 is near none (0.05)
> corr_matrix <- matrix(c(1, 0.05, 0.4,
+ 0.05, 1, 0.2,
+ 0.4, 0.2, 1),
+ nrow = 3) # Target correlation matrix
> # 2. covariance matrix
> # diagonal matrix of SDs
> sd_diag <- diag(sds)
> sd_diag
[,1] [,2] [,3]
[1,] 11 0 0.0
[2,] 0 8 0.0
[3,] 0 0 0.4
>
> # Calculate covariance matrix
> cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
>
> # 3. population data
> # set.seed(32)
> set.seed(432)
> # n_pop <- 100000
> # pop <- mvrnorm(n = n_pop, mu = means, Sigma = cov_matrix)
> # pop <- as.data.frame(pop)
> # colnames(pop) <- c("x", "x2", "y")
> # cor(pop)
>
> # n_sam <- 100
> # sam <- pop |> slice_sample(n = n_sam, replace = T)
> # head(sam)
> # cor(sam)
>
> n_sam <- 100
> sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
> sam <- as.data.frame(sam)
> colnames(sam) <- c("x", "x2", "y")
> cor(sam)
x x2 y
x 1.0000000 0.1880415 0.4143930
x2 0.1880415 1.0000000 0.3373768
y 0.4143930 0.3373768 1.0000000
>
>
> y <- sam$y
> x <- sam$x
> x2 <- sam$x2
>
> # check with regression
> mo <- lm(y ~ x, data = sam)
> summary(mo)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
>
> # graph
> ggplot(data = sam, aes(x = x, y = y)) +
+ geom_point() +
+ stat_poly_line() +
+ stat_poly_eq(use_label(c("eq", "R2"))) +
+ theme_classic()
>
> a <- summary(mo)$coefficient[1]
> b <- summary(mo)$coefficient[2]
> cat(a, b)
1.614428 0.01447618> # We know what b and a are from the proof
> bc <- sp(y, x)/ss(x)
> ac <- mean(y) - bc * mean(x)
> cat(ac, bc)
1.614428 0.01447618>
> # but,
> # finding a with b being fixed
> # set.seed(191)
> # Initialize random betas
> # 우선 b를 고정하고 a만
> # 변화시켜서 이해
> b <- summary(mo)$coefficients[2] # 원래 b값
> a <- 0
>
> b.init <- b
> a.init <- a
> b.init
[1] 0.01447618
> a.init
[1] 0
>
> # Predict function:
> predict <- function(x, a, b){
+ return (a + b * x)
+ }
>
> # we use sum of square of error which oftentimes become big
> sseloss<- function(predictions, y) {
+ residuals <- (y - predictions)
+ return(sum(residuals^2))
+ }
>
> sses <- c() # for sum of square residuals
> as <- c() # for as (intercepts)
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
87.12 98.44 106.27 106.80 114.62 133.33
> for (i in seq(from = -40, to = 40, by = 0.01)) {
+ pred <- predict(x, i, b)
+ sse <- sseloss(pred, y)
+ sses <- append(sses, sse)
+ as <- append(as, i)
+ }
> length(sses)
[1] 8001
> length(as)
[1] 8001
>
> min(sses)
[1] 10.97429
> max(sses)
[1] 173187
> min.pos.sses <- which(sses == min(sses))
> min.pos.sses
[1] 4162
> print(as[min.pos.sses])
[1] 1.61
> summary(mo)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
> plot(as, sses, type='l', lwd=2)
> abline(v=as[min.pos.sses])
> text(x = as[min.pos.sses], y = median(sses), col='red',
+ labels = paste(" sse = ", round(min(sses),4),
+ "\n is minimum value
+ when a =", as[min.pos.sses]),
+ pos=4)
>
위의 두번째 그래프는 a값이 -40 에서 40 까지 0.1 단위로 변할 때의 sse값을 구한 후에 이를 기록하여 그래프로 만든 것이다. 이렇게 구한 sse 값들 중 최소값일 때의 a값을 regression의 a값으로 추정한 것이다. 이렇게 해서 구한 값은 summary(mo)에서의 a와 값과 같다.
SSE 대신에 MSE를 쓰기
위에서 sse값을 저장해서 그 중에 최소값이 13.98253 임을 알게 된것이다. 그러나 이 계산은 각각의 sse값이 아주 크게 되는 경우가 많다.
> min(sses)
[1] 13.98253
> max(sses)
[1] 180081.3
따라서 sse 대신에 mse 값을 구해서 저장한다. mse = sse 를 n으로 나눈 값 (mean square error (residual)).
rs.01
#####
# with mean square error (mse) instead of sse
b <- summary(mo)$coefficients[2]
a <- 0
# we use sum of square of error which oftentimes become big
mseloss <- function(predictions, y) {
residuals <- (y - predictions)
return(mean(residuals^2))
}
mses <- c() # for sum of square residuals
as <- c() # for as (intercepts)
for (i in seq(from = -40, to = 40, by = 0.01)) {
pred <- predict(x, i, b)
mse <- mseloss(pred, y)
mses <- append(mses, mse)
as <- append(as, i)
}
length(mses)
length(as)
min(mses)
max(mses)
min.pos.mses <- which(mses == min(mses))
min.pos.mses
print(as[min.pos.mses])
summary(mo)
plot(as, mses, type='l', lwd=2)
abline(v=as[min.pos.mses])
text(x = as[min.pos.mses], y = median(mses), col='red',
labels = paste(" mse = ", round(min(mses),4),
"\n is minimum value
when a =", as[min.pos.mses]),
pos=4)
ro.02
> #####
> # with mean square error (mse) instead of sse
>
> b <- summary(mo)$coefficients[2]
> a <- 0
>
> # we use sum of square of error which oftentimes become big
> mseloss <- function(predictions, y) {
+ residuals <- (y - predictions)
+ return(mean(residuals^2))
+ }
>
> mses <- c() # for sum of square residuals
> as <- c() # for as (intercepts)
>
> for (i in seq(from = -40, to = 40, by = 0.01)) {
+ pred <- predict(x, i, b)
+ mse <- mseloss(pred, y)
+ mses <- append(mses, mse)
+ as <- append(as, i)
+ }
> length(mses)
[1] 8001
> length(as)
[1] 8001
>
> min(mses)
[1] 0.1097429
> max(mses)
[1] 1731.87
> min.pos.mses <- which(mses == min(mses))
> min.pos.mses
[1] 4162
> print(as[min.pos.mses])
[1] 1.61
> summary(mo)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
> plot(as, mses, type='l', lwd=2)
> abline(v=as[min.pos.mses])
> text(x = as[min.pos.mses], y = median(mses), col='red',
+ labels = paste(" mse = ", round(min(mses),4),
+ "\n is minimum value
+ when a =", as[min.pos.mses]),
+ pos=4)
>
>
b값 구하기
이제는 a값을 고정하고 b값도 같은 방식으로 구해볼 수 있다
rs01
##############################################
# b값도 범위를 추측한 후에 0.01씩 증가시키면서
# 각 b값에서 mse값을 구해본후 가장 작은 값을
# 가질 때의 b값을 구하면 된다.
# 그러나 b값의 적절한 구간을 예측하는 것이
# 불가능하다 (그냥 추측뿐)
# 위의 y 데이터에서 y = 314*x + rnorm(. . .)
# 이라면 -30-30 구간은 적절하지 않은 구간이 된다.
# 더우기 a값을 정확히 알아야 b값을 추출할 수 있다.
# 이는 적절한 방법이 아니다.
b <- 1
a <- summary(mo)$coefficients[1]
b.init <- b
a.init <- a
# Predict function:
predict <- function(x, a, b){
return (a + b * x)
}
# we use sum of square of error which oftentimes become big
mseloss <- function(predictions, y) {
residuals <- (y - predictions)
return(mean(residuals^2))
}
mses <- c()
bs <- c()
for (i in seq(from = -40, to = 40, by = 0.01)) {
pred <- predict(x, b, i)
mse <- mseloss(pred, y)
mses <- append(mses, mse)
bs <- append(bs,i)
}
min(mses)
max(mses)
min.pos.mses <- which(mses == min(mses))
print(bs[min.pos.mses])
summary(mo)
plot(bs, mses, type='l', lwd=2)
abline(v=bs[min.pos.mses])
text(x = bs[min.pos.mses], y = median(mses), col='red',
labels = paste(" mse = ", round(min(mses),4),
"\n is minimum value
when b =", round(bs[min.pos.mses],4)),
pos=4)
ro01
> ##############################################
> # b값도 범위를 추측한 후에 0.01씩 증가시키면서
> # 각 b값에서 mse값을 구해본후 가장 작은 값을
> # 가질 때의 b값을 구하면 된다.
> # 그러나 b값의 적절한 구간을 예측하는 것이
> # 불가능하다 (그냥 추측뿐)
> # 위의 y 데이터에서 y = 314*x + rnorm(. . .)
> # 이라면 -30-30 구간은 적절하지 않은 구간이 된다.
> # 더우기 a값을 정확히 알아야 b값을 추출할 수 있다.
> # 이는 적절한 방법이 아니다.
>
> b <- 1
> a <- summary(mo)$coefficients[1]
>
> b.init <- b
> a.init <- a
>
> # Predict function:
> predict <- function(x, a, b){
+ return (a + b * x)
+ }
>
> # we use sum of square of error which oftentimes become big
> mseloss <- function(predictions, y) {
+ residuals <- (y - predictions)
+ return(mean(residuals^2))
+ }
>
> mses <- c()
> bs <- c()
>
> for (i in seq(from = -40, to = 40, by = 0.01)) {
+ pred <- predict(x, b, i)
+ mse <- mseloss(pred, y)
+ mses <- append(mses, mse)
+ bs <- append(bs,i)
+ }
>
> min(mses)
[1] 0.1136351
> max(mses)
[1] 18442152
> min.pos.mses <- which(mses == min(mses))
> print(bs[min.pos.mses])
[1] 0.02
> summary(mo)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
> plot(bs, mses, type='l', lwd=2)
> abline(v=bs[min.pos.mses])
> text(x = bs[min.pos.mses], y = median(mses), col='red',
+ labels = paste(" mse = ", round(min(mses),4),
+ "\n is minimum value
+ when b =", round(bs[min.pos.mses],4)),
+ pos=4)
>
a와 b를 동시에 구할 수 있는 방법은 없을까? 위의 방법으로는 어렵다. 일반적으로 우리는 a와 b값이 무엇이되는가를 미분을 이용해서 구할 수 있었다. R에서 미분의 해를 구하기 보다는 해에 접근하도록 하는 프로그래밍을 써서 a와 b의 근사값을 구한다. 이것을 gradient descent라고 부른다.
Gradient descend
위에서 a값이 무엇일 때 mse값이 최소가 될까를 봐서 이 때의 a값을 실제 $y = a + bx$ 의 a 값으로 삼았다. 이 때 a값의 추정을 위해서
seq(-10, 10, 0.01) 의 범위와 증가값을 가지고 일일이 대입하여 mse를 구아였다.
위에서의 0.01씩 증가시켜 대입하는 것이 아니라 처음 한 숫자에서 시작한 후 그 다음 숫자를 정한 후에 점진적으로 그 숫자 간격을 줄여가면서 보면 정율적으로 0.01씩 증가시키는 것 보다 효율적일 것 같다. 이 증가분을 구하기 위해서 미분을 사용한다.
점차하강 = 조금씩 깍아서 원하는 기울기 (미분값) 찾기
prerequisite:
표준편차 추론에서 평균을 사용하는 이유: 실험적_수학적_이해
deriviation of a and b in a simple regression
위의 문서는 a, b에 대한 값을 미분법을 이용해서 직접 구하였다. 컴퓨터로는 이렇게 하기가 쉽지 않다. 그렇다면 이 값을 반복계산을 이용해서 추출하는 방법은 없을까? gradient descent
\begin{eqnarray*}
\text{for a (constant)} \\
\\
\text{SSE} & = & \text{Sum of Square Residuals} \\
\text{Residual} & = & (Y_i - (a + bX_i)) \\
\\
\frac{\text{dSSE}}{\text{da}}
& = & \frac{\text{dResidual^2}}{\text{dResidual}} * \frac{\text{dResidual}}{\text{da}} \\
& = & 2 * \text{Residual} * \dfrac{\text{d}}{\text{da}} (Y_i - (a+bX_i)) \\
& \because & \dfrac{\text{d}}{\text{da}} (Y_i - (a+bX_i)) = -1 \\
& = & 2 * \sum{(Y_i - (a + bX_i))} * -1 \\
& = & -2 *\sum{\text{residual}} \\
& .. & -2 \frac{\sum{\text{residual}}}{n} \\
& = & -2 * \overline{\text{residual}} \\
\end{eqnarray*}
위 미분결과는 주어진 Xi값에서 얻는 a값을 의미한다.
아래 R code에서 gradient function을 참조.
\begin{eqnarray*} \text{for b, (coefficient)} \\ \\ \dfrac{\text{d}}{\text{db}} \sum{(Y_i - (a + bX_i))^2} & = & \sum{\dfrac{\text{dResidual}^2}{\text{db}}} \\ & = & \sum{\dfrac{\text{dResidual}^2}{\text{dResidual}}*\dfrac{\text{dResidual}}{\text{db}} } \\ & = & \sum{2*\text{Residual} * \dfrac{\text{dResidual}}{\text{db}} } \\ & = & \sum{2*\text{Residual} * (-X_i) } \;\;\;\; \\ & \because & \dfrac{\text{dResidual}}{\text{db}} = (Y_i - (a+bX_i)) = -X_i \\ & = & -2 X_i \sum{(Y_i - (a + bX_i))} \\ & = & -2 * X_i * \sum{\text{residual}} \\ & .. & -2 * X_i * \frac{\sum{\text{residual}}}{n} \\ & = & -2 * \overline{X_i * \text{residual}} \\ \end{eqnarray*}
위의 미분결과는 주어진 Xi 값에서의 b값을 즉 기울기 값을 구하것을 의미한다.
위의 설명은 Sum of Square값을 미분하는 것을 전제로 하였지만, Mean Square 값을 (Sum of Square값을 N으로 나눈 것) 대용해서 이해할 수도 있다. 아래의 code는 (미분을 이해한다는 것을 전제로) b값과 a값이 변할 때 msr (mean square residual) 값이 어떻게 변하는가를 알려주는 것이다.
gradient <- function(x, y, predictions){
error = y - predictions
db = -2 * mean(x * error)
da = -2 * mean(error)
return(list("b" = db, "a" = da))
}
위 펑션으로 얻은 da와 db값을 초기에 설정한 a, b 값에 더해 준 값을 다시 a, b값으로 하여 gradient 펑션을 통해서 다시 db, da값을 구하고 이를 다시 이전 단계에서 구한 a, b값에 더하여 그 값을 다시 a, b값을 하여 . . . .
위를 반복한다. 단, db값과 da값을 그냥 대입하기 보다는 초기에 설정한 learning.rate값을 (0.01 예를 들면) 곱하여 구한 값을 더하게 된다. 이것이 아래의 code이다.
a 와 b 값을 gradient descent 방법을 이용하여 한꺼번에 구하기
rs01
# 아래는 a 와 b 를 랜덤하게 정한 후 (rnorm)
# 두 가지 일을 하게 된다. 첫 째는
# 이 때의 mse값을 저장한다. 둘째는
# a, b 값을 가질 때의 a 와 b에 대한 미분 값을
# 구한 후, 이 값에 1보다 작은 값을 곱하여 (learning rate)
# 구한 값을 다음의 a 와 b 값으로 써서 다시
# 이를 대입하여 mse값을 구하고 이 때의 미분 값을
# 구한다. 이 미분하여 구한 a, b값을 저장한 후에
# 다시 learning rate값을 곱하여 이를 다음의 a, b
# 값으로 사용한다. . . . 이를 반복한다.
a <- rnorm(1)
b <- rnorm(1)
a.init <- a
b.init <- b
# 이 때의 db, da는 무슨 값인가?
# 미분결과값이므로 db는 기울기값
# da는 그 당시의 a값
gradient <- function(x, y, predictions){
error = y - predictions
db = -2 * mean(x * error)
da = -2 * mean(error)
return(list("b" = db, "a" = da))
# return(c(db, da))
}
mseloss <- function(predictions, y) {
residuals <- (y - predictions)
return(mean(residuals^2))
}
# Train the model with scaled features
learning.rate = 0.1
# Record Loss for each epoch:
as = c()
bs = c()
das = c()
dbs = c()
ep <- c()
mses = c()
zx <- (x-mean(x))/sd(x)
nlen <- 75
for (epoch in 1:nlen) {
predictions <- predict(zx, a, b)
#residual <- residuals(predictions, y)
loss <- mseloss(predictions, y)
mses <- append(mses, loss)
grad <- gradient(zx, y, predictions)
db <- grad$b
da <- grad$a
dbs <- append(dbs, db)
das <- append(das, da)
ep <- append(ep, epoch)
step.b <- grad$b * learning.rate
step.a <- grad$a * learning.rate
b <- b-step.b
a <- a-step.a
as <- append(as, a)
bs <- append(bs, b)
}
tmp <- data.frame(ep, das, as, dbs, bs, mse)
tmp
# 위의 루프를 이해했는지 체크하기
a.init
b.init
tmp.p <- predict(zx, a.init, b.init)
# a.init b.init 일 때의 미분 결과
# 즉, 기울기와 절편 값
tmp.g <- gradient(zx, y, tmp.p)
tmp.g
# 이 값에 learning rate 값을 곱한 후에
a.step <- tmp.g$a*learning.rate
b.step <- tmp.g$b*learning.rate
# 이를 원래 a값에서 빼준 값을 다음 순서의
a.init <- a.init-a.step
b.init <- b.init-b.step
# a, b값으로 쓰고, 이 때의 미분 값을 구한 후에
# (루프 안에서), 미분 결과를 다시 learning rate
# 로 곱해 준 값을 (두번째 단계의) a, b 값에서
# 빼 준 값을 다음 단계의 a, b값으로 쓰고, 다시
# . . . . . . 위의 loop문은 이것을 75번 반복한 것
# scaled
# 이렇게 해서 구한 제일 마지막 a, b 값은 x 변인 대신에
# 표준화한 zx 변인을 사용해서 구한 값이므로 이 값을
# 다시 x 를 사용했을 때의 값으로 환원을 해야 한다 (unscale).
data.frame(a, b)
# unscale coefficients to make them comprehensible
# see http://commres.net/wiki/gradient_descent#why_normalize_scale_or_make_z-score_xi
# and
# http://commres.net/wiki/gradient_descent#how_to_unnormalize_unscale_a_and_b
a = a - (mean(x) / sd(x)) * b
b = b / sd(x)
cat(" unscaled a:", a, "\n unscaled b:", b, "\n")
# 아래 lm 결과와 확인
summary(mo)
# 아래는 a와 b를 구하는 과정에서 저장했던 temporary
# a, b 값들은 (as, bs 변인 내의) 모두 unscale한 것
as <- as - (mean(x) /sd(x)) * bs
bs <- bs / sd(x)
as
bs
# 이것을 그래프로 나타내기 위해서 parameters 라는
# 변인으로 저장
# 이하 if 문들은 그래프에 처음 a, b값과 루프 문에서의
# 그 값들이 변화하는 것을 모두 나타내기 위해서 사용한
# 것으로 결과와 상관없다.
parameters <- data.frame(as, bs, mses)
cat(paste0("Intercept: ", a, "\n", "Slope: ", b, "\n"))
summary(lm(y~x, data=sam))$coefficients
if (as[nlen] > as[1]) {
amax <- as[nlen]
amin <- as[1]
} else {
amax <- as[1]
amin <- as[nlen]
}
if (max(y) > amax) {
y.max <- max(y)
} else {
y.max <- amax
}
if (min(y) > amin) {
y.min <- amin
} else {
y.min <- min(y)
}
max(y)
min(y)
y.max
y.min
ggplot(sam, aes(x = x, y = y)) +
geom_point(size = 2) +
coord_cartesian(xlim = c(0, max(x)+30), ylim = c(y.min-2,y.max+2)) +
geom_abline(aes(intercept = as, slope = bs),
data = parameters, linewidth = 0.5,
color = 'green') +
geom_abline(aes(intercept = as, slope = bs),
data = parameters %>% slice_head(),
linewidth = 1, color = 'blue') +
geom_abline(aes(intercept = as, slope = bs),
data = parameters %>% slice_tail(),
linewidth = 1, color = 'red') +
# stat_poly_line() +
stat_poly_eq(use_label(c("eq", "R2"))) +
theme_classic() +
labs(title = 'Gradient descent. blue: start, red: end, green: gradients')
summary(lm(y~x))
data.frame(head(as), head(bs))
data.frame(tail(as), tail(bs))
a
b
ggplot(sam, aes(x = x, y = y)) +
geom_point(size = 2) +
geom_abline(aes(intercept = as, slope = bs),
data = parameters %>% slice_tail(),
linewidth = 1, color = 'red') +
stat_poly_line() +
stat_poly_eq(use_label(c("eq", "R2"))) +
theme_classic() +
labs(title = 'Gradient descent. blue: start, red: end, green: gradients')
ro01
>
> # 아래는 a 와 b 를 랜덤하게 정한 후 (rnorm)
> # 두 가지 일을 하게 된다. 첫 째는
> # 이 때의 mse값을 저장한다. 둘째는
> # a, b 값을 가질 때의 a 와 b에 대한 미분 값을
> # 구한 후, 이 값에 1보다 작은 값을 곱하여 (learning rate)
> # 구한 값을 다음의 a 와 b 값으로 써서 다시
> # 이를 대입하여 mse값을 구하고 이 때의 미분 값을
> # 구한다. 이 미분하여 구한 a, b값을 저장한 후에
> # 다시 learning rate값을 곱하여 이를 다음의 a, b
> # 값으로 사용한다. . . . 이를 반복한다.
> a <- rnorm(1)
> b <- rnorm(1)
> a.init <- a
> b.init <- b
>
> # 이 때의 db, da는 무슨 값인가?
> # 미분결과값이므로 db는 기울기값
> # da는 그 당시의 a값
> gradient <- function(x, y, predictions){
+ error = y - predictions
+ db = -2 * mean(x * error)
+ da = -2 * mean(error)
+ return(list("b" = db, "a" = da))
+ # return(c(db, da))
+ }
>
> mseloss <- function(predictions, y) {
+ residuals <- (y - predictions)
+ return(mean(residuals^2))
+ }
>
> # Train the model with scaled features
> learning.rate = 0.1
>
> # Record Loss for each epoch:
> as = c()
> bs = c()
> das = c()
> dbs = c()
> ep <- c()
>
> mses = c()
> zx <- (x-mean(x))/sd(x)
>
>
> nlen <- 75
> for (epoch in 1:nlen) {
+ predictions <- predict(zx, a, b)
+ #residual <- residuals(predictions, y)
+ loss <- mseloss(predictions, y)
+ mses <- append(mses, loss)
+
+ grad <- gradient(zx, y, predictions)
+ db <- grad$b
+ da <- grad$a
+
+ dbs <- append(dbs, db)
+ das <- append(das, da)
+ ep <- append(ep, epoch)
+
+ step.b <- grad$b * learning.rate
+ step.a <- grad$a * learning.rate
+ b <- b-step.b
+ a <- a-step.a
+
+ as <- append(as, a)
+ bs <- append(bs, b)
+ }
>
> mses
[1] 18.5823236 11.9374649 7.6828723 4.9587220 3.2144868 2.0976752 1.3825935 0.9247339 0.6315705
[10] 0.4438601 0.3236703 0.2467134 0.1974381 0.1658872 0.1456852 0.1327499 0.1244674 0.1191640
[19] 0.1157683 0.1135940 0.1122018 0.1113103 0.1107395 0.1103740 0.1101400 0.1099901 0.1098942
[28] 0.1098327 0.1097934 0.1097682 0.1097520 0.1097417 0.1097351 0.1097309 0.1097282 0.1097264
[37] 0.1097253 0.1097246 0.1097241 0.1097238 0.1097237 0.1097235 0.1097235 0.1097234 0.1097234
[46] 0.1097234 0.1097234 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233
[55] 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233
[64] 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233 0.1097233
[73] 0.1097233 0.1097233 0.1097233
> as
[1] -0.1210308 0.5352713 1.0603129 1.4803462 1.8163729 2.0851942 2.3002513 2.4722969 2.6099334
[10] 2.7200426 2.8081300 2.8785999 2.9349758 2.9800766 3.0161571 3.0450216 3.0681132 3.0865865
[19] 3.1013651 3.1131880 3.1226463 3.1302129 3.1362662 3.1411089 3.1449830 3.1480823 3.1505617
[28] 3.1525453 3.1541321 3.1554016 3.1564172 3.1572296 3.1578796 3.1583996 3.1588156 3.1591484
[37] 3.1594146 3.1596276 3.1597980 3.1599343 3.1600433 3.1601305 3.1602003 3.1602562 3.1603008
[46] 3.1603366 3.1603652 3.1603880 3.1604063 3.1604210 3.1604327 3.1604420 3.1604495 3.1604555
[55] 3.1604603 3.1604642 3.1604672 3.1604697 3.1604716 3.1604732 3.1604745 3.1604755 3.1604763
[64] 3.1604769 3.1604774 3.1604779 3.1604782 3.1604784 3.1604787 3.1604788 3.1604790 3.1604791
[73] 3.1604792 3.1604792 3.1604793
> bs
[1] -0.882887399 -0.678061874 -0.513791804 -0.382047207 -0.276388041 -0.191649390 -0.123688991
[8] -0.069184752 -0.025472351 0.009584993 0.037700984 0.060250009 0.078334326 0.092837949
[15] 0.104469854 0.113798642 0.121280331 0.127280644 0.132092896 0.135952322 0.139047582
[22] 0.141529980 0.143520863 0.145117551 0.146398096 0.147425092 0.148248743 0.148909311
[29] 0.149439087 0.149863967 0.150204721 0.150478005 0.150697180 0.150872957 0.151013931
[36] 0.151126992 0.151217667 0.151290388 0.151348711 0.151395485 0.151432999 0.151463084
[43] 0.151487213 0.151506564 0.151522084 0.151534530 0.151544513 0.151552518 0.151558939
[50] 0.151564088 0.151568218 0.151571530 0.151574187 0.151576317 0.151578026 0.151579396
[57] 0.151580495 0.151581376 0.151582083 0.151582650 0.151583104 0.151583469 0.151583762
[64] 0.151583996 0.151584184 0.151584335 0.151584456 0.151584553 0.151584631 0.151584693
[71] 0.151584743 0.151584783 0.151584816 0.151584841 0.151584862
>
> ep
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
[34] 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
[67] 67 68 69 70 71 72 73 74 75
> das
[1] -8.203776e+00 -6.563021e+00 -5.250416e+00 -4.200333e+00 -3.360267e+00 -2.688213e+00 -2.150571e+00
[8] -1.720456e+00 -1.376365e+00 -1.101092e+00 -8.808737e-01 -7.046990e-01 -5.637592e-01 -4.510073e-01
[15] -3.608059e-01 -2.886447e-01 -2.309158e-01 -1.847326e-01 -1.477861e-01 -1.182289e-01 -9.458309e-02
[22] -7.566648e-02 -6.053318e-02 -4.842654e-02 -3.874124e-02 -3.099299e-02 -2.479439e-02 -1.983551e-02
[29] -1.586841e-02 -1.269473e-02 -1.015578e-02 -8.124626e-03 -6.499701e-03 -5.199761e-03 -4.159809e-03
[36] -3.327847e-03 -2.662277e-03 -2.129822e-03 -1.703858e-03 -1.363086e-03 -1.090469e-03 -8.723751e-04
[43] -6.979001e-04 -5.583200e-04 -4.466560e-04 -3.573248e-04 -2.858599e-04 -2.286879e-04 -1.829503e-04
[50] -1.463603e-04 -1.170882e-04 -9.367056e-05 -7.493645e-05 -5.994916e-05 -4.795933e-05 -3.836746e-05
[57] -3.069397e-05 -2.455518e-05 -1.964414e-05 -1.571531e-05 -1.257225e-05 -1.005780e-05 -8.046240e-06
[64] -6.436992e-06 -5.149594e-06 -4.119675e-06 -3.295740e-06 -2.636592e-06 -2.109273e-06 -1.687419e-06
[71] -1.349935e-06 -1.079948e-06 -8.639584e-07 -6.911667e-07 -5.529334e-07
> dbs
[1] -2.553934e+00 -2.048255e+00 -1.642701e+00 -1.317446e+00 -1.056592e+00 -8.473865e-01 -6.796040e-01
[8] -5.450424e-01 -4.371240e-01 -3.505734e-01 -2.811599e-01 -2.254902e-01 -1.808432e-01 -1.450362e-01
[15] -1.163191e-01 -9.328788e-02 -7.481688e-02 -6.000314e-02 -4.812252e-02 -3.859426e-02 -3.095260e-02
[22] -2.482398e-02 -1.990883e-02 -1.596688e-02 -1.280544e-02 -1.026996e-02 -8.236511e-03 -6.605682e-03
[29] -5.297757e-03 -4.248801e-03 -3.407538e-03 -2.732846e-03 -2.191742e-03 -1.757777e-03 -1.409737e-03
[36] -1.130609e-03 -9.067487e-04 -7.272125e-04 -5.832244e-04 -4.677460e-04 -3.751323e-04 -3.008561e-04
[43] -2.412866e-04 -1.935118e-04 -1.551965e-04 -1.244676e-04 -9.982301e-05 -8.005805e-05 -6.420656e-05
[50] -5.149366e-05 -4.129791e-05 -3.312093e-05 -2.656298e-05 -2.130351e-05 -1.708542e-05 -1.370250e-05
[57] -1.098941e-05 -8.813506e-06 -7.068432e-06 -5.668882e-06 -4.546444e-06 -3.646248e-06 -2.924291e-06
[64] -2.345281e-06 -1.880915e-06 -1.508494e-06 -1.209812e-06 -9.702695e-07 -7.781561e-07 -6.240812e-07
[71] -5.005131e-07 -4.014115e-07 -3.219321e-07 -2.581895e-07 -2.070680e-07
> tmp <- data.frame(ep, das, as, dbs, bs, mse)
> tmp
ep das as dbs bs mse
1 1 -8.203776e+00 -0.1210308 -2.553934e+00 -0.882887399 18404982
2 2 -6.563021e+00 0.5352713 -2.048255e+00 -0.678061874 18404982
3 3 -5.250416e+00 1.0603129 -1.642701e+00 -0.513791804 18404982
4 4 -4.200333e+00 1.4803462 -1.317446e+00 -0.382047207 18404982
5 5 -3.360267e+00 1.8163729 -1.056592e+00 -0.276388041 18404982
6 6 -2.688213e+00 2.0851942 -8.473865e-01 -0.191649390 18404982
7 7 -2.150571e+00 2.3002513 -6.796040e-01 -0.123688991 18404982
8 8 -1.720456e+00 2.4722969 -5.450424e-01 -0.069184752 18404982
9 9 -1.376365e+00 2.6099334 -4.371240e-01 -0.025472351 18404982
10 10 -1.101092e+00 2.7200426 -3.505734e-01 0.009584993 18404982
11 11 -8.808737e-01 2.8081300 -2.811599e-01 0.037700984 18404982
12 12 -7.046990e-01 2.8785999 -2.254902e-01 0.060250009 18404982
13 13 -5.637592e-01 2.9349758 -1.808432e-01 0.078334326 18404982
14 14 -4.510073e-01 2.9800766 -1.450362e-01 0.092837949 18404982
15 15 -3.608059e-01 3.0161571 -1.163191e-01 0.104469854 18404982
16 16 -2.886447e-01 3.0450216 -9.328788e-02 0.113798642 18404982
17 17 -2.309158e-01 3.0681132 -7.481688e-02 0.121280331 18404982
18 18 -1.847326e-01 3.0865865 -6.000314e-02 0.127280644 18404982
19 19 -1.477861e-01 3.1013651 -4.812252e-02 0.132092896 18404982
20 20 -1.182289e-01 3.1131880 -3.859426e-02 0.135952322 18404982
21 21 -9.458309e-02 3.1226463 -3.095260e-02 0.139047582 18404982
22 22 -7.566648e-02 3.1302129 -2.482398e-02 0.141529980 18404982
23 23 -6.053318e-02 3.1362662 -1.990883e-02 0.143520863 18404982
24 24 -4.842654e-02 3.1411089 -1.596688e-02 0.145117551 18404982
25 25 -3.874124e-02 3.1449830 -1.280544e-02 0.146398096 18404982
26 26 -3.099299e-02 3.1480823 -1.026996e-02 0.147425092 18404982
27 27 -2.479439e-02 3.1505617 -8.236511e-03 0.148248743 18404982
28 28 -1.983551e-02 3.1525453 -6.605682e-03 0.148909311 18404982
29 29 -1.586841e-02 3.1541321 -5.297757e-03 0.149439087 18404982
30 30 -1.269473e-02 3.1554016 -4.248801e-03 0.149863967 18404982
31 31 -1.015578e-02 3.1564172 -3.407538e-03 0.150204721 18404982
32 32 -8.124626e-03 3.1572296 -2.732846e-03 0.150478005 18404982
33 33 -6.499701e-03 3.1578796 -2.191742e-03 0.150697180 18404982
34 34 -5.199761e-03 3.1583996 -1.757777e-03 0.150872957 18404982
35 35 -4.159809e-03 3.1588156 -1.409737e-03 0.151013931 18404982
36 36 -3.327847e-03 3.1591484 -1.130609e-03 0.151126992 18404982
37 37 -2.662277e-03 3.1594146 -9.067487e-04 0.151217667 18404982
38 38 -2.129822e-03 3.1596276 -7.272125e-04 0.151290388 18404982
39 39 -1.703858e-03 3.1597980 -5.832244e-04 0.151348711 18404982
40 40 -1.363086e-03 3.1599343 -4.677460e-04 0.151395485 18404982
41 41 -1.090469e-03 3.1600433 -3.751323e-04 0.151432999 18404982
42 42 -8.723751e-04 3.1601305 -3.008561e-04 0.151463084 18404982
43 43 -6.979001e-04 3.1602003 -2.412866e-04 0.151487213 18404982
44 44 -5.583200e-04 3.1602562 -1.935118e-04 0.151506564 18404982
45 45 -4.466560e-04 3.1603008 -1.551965e-04 0.151522084 18404982
46 46 -3.573248e-04 3.1603366 -1.244676e-04 0.151534530 18404982
47 47 -2.858599e-04 3.1603652 -9.982301e-05 0.151544513 18404982
48 48 -2.286879e-04 3.1603880 -8.005805e-05 0.151552518 18404982
49 49 -1.829503e-04 3.1604063 -6.420656e-05 0.151558939 18404982
50 50 -1.463603e-04 3.1604210 -5.149366e-05 0.151564088 18404982
51 51 -1.170882e-04 3.1604327 -4.129791e-05 0.151568218 18404982
52 52 -9.367056e-05 3.1604420 -3.312093e-05 0.151571530 18404982
53 53 -7.493645e-05 3.1604495 -2.656298e-05 0.151574187 18404982
54 54 -5.994916e-05 3.1604555 -2.130351e-05 0.151576317 18404982
55 55 -4.795933e-05 3.1604603 -1.708542e-05 0.151578026 18404982
56 56 -3.836746e-05 3.1604642 -1.370250e-05 0.151579396 18404982
57 57 -3.069397e-05 3.1604672 -1.098941e-05 0.151580495 18404982
58 58 -2.455518e-05 3.1604697 -8.813506e-06 0.151581376 18404982
59 59 -1.964414e-05 3.1604716 -7.068432e-06 0.151582083 18404982
60 60 -1.571531e-05 3.1604732 -5.668882e-06 0.151582650 18404982
61 61 -1.257225e-05 3.1604745 -4.546444e-06 0.151583104 18404982
62 62 -1.005780e-05 3.1604755 -3.646248e-06 0.151583469 18404982
63 63 -8.046240e-06 3.1604763 -2.924291e-06 0.151583762 18404982
64 64 -6.436992e-06 3.1604769 -2.345281e-06 0.151583996 18404982
65 65 -5.149594e-06 3.1604774 -1.880915e-06 0.151584184 18404982
66 66 -4.119675e-06 3.1604779 -1.508494e-06 0.151584335 18404982
67 67 -3.295740e-06 3.1604782 -1.209812e-06 0.151584456 18404982
68 68 -2.636592e-06 3.1604784 -9.702695e-07 0.151584553 18404982
69 69 -2.109273e-06 3.1604787 -7.781561e-07 0.151584631 18404982
70 70 -1.687419e-06 3.1604788 -6.240812e-07 0.151584693 18404982
71 71 -1.349935e-06 3.1604790 -5.005131e-07 0.151584743 18404982
72 72 -1.079948e-06 3.1604791 -4.014115e-07 0.151584783 18404982
73 73 -8.639584e-07 3.1604792 -3.219321e-07 0.151584816 18404982
74 74 -6.911667e-07 3.1604792 -2.581895e-07 0.151584841 18404982
75 75 -5.529334e-07 3.1604793 -2.070680e-07 0.151584862 18404982
>
> a.init
[1] -0.9414084
> b.init
[1] -1.138281
> tmp.p <- predict(zx, a.init, b.init)
> # a.init b.init 일 때의 미분 결과
> # 즉, 기울기와 절편 값
> tmp.g <- gradient(zx, y, tmp.p)
> tmp.g
$b
[1] -2.553934
$a
[1] -8.203776
>
> a.step <- tmp.g$a*learning.rate
> b.step <- tmp.g$b*learning.rate
>
> a.init <- a.init-a.step
> b.init <- b.init-b.step
>
> # scaled
> data.frame(a, b)
a b
1 3.160479 0.1515849
>
> # unscale coefficients to make them comprehensible
> # see http://commres.net/wiki/gradient_descent#why_normalize_scale_or_make_z-score_xi
> # and
> # http://commres.net/wiki/gradient_descent#how_to_unnormalize_unscale_a_and_b
> #
> a = a - (mean(x) / sd(x)) * b
> b = b / sd(x)
> cat(" unscaled a:", a, "\n unscaled b:", b, "\n")
unscaled a: 1.614429
unscaled b: 0.01447617
> summary(mo)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
>
> # changes of estimators
> as <- as - (mean(x) /sd(x)) * bs
> bs <- bs / sd(x)
>
> as
[1] 8.883749 7.450986 6.300598 5.376937 4.635321 4.039872 3.561785 3.177929 2.869732 2.622283 2.423609
[12] 2.264096 2.136025 2.033200 1.950644 1.884362 1.831146 1.788421 1.754118 1.726578 1.704467 1.686715
[23] 1.672463 1.661020 1.651834 1.644459 1.638537 1.633784 1.629967 1.626903 1.624443 1.622468 1.620883
[34] 1.619610 1.618588 1.617768 1.617109 1.616581 1.616156 1.615816 1.615542 1.615322 1.615146 1.615004
[45] 1.614891 1.614800 1.614726 1.614668 1.614620 1.614583 1.614552 1.614528 1.614508 1.614492 1.614480
[56] 1.614470 1.614462 1.614455 1.614450 1.614446 1.614442 1.614439 1.614437 1.614436 1.614434 1.614433
[67] 1.614432 1.614431 1.614431 1.614430 1.614430 1.614430 1.614429 1.614429 1.614429
> bs
[1] -0.0843146613 -0.0647540755 -0.0490664857 -0.0364850387 -0.0263947182 -0.0183022811 -0.0118121466
[8] -0.0066070587 -0.0024325783 0.0009153551 0.0036003976 0.0057538017 0.0074808319 0.0088659100
[15] 0.0099767427 0.0108676305 0.0115821225 0.0121551451 0.0126147092 0.0129832796 0.0132788731
[22] 0.0135159391 0.0137060660 0.0138585478 0.0139808382 0.0140789151 0.0141575727 0.0142206562
[29] 0.0142712491 0.0143118246 0.0143443662 0.0143704645 0.0143913954 0.0144081820 0.0144216448
[36] 0.0144324420 0.0144411013 0.0144480461 0.0144536158 0.0144580827 0.0144616652 0.0144645383
[43] 0.0144668426 0.0144686906 0.0144701727 0.0144713614 0.0144723147 0.0144730792 0.0144736924
[50] 0.0144741841 0.0144745785 0.0144748948 0.0144751485 0.0144753519 0.0144755151 0.0144756460
[57] 0.0144757509 0.0144758351 0.0144759026 0.0144759567 0.0144760001 0.0144760350 0.0144760629
[64] 0.0144760853 0.0144761032 0.0144761177 0.0144761292 0.0144761385 0.0144761459 0.0144761519
[71] 0.0144761566 0.0144761605 0.0144761636 0.0144761660 0.0144761680
>
> parameters <- data.frame(as, bs, mses)
>
> cat(paste0("Intercept: ", a, "\n", "Slope: ", b, "\n"))
Intercept: 1.61442902414182
Slope: 0.0144761679939984
> summary(lm(y~x, data=sam))$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.61442839 0.344622352 4.684631 9.041202e-06
x 0.01447618 0.003211564 4.507515 1.817763e-05
>
> if (as[nlen] > as[1]) {
+ amax <- as[nlen]
+ amin <- as[1]
+ } else {
+ amax <- as[1]
+ amin <- as[nlen]
+ }
> if (max(y) > amax) {
+ y.max <- max(y)
+ } else {
+ y.max <- amax
+ }
> if (min(y) > amin) {
+ y.min <- amin
+ } else {
+ y.min <- min(y)
+ }
> max(y)
[1] 3.878048
> min(y)
[1] 2.092215
> y.max
[1] 8.883749
> y.min
[1] 1.614429
>
> ggplot(sam, aes(x = x, y = y)) +
+ geom_point(size = 2) +
+ coord_cartesian(xlim = c(0, max(x)+30), ylim = c(y.min-2,y.max+2)) +
+ geom_abline(aes(intercept = as, slope = bs),
+ data = parameters, linewidth = 0.5,
+ color = 'green') +
+ geom_abline(aes(intercept = as, slope = bs),
+ data = parameters %>% slice_head(),
+ linewidth = 1, color = 'blue') +
+ geom_abline(aes(intercept = as, slope = bs),
+ data = parameters %>% slice_tail(),
+ linewidth = 1, color = 'red') +
+ # stat_poly_line() +
+ stat_poly_eq(use_label(c("eq", "R2"))) +
+ theme_classic() +
+ labs(title = 'Gradient descent. blue: start, red: end, green: gradients')
>
> summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
>
> data.frame(head(as), head(bs))
head.as. head.bs.
1 8.883749 -0.08431466
2 7.450986 -0.06475408
3 6.300598 -0.04906649
4 5.376937 -0.03648504
5 4.635321 -0.02639472
6 4.039872 -0.01830228
> data.frame(tail(as), tail(bs))
tail.as. tail.bs.
1 1.614430 0.01447615
2 1.614430 0.01447616
3 1.614430 0.01447616
4 1.614429 0.01447616
5 1.614429 0.01447617
6 1.614429 0.01447617
> a
[1] 1.614429
> b
[1] 0.01447617
>
> ggplot(sam, aes(x = x, y = y)) +
+ geom_point(size = 2) +
+ geom_abline(aes(intercept = as, slope = bs),
+ data = parameters %>% slice_tail(),
+ linewidth = 1, color = 'red') +
+ stat_poly_line() +
+ stat_poly_eq(use_label(c("eq", "R2"))) +
+ theme_classic() +
+ labs(title = 'Gradient descent. blue: start, red: end, green: gradients')
>
Why normalize (scale or make z-score) xi
x 변인의 측정단위로 인해서 b 값이 결정되게 되는데 이 때의 b값은 상당하고 다양한 범위를 가질 수 있다. 가령 월 수입이 (인컴) X 라고 한다면 우리가 추정해야 (추적해야) 할 b값은 수백만이 될 수도 있다.이 값을 gradient로 추적하게 된다면 너무도 많은 iteration을 거쳐야 할 수 있다. 변인이 바뀌면 이 b의 추적범위도 드라마틱하게 바뀌게 된다. 이를 표준화한 x 점수를 사용하게 된다면 일정한 learning rate와 iteration만으로도 정확한 a와 b를 추적할 수 있게 된다.
How to unnormalize (unscale) a and b
\begin{eqnarray*} y & = & a + b * x \\ & & \text{we use z instead of x} \\ & & \text{and } \\ & & z = \frac{(x - \mu)}{\sigma} \\ & & \text{suppose that the result after calculation be } \\ y & = & k + m * z \\ & = & k + m * \frac{(x - \mu)}{\sigma} \\ & = & k + \frac{m * x}{\sigma} - \frac{m * \mu}{\sigma} \\ & = & k - \frac{m * \mu}{\sigma} + \frac{m * x}{\sigma} \\ & = & \underbrace{k - \frac{\mu}{\sigma} * m}_\text{ 1 } + \underbrace{\frac{m}{\sigma}}_\text{ 2 } * x \\ & & \text{therefore, a and be that we try to get are } \\ a & = & k - \frac{\mu}{\sigma} * m \\ b & = & \frac{m}{\sigma} \\ \end{eqnarray*}
gradient_descent.1778623341.txt.gz · Last modified: by hkimscil