y <- c( 0, rep( 1, 5 ), rep( 2, 4 ), rep( 3, 2 ), 4, 6 )

print( lambda.hat.mle <- mean( y ) )

#[1] 2.071429

lambda.grid <- seq( 0.001, 4.2, length = 500 )

poisson.log.likelihood <- function( lambda, s, n ) {

  ll <- s * log( lambda ) - n * lambda 

  return( ll )

}

poisson.likelihood <- function( lambda, s, n ) {

  return( exp( poisson.log.likelihood( lambda, s, n ) ) )

}

print( s <- sum( y ) )

[1] 29

print( n <- length( y ) )

[1] 14

par( mfrow = c( 1, 2 ) )

plot( lambda.grid, poisson.likelihood( lambda.grid, s, n ),
  type = 'l', lwd = 2, xlab = 'lambda',
  ylab = 'likelihood' )

abline( v = lambda.hat.mle, lty = 2, lwd = 2, col = 'red' )

plot( lambda.grid, poisson.log.likelihood( lambda.grid, s, n ),
  type = 'l', lwd = 2, xlab = 'lambda',
  ylab = 'log likelihood' )

abline( v = lambda.hat.mle, lty = 2, lwd = 2, col = 'red' )

print( se.hat <- sqrt( lambda.hat.mle / n ) )

[1] 0.3846546

alpha <- 0.05

print( mle.interval <- c( lambda.hat.mle - qnorm( 1 - alpha / 2 ) *
  se.hat, lambda.hat.mle + qnorm( 1 - alpha / 2 ) * se.hat ) )

[1] 1.317519 2.825338