x1<- c(19, 11, 15, 17, 10, 14, 17, 11, 13, 18, 
       16, 18, 10, 12, 17, 15, 15, 12, 24, 17)

loglike1<- function(par, data){
  
  lambda <- par
  n<- length(data)
  like1<- log(lambda) * sum(data) - 2*n*lambda - sum(log(factorial(data)))
         
  return(like1)
}


results1<- optim(par = c(1), data = x1, fn = loglike1, method = c("L-BFGS-B"),
                lower = 0.00000001, control = list(fnscale = -1), 
                hessian = T)

lam <- results1$par
lam
## [1] 7.525
cov.est1 <- -1* solve(results1$hessian)
cov.est1
##          [,1]
## [1,] 0.188125
results1$par - qnorm(1 - 0.05/2) * sqrt(cov.est1[1, 1])
## [1] 6.674897
results1$par + qnorm(1 - 0.05/2) * sqrt(cov.est1[1, 1])
## [1] 8.375103
score1 <- function(par, data){
  lambda <- par
  
  n <- length(data)
  lambda.score <- sum(data)/lambda - 2*n
  
  return(lambda.score)
  
}

score1(results1$par, data = x1)
## [1] -2.22354e-07
lambda_pts <- seq(6.67,8.37, by=0.05)
loglike_pts <- loglike1(lambda_pts, data = x1)
plot(lambda_pts, loglike_pts)