logist.llik <- function(y,p.hat) {
  return(sum(y*log(p.hat/(1-p.hat)))+sum(log(1-p.hat)))
}

logist.phat <- function(b0,b1,x) {
  return(exp(b0+b1*x)/(1+exp(b0+b1*x)))
}

logist.score <- function(y,p.hat,x) {
  return(c(sum(y-p.hat),sum(x*(y-p.hat))))
} 

logist.hess <- function(p.hat,x) {
  h00 <- -1*sum(p.hat*(1-p.hat))
  h01 <- -1*sum(x*p.hat*(1-p.hat))
  h11 <- -1*sum(x^2*p.hat*(1-p.hat))
  return(matrix(c(h00,h01,h01,h11),ncol=2))
}

x <- c(rep(-3,25),rep(-2,25),rep(-1,25),rep(0,25),rep(1,25),rep(2,25),rep(3,25))

y <- c(rep(0,49),1,rep(0,20),rep(1,5),rep(0,14),rep(1,11),rep(0,6),rep(1,19),
rep(0,3),rep(1,47))

table(x,y)

b0.0 <- 0
b1.0 <- 0

p0 <- logist.phat(b0.0,b1.0,x)

table(p0)

logist.llik(y,p0)

# Iteration 1

sc <- logist.score(y,p0,x)

hess <- logist.hess(p0,x)

b.1 <- c(b0.0,b1.0)-solve(hess,sc)

b.1

p1 <- logist.phat(b.1[1],b.1[2],x)

plot(x,p1)

logist.llik(y,p1)

# Iteration 2

sc <- logist.score(y,p1,x)

hess <- logist.hess(p1,x)

b.2 <- b.1-solve(hess,sc)

b.2

p2 <- logist.phat(b.2[1],b.2[2],x)

plot(x,p2)

logist.llik(y,p2)

# Iteration 3

sc <- logist.score(y,p2,x)

hess <- logist.hess(p2,x)

b.3 <- b.2-solve(hess,sc)

b.3

p3 <- logist.phat(b.3[1],b.3[2],x)

plot(x,p3)

logist.llik(y,p3)

# Iteration 4

sc <- logist.score(y,p3,x)

hess <- logist.hess(p3,x)

b.4 <- b.3-solve(hess,sc)

b.4

p4 <- logist.phat(b.4[1],b.4[2],x)

plot(x,p4)

logist.llik(y,p4)