loglikelihood <- function(lambda.hat){
        S <- dpois(y, lambda.hat, log=TRUE)
        sum(S)
}
plot(x=lambda, loglikelihood, type='o', main = "Estimate of lambda",
     xlab="lambda", ylab = "log-likelihood")
abline(v=lambda[which.max(loglikelihood)], col="red")
text(1.8, -7.0, labels=lambda[which.max(loglikelihood)])