# likelihood analysis of case study 2: AMI mortality

# unbuffer the output (if necessary)

# change the next line to the directory path of your choice
# (where you want any output, such as plots, to go)

# setwd( 'C:/DD/Teaching/AMS-STAT-206/Winter-2020' )

n <- 403

s <- 72

bernoulli.likelihood <- function( theta, n, s ) {

  likelihood <- theta^s * ( 1 - theta )^( n - s )

  return( likelihood )

}

bernoulli.likelihood( 0.2, n, s )

# [1] 3.953169e-83

bernoulli.likelihood( 0.2, 10 * n, 10 * s )

# [1] 0

# R has just lied to us: the actual value is 
#   9.32081356715980 × 10^( -825 )

# wolfram alpha code:

#   evaluate theta^72 * ( 1 - theta )^( 403 - 72 ) at theta = 0.2

# in R we've just experienced an instance of underflow:
# R considers any value between -2.225074e-308
# and 2.225074e-308 to be 0

# i got the 2.225074e-308 value by asking R about
# the built-in object .Machine

.Machine

# $double.eps
# [1] 2.220446e-16

# $double.neg.eps
# [1] 1.110223e-16

# $double.xmin
# [1] 2.225074e-308

# $double.xmax
# [1] 1.797693e+308

#  .
#  .
#  .

# with ( n, s ) = ( 403, 72 ) we're fine in R (no underflow)

theta.grid.1 <- seq( 0, 1, length = 500 )

plot( theta.grid.1, bernoulli.likelihood( theta.grid.1, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'theta',
  ylab = 'Bernoulli likelihood' )

theta.grid.2 <- seq( 0.1, 0.3, length = 500 )

plot( theta.grid.2, bernoulli.likelihood( theta.grid.2, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'theta',
  ylab = 'Bernoulli likelihood' )

bernoulli.log.likelihood <- function( theta, n, s ) {

  log.likelihood <- s * log( theta ) + ( n - s ) * log( 1 - theta )

  return( log.likelihood )

}

plot( theta.grid.1, bernoulli.log.likelihood( theta.grid.1, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'theta',
  ylab = 'Bernoulli log likelihood' )

plot( theta.grid.2, bernoulli.log.likelihood( theta.grid.2, n, s ),
  type = 'l', lwd = 2, col = 'red', xlab = 'theta',
  ylab = 'Bernoulli log likelihood' )

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

# [1] 0.17866

# comparing log likelihood functions at the MLE for increasing n,
# keeping s / n fixed at theta.hat.mle:

print( maximum.log.likelihood.for.n.of.403 <- 
  bernoulli.log.likelihood( theta.hat.mle, n, s ) )

# [1] -189.1503

plot( theta.grid.2, bernoulli.log.likelihood( theta.grid.2, n, s ) -
  maximum.log.likelihood.for.n.of.403, type = 'l', lwd = 2, 
  col = 'red', xlab = 'theta', ylab = 'Bernoulli log likelihood' )

text( 0.25, -1.0, 'n = 403', col = 'red' )

print( maximum.log.likelihood.for.n.of.4030 <- 
  bernoulli.log.likelihood( theta.hat.mle, 10 * n, 10 * s ) )

# [1] -1891.503

lines( theta.grid.2, 
  bernoulli.log.likelihood( theta.grid.2, 10 * n, 10 * s ) -
  maximum.log.likelihood.for.n.of.4030, lwd = 2, 
  col = 'blue', lty = 2 )

text( 0.23, -13.0, 'n = 4030', col = 'blue' )