garybaylor / algorithms Goto Github PK

######################################################
## gradient descent algorithm for logistic regression
######################################################

## ==========
## FUNCTIONS
## ==========

## 1. function to perform one GD update
# INPUT: 
#   par -->> a vector of parameters, in the order of intercept, slope_1, slope_2, ..., slope_p
#   X   -->> the design matrix of dimension n * p 
#   y   -->> the response vector of length n
#   alpha->> learning rate
# OUTPUT:
#   a vector of updated parameters and the value of the cost
one_update <- function(par, X, y, alpha = 0.1) {
	nr <- nrow(X)
	X <- cbind(1, X)
	a <- sigma(X %*% par)
	dz <- a - y
	dpar <-(t(X) %*% dz) / nr
	J <- mean(L(a, y))
	
	list(par = par - alpha * dpar, cost = J)
}

## 2. function to estimate parameters for LR using GD algorithm
# INPUT:
#   X   -->> the design matrix of dimension n * p 
#   y   -->> the response vector of length n
#   alpha->> learning rate
#   max_iter  -->> maximum number of iterations that can be performed
#   tol       -->> the tolerance
# OUTPUT:
#   a matrix with each row representing the estimated parameters and the
#   cost at each iteration
GD <- function(X, y, max_iter = 1000, tol = 1e-8, alpha = 0.1) {
	nr <- nrow(X)
	nc <- ncol(X)
	par <- runif(nc + 1)
	J <- 0
    res <- matrix(NA, ncol = nc + 2, nrow = max_iter)
    colnames(res) <- c("b", paste0("w", 1:nc), "J")
    
	for(i in 1:max_iter) {
		tmp <- one_update(par, X, y, alpha)
		par_new <- tmp[[1]]
		J_new <- tmp[[2]]
	    res[i, ] <- c(par_new, J_new)
		if(abs(J_new - J) < tol) break
		par <- par_new
		J <- J_new
	}
	res <- res[1:i, ]
	class(res) <- "GD"
	res
}

## the print function for the LR model
print.GD <- function(x) {
	nr <- nrow(x)
	nc <- ncol(x)
    print(x[nr, -nc])
}

sigma <- function(x) {
	1 / (1 + exp(-x))
}

# the LOSS function
L <- function(a, y) {
	-(y * log(a) + (1 - y) * log(1 - a))
}

## ==========
## EXAMPLES 1
## ==========

set.seed(101)
X <- rnorm(200)
X <- matrix(X, ncol = 2)
set.seed(202)
y <- rbinom(100, 1, 0.3)

## glm
mod1 <- glm(y ~ X, family = "binomial")

## gradient descent algorithm
mod2 <- GD(X, y, alpha = 0.1)



## ==========
## EXAMPLES 1
## ==========

set.seed(101)
X <- rnorm(200)
X <- matrix(X, ncol = 2)
set.seed(202)
y <- rbinom(100, 1, 0.3)

## glm
mod1 <- glm(y ~ X, family = "binomial")

## gradient descent algorithm
mod2 <- GD(X, y, alpha = 0.1)
plot(mod2[, 4], type = "l", main = "Value of the Cost Function at each iteration",
     xlab = "Iteration", ylab = "Cost")
> mod1

Call:  glm(formula = y ~ X, family = "binomial")

Coefficients:
(Intercept)           X1           X2  
    -1.1036      -0.1134       0.1334  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    112.5 
Residual Deviance: 112 	AIC: 118
> mod2
         b         w1         w2 
-1.1025514 -0.1117739  0.1330478 

## ==========
## EXAMPLES 2
## ==========

set.seed(303)
X <- rnorm(1000)
X <- matrix(X, ncol = 10)
set.seed(404)
y <- rbinom(100, 1, 0.3)

## glm
mod1 <- glm(y ~ X, family = "binomial")

## gradient descent algorithm
mod2 <- GD(X, y, alpha = 0.1)

plot(mod2[, 12], type = "l", main = "Value of the Cost Function at each iteration",
     xlab = "Iteration", ylab = "Cost")

> mod1

Call:  glm(formula = y ~ X, family = "binomial")

Coefficients:
(Intercept)           X1           X2           X3           X4  
  -0.672672    -0.004971     0.175255    -0.164195    -0.088264  
         X5           X6           X7           X8           X9  
  -0.084391    -0.208116     0.167310     0.598571    -0.166198  
        X10  
   0.354586  

Degrees of Freedom: 99 Total (i.e. Null);  89 Residual
Null Deviance:	    126.8 
Residual Deviance: 117.2 	AIC: 139.2
> mod2
           b           w1           w2           w3           w4 
-0.672770788 -0.004100684  0.175682986 -0.164154564 -0.087450584 
          w5           w6           w7           w8           w9 
-0.082696729 -0.207435362  0.165922662  0.598113358 -0.165788577 
         w10 
 0.353785814