################################################################################
# logistic regression example
mydata = read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
model = glm(admit ~ ., data = mydata, family = binomial)
model$coefficients

X = as.matrix(cbind(rep(1, nrow(mydata)), mydata[,2:ncol(mydata)]))
y = mydata$admit

# logistic regression using Iteratively Reweighted Least Squares (IRLS)
beta = as.vector(array(0, ncol(X)))
for (i in 1:25) {
    predictions = 1 / (1 + exp(- (X %*% beta)))
    gradient = t(X) %*% (predictions - y)
    Hessian = t(X) %*% diag(as.vector(predictions * (1 - predictions))) %*% X
    beta = beta - solve(Hessian, diag(ncol(X))) %*% gradient
}
beta

summary(model)
predictions = predict(model, mydata, type="response")
C = solve(t(X) %*% diag(predictions * (1 - predictions)) %*% X, diag(ncol(X)))
sqrt(diag(C))

predictions = predict(model, mydata, type="response", se.fit = T)
stderr = predictions$fit * (1 - predictions$fit) * sqrt(diag(X %*% C %*% t(X)))


################################################################################
# gradient descent for logistic regression
n.pos = 30
n.neg = 30
n = n.pos + n.neg
mu.pos = c(3, 3)
mu.neg = c(-3, -3)
library(MASS)
X.pos = mvrnorm(n.pos, mu = mu.pos, Sigma = diag(length(mu.pos)))
X.neg = mvrnorm(n.neg, mu = mu.neg, Sigma = diag(length(mu.neg)))
X = cbind(rep(1, n), rbind(X.pos, X.neg))
y = c(rep(1, n.pos), rep(0, n.neg))
plot(X[,2:3], pch=y+1)
learningRate = 0.01
beta = c(0, 0, 0)

predictions = 1 / (1 + exp(-(X %*% beta)))
-mean(log(y * predictions + (1 - y) * (1 - predictions)))
exp(-0.6931472)    # exp(-log(0.5)) == 0.5

negative.gradient = c(0, 0, 0)
for (i in 1:n) {
    negative.gradient[1] = negative.gradient[1] + ((y[i] - predictions[i]) * X[i,1]) / n
    negative.gradient[2] = negative.gradient[2] + ((y[i] - predictions[i]) * X[i,2]) / n
    negative.gradient[3] = negative.gradient[3] + ((y[i] - predictions[i]) * X[i,3]) / n
}
beta[1] = beta[1] + learningRate * negative.gradient[1]
beta[2] = beta[2] + learningRate * negative.gradient[2]
beta[3] = beta[3] + learningRate * negative.gradient[3]

predictions = 1 / (1 + exp(-(X %*% beta)))
-mean(log(y * predictions + (1 - y) * (1 - predictions)))    # error is reduced



################################################################################
# logistic regression example for cancer diagnosis
start.time = Sys.time()

trn_X = read.csv("C:/Data/cancer/trn_X.csv", header = F)
trn_y = factor(scan("C:/Data/cancer/trn_y.txt"), levels = c(0, 1), labels = c("no", "yes"))
tst_X = read.csv("C:/Data/cancer/tst_X.csv", header = F)

trn = data.frame(y = trn_y, X = trn_X)
tst = data.frame(X = tst_X)
model = glm(y ~ ., data = trn, family = binomial)
predictions = predict(model, newdata = tst, type = "response")

output = data.frame(Index = 1:length(predictions), Prediction = predictions)
write.csv(output, "C:/Data/cancer/predictions.csv", quote=F, row.names = F)
Sys.time() - start.time