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

hist( y, breaks = ( -0.5 ):6.5, prob = T )

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

  likelihood <- lambda^s * exp( - n * lambda )

  return( likelihood )

}

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

  log.likelihood <- s * log( lambda ) - n * lambda

  return( log.likelihood )

}

print( n <- length( y ) )

# [1] 14

print( s <- sum( y ) )

# [1] 29

lambda.grid <- seq( 0.5, 4, length = 500 )

plot( lambda.grid, poisson.likelihood( lambda.grid, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'lambda',
  ylab = 'likelihood', main = 'Poisson sampling model' )

plot( lambda.grid, poisson.log.likelihood( lambda.grid, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'lambda',
  ylab = 'log likelihood', main = 'Poisson sampling model' )

print( lambda.hat.mle <- s / n )

# [1] 2.071429
 
print( se.hat.lambda.hat.mle <- sqrt( lambda.hat.mle / n ) )

# [1] 0.3846546

# for 99.9% confidence use 3.2 as the z-value,
# not 1.96 as with 95% confidence

# 0.999 probability in the middle means 
# 0.0005 in the right tail and therefore 0.9995
# to the left:

qnorm( .9995 )

# [1] 3.290527

print( lambda.interval <- c( lambda.hat.mle - qnorm( 0.9995 ) *
  se.hat.lambda.hat.mle, lambda.hat.mle + qnorm( 0.9995 ) *
  se.hat.lambda.hat.mle ) )

# [1] 0.8057122 3.3371449

# we're 99.9% confident that the population mean
# length of stay is between about 0.8 days and 3.3 days