y1 <- rpois(10, 5)
y2 <- rep(y1, 3)

llik.pois <- function(lambda, y) log(lambda)*sum(y) - length(y)*lambda 

lambda.seq <- seq(0,10, by = 0.1)

llik1 <- sapply(lambda.seq, function(x) llik.pois(x, y1))
llik2 <- sapply(lambda.seq, function(x) llik.pois(x, y2))

plot(lambda.seq, llik1 - llik.pois(mean(y1),y1), type = "l")
lines(lambda.seq, llik2 - llik.pois(mean(y2),y2), col = 2)

sum(y1)/mean(y1) - length(y1)
sum(y1)/mean(y1)
sum(y2)/mean(y2)

Nsim <- 10000
score <- J <- lambda.ML <- c()
llik.p <- matrix(NA, length(lambda.seq), Nsim)
for(i in 1:Nsim){
  y <- rpois(10, 5)
  score[i] <- sum(y)/5 - 30
  J[i] <- sum(y)/5^2
  lambda.ML[i] <- mean(y)
  llik <- sapply(lambda.seq, function(x) llik.pois(x, y))
  llik.p[, i] <- llik - llik.pois(mean(y), y)
}

matplot(lambda.seq, llik.p[,1:5], type = "l")
abline(v = 5)

mean(score)
var(score)
mean(J)
mean(score^2)

30/5

hist(lambda.ML, freq = F) 
curve(dnorm(x, 5, sqrt(5/10)), from = 3.5, to = 6.5, 
      add = T, col = 2)