n <- 100

y <- rpois(n,0.8)

table(y)

y0 <- sum(y==0)
y1 <- sum(y==1)
y2 <- n-y0-y1

llik <- function(y0,y1,y2,theta){
  return(-1*theta*(y0+y1)+y1*log(theta)+y2*log(1-exp(-1*theta)*(1+theta)))
}

dl <- function(y0,y1,y2,theta){
  return(-1*(y0+y1)+y1/theta+y2*theta*exp(-1*theta)/(1-exp(-1*theta)*(1+theta)))
}

theta <- seq(0.2,1,0.01)
plot(theta,dl(y0,y1,y2,theta))
plot(theta,llik(y0,y1,y2,theta))

nr.adjustment <- function(y0,y1,y2,theta){
  numr <- dl(y0,y1,y2,theta)
  a1 <- exp(-1*theta)*((1-theta)*(1-exp(-1*theta)*(1+theta))-theta^2*exp(-1*theta))
  a2 <- (1-exp(-1*theta)*(1+theta))^2
  denom <- -1*y1/theta^2+y2*a1/a2
  return(-1*numr/denom)
}

theta0 <- (y1+2*y2)/n
print(theta0)
print(llik(y0,y1,y2,theta0))

theta1 <- theta0+nr.adjustment(y0,y1,y2,theta0)
print(theta1)
print(llik(y0,y1,y2,theta1))

theta2 <- theta1+nr.adjustment(y0,y1,y2,theta1)
print(theta2)
print(llik(y0,y1,y2,theta2))

theta3 <- theta2+nr.adjustment(y0,y1,y2,theta2)
print(theta3)
print(llik(y0,y1,y2,theta3))


th1 <- seq(0.71,0.72,0.001)
print(cbind(th1,llik(y0,y1,y2,th1)))




f.score.adjustment <- function(y0,y1,y2,theta){
  numr <- (1-exp(-1*theta)*(1+theta))*(1-theta+theta^2)+(theta^3)*exp(-1*theta)
  denom <- theta*(1-exp(-1*theta)*(1+theta))
  return(dl(y0,y1,y2,theta)/(n*exp(-1*theta)*numr/denom))
}

theta0 <- (y1+2*y2)/n
print(theta0)
print(llik(y0,y1,y2,theta0))

theta1 <- theta0+f.score.adjustment(y0,y1,y2,theta0)
print(theta1)
print(llik(y0,y1,y2,theta1))

theta2 <- theta1+f.score.adjustment(y0,y1,y2,theta1)
print(theta2)
print(llik(y0,y1,y2,theta2))

theta3 <- theta2+f.score.adjustment(y0,y1,y2,theta2)
print(theta3)
print(llik(y0,y1,y2,theta3))