logistic reg.

df = read.table(textConnection('
ID      admit   gre     gpa
1       0       380     3.61
2       1       660     3.67
3       1       800     4
4       1       640     3.19
5       0       520     2.93
6       1       760     3
7       1       560     2.98
8       0       400     3.08
9       1       540     3.39
10      0       700     3.92
11      0       800     4
12      0       440     3.22
13      1       760     4
14      0       700     3.08
15      1       700     4
16      0       480     3.44
17      0       780     3.87
18      0       360     2.56
19      0       800     3.75
20      1       540     3.81
21      0       500     3.17
22      1       660     3.63
23      0       600     2.82
24      0       680     3.19
25      1       760     3.35
26      1       800     3.66
27      1       620     3.61
28      1       520     3.74
29      1       780     3.22
30      0       520     3.29
31      0       540     3.78
32      0       760     3.35
33      0       600     3.4
34      1       800     4
35      0       360     3.14
36      0       400     3.05
37      0       580     3.25
38      0       520     2.9
39      1       500     3.13
40      1       520     2.68
41      0       560     2.42
42      1       580     3.32
43      1       600     3.15
44      0       500     3.31
45      0       700     2.94
46      1       460     3.45
47      1       580     3.46
48      0       500     2.97
49      0       440     2.48
50      0       400     3.35'), header = TRUE)
X = as.matrix(cbind(1, df[,3:4]))
y = df[,2]
f3<-function(beta){
  xb = X %*% beta
  pi_est = exp(xb)/(1+exp(xb))
  ypi = y-pi_est
  gradient = t(X) %*% ypi
  hessian = -t(X) %*% diag(as.vector(pi_est * (1-pi_est))) %*% X
  return(list(gradient, hessian))
}

newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
  x <- x0
  f3.x <- f3(x)
  n <- 0
  while ((max(abs(f3.x[[1]])) > tol) & (n < n.max)) {
    x <- x - solve(f3.x[[2]]) %*% f3.x[[1]]
    f3.x <- f3(x)
    n <- n + 1
  }
  if (n == n.max) {
    cat('newton failed to converge\n')
  } else {
    return(x)
  }
}

b1 = newton(f3,c(0, 0, 0),1e-6)
b1
coef(glm(admit ~ gre + gpa, df, family = 'binomial'))
-solve(f3(b1)[[2]])
vcov(glm(admit ~ gre + gpa, df, family = 'binomial'))