# contour and perspective plots of the likelihood
# and log likelihood functions in take-home test 2,
# problem 2(B)

# unbuffer the output

# set Your working directory to where You want
# the PDF file of the plots below to be stored

# enter the data:

y <- c( 2.77, 2.50, 1.84, 2.56, 2.32, -1.15 )

V <- c( 1.65, 1.31, 2.34, 1.67, 1.98, 0.90 )^2

# here's a function that computes the log likelihood:

aspirin.mortality.log.likelihood <- function( mu, sigma, y, V ) {

  ll <- ( - 1 / 2 ) * sum( log( V + sigma^2 ) + ( y - mu )^2 /
    ( V + sigma^2 ) )

  return( ll )

} 

n.grid <- 50

# left and right enpoints of mu and sigma grids (see rationale
# for these choices in take-home test 2):

mu.grid <- seq( 0.86 - 4 * 0.59, 0.86 + 4 * 0.59, length = n.grid )

sigma.grid <- seq( 0, 3 * sd( y ), length = n.grid )

ll.grid <- matrix( NA, n.grid, n.grid )

# the next 'for' loops fill in the log-likelihood grid:

for ( i in 1:n.grid ) {

  for ( j in 1:n.grid ) {

    ll.grid[ i, j ] <- aspirin.mortality.log.likelihood( mu.grid[ i ], 
      sigma.grid[ j ], y, V )

  }

}

# You can of course choose Your favorite colors in what follows,
# remembering to ensure that the resulting plots are readable;
# black (the default color) is always a good choice

# the next command makes a contour plot of the log-likelihood function:

contour( mu.grid, sigma.grid, ll.grid, nlevels = 50,
  xlab = 'mu', ylab = 'sigma', main = 'Log Likelihood', col = 'red' )

# stop here and study this plot; when you're ready,
# go on to the next plot, which examines the same function
# with a perspective plot:

persp( mu.grid, sigma.grid, ll.grid, xlab = 'mu', 
  ylab = 'sigma', zlab = 'Log Likelihood', theta = 30, phi = 30,
  col = 'blue' )

# stop here and study this plot

# theta (different meaning than in the aspirin case study) and phi
# control the angles that define Your point of view
# when You look at the perspective plot

# You can experiment with theta and phi (they're in degrees, 
# thinking of a circle as defining 360 degrees) to control Your
# vantage point; i recommend moving them up and down
# in increments of 15 degrees, one by one, to see what they do

# since likelihood = exp( log( likelihood ) ), we don't need
# to write a new function to make contour and perspective plots
# of the likelihood function:

contour( mu.grid, sigma.grid, exp( ll.grid ), nlevels = 50,
  xlab = 'mu', ylab = 'sigma', main = 'Likelihood', col = 'red' )

# stop here and study this plot; when you're ready,
# go on to the next plot, which examines the likelihood function
# with a perspective plot:

persp( mu.grid, sigma.grid, exp( ll.grid ), xlab = 'mu', 
  ylab = 'sigma', zlab = 'Likelihood', theta = 30, phi = 30, 
  col =  'blue' )

# now, to save space, put all 4 plots into a single graph
# for inclusion in Your solutions:

# insert Your choice of file name in place of 'foo' in the next command:

pdf( file = 'foo.pdf' )

par( mfrow = c( 2, 2 ) )

contour( mu.grid, sigma.grid, ll.grid, nlevels = 50,
  xlab = 'mu', ylab = 'sigma', main = 'Log Likelihood', col = 'red' )

# here and in the second 'persp' command below, You may wish
# to change the theta and phi values to ones You liked better
# in Your earlier exploration

persp( mu.grid, sigma.grid, ll.grid, xlab = 'mu', 
  ylab = 'sigma', zlab = 'Log Likelihood', theta = 30, phi = 30,
  col = 'blue' )

contour( mu.grid, sigma.grid, exp( ll.grid ), nlevels = 50,
  xlab = 'mu', ylab = 'sigma', main = 'Likelihood', col = 'red' )

persp( mu.grid, sigma.grid, exp( ll.grid ), xlab = 'mu', 
  ylab = 'sigma', zlab = 'Likelihood', theta = 30, phi = 30,
  col = 'blue' )

par( mfrow = c( 1, 1 ) )

dev.off( )