2-grad.R

# Sun Nov 14 20:13:12 2021 ------------------------------

# Logistic regression mtcars
library(RColorBrewer)

data("mtcars")

# Determine number of iterations
niter <- 100
# Determine learning rate / step size
alpha <- 1e-5
set.seed(1)
b0 <- rnorm(1) # intercept
b1 <- rnorm(1) # slope1
b2 <- rnorm(1) # slope2

# Set palette
cols <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(niter)

# Plot
layout(matrix(c(1,1,2,2,3), nrow = 1, ncol = 5))
plot(mpg ~ hp, data = mtcars, pch = 16,
     xlab = "Horsepower", ylab = "Miles per gallon (mpg)",
     col = ifelse(mtcars$vs, rgb(0,0,1,.5), rgb(1,0,0,.5)))
legend("topright", legend = c("Straight", "V-shaped"),
       col = c(rgb(0,0,1,.5), rgb(1,0,0,.5)), pch = 16,
       bty = "n")
# Perform gradient descent
slopes1 <- rep(NA, niter)
slopes2 <- rep(NA, niter)
intercepts <- rep(NA, niter)
for(i in 1:niter){
   # prediction
   logit <- b0 + b1 * mtcars$mpg + b2 * mtcars$hp
   y <- 1 / (1 + exp(-logit))
   # b0 = b0 - dJ/da * alpha
   b0 <- b0 - sum(y - mtcars$vs) * alpha
   # b1 = b1 - dJ/db * alpha
   b1 <- b1 - sum((y - mtcars$vs) * mtcars$mpg) * alpha
   # b2 = b2 - dJ/db * alpha      
   b2 <- b2 - sum((y - mtcars$vs) * mtcars$hp) * alpha
   # Solve for b0 + b1 * mtcars$mpg + b2 * mtcars$hp = 0
   # mtcars$mpg = -(b2 * mtcars$hp + b0) / b1
   db_points <- -(b2 * c(40, 400) + b0) / b1
   # Plot decision boundary
   lines(c(40, 400), db_points, col = cols[i], lty = 2)
   # Save estimates over all iterations
   intercepts[i] <- b0
   slopes1[i] <- b1
   slopes2[i] <- b2
}
title("Decision Boundary")

plot(x=1:niter, intercepts, type = "l", lty=1,
     ylim = c(-1, 1), xlab = "Iteration No.", ylab = "Coefficient value")
abline(h = 0, lty = 2, col="grey")
lines(x=1:niter, slopes1, lty=3)
lines(x=1:niter, slopes2, lty=4)
legend("topright",
       legend = c(expression(beta[0]),
                  expression(beta[1]),
                  expression(beta[2])),
       lty = c(1, 3, 4),
       bty = "n")
title("Coefficient Estimation")

# Add colorbar
z = matrix(1:niter, nrow = 1)
image(1, 1:niter, z,
      col = cols, axes = FALSE, xlab = "", ylab = "")
title("Iteration No.")
axis(2, at = c(1, seq(20, niter, by=20)))