# the 'plotly' code below was written by erdong guo

# some code in this file is from Professor Draper's contour and 
# perspective plots of the likelihood and log likelihood functions 
# in take-home test 2, problem 2(B), available on the course web page

# reference for package 'plotly': https://plot.ly/r/#fundamentals

# this is our first introduction to the current best practice
# in making plots in R: an approach created by hadley wickham
# based on a CRAN package called 'ggplot2'

install.packages( 'plotly' )

# a window will appear; select a mirror site:
# 
#   USA (CA 1) [https]
#
# is a good choice for the mirror site; and 'plotly' will be installed

library( plotly )

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

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 two nested '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 )

  }

}

# we transpose the log likelihood matrix to get ready for 'plot_ly'

ll.grid.t <- t( ll.grid )

# code to create interactive log likelihood contour plot

interactive.log.likelihood.contour.plot <- 
  plot_ly( x = mu.grid, y = sigma.grid, z = ll.grid.t, 
  type = 'contour', contours = list( showlabels = T ), 
  colorscale = 'Viridis' ) %>% colorbar( title = 'Log Likelihood') %>%
  layout( xaxis = list( title = 'mu' ), yaxis = list( title = 'sigma' ),
  title = 'Contours (level sets) of the Log likelihood function' )

# notice how this approach to plotting breaks the different parts
# of the final plot into separate steps (in this case 'plot_ly'
# and 'layout', connected together with R's pipe operator '%>%',
# which behaves like the pipe '>' in unix and linux

# if you get into 'plot_ly', here's a list of its 'colorscales'
# for you to play with: ‘Blackbody’, ‘Bluered’, ‘Blues’, ‘Earth’,
# ‘Electric’, ‘Greens’, ‘Greys’, ‘Hot’, ‘Jet’, ‘Picnic’, ‘Portland’,
# ‘Rainbow’, ‘RdBu’, ‘Reds’, ‘Viridis’, ‘YlGnBu’, and ‘YlOrRd’

# although when i played with this in the contour plot above
# it didn't seem to make any difference

# the next command will show the log likelihood contour plot in your browser;
# please allow the blocked content if it's blocked

interactive.log.likelihood.contour.plot

# you can try clicking inside the plot, dragging the mouse arrow, ...

# to check this result, we can also make a non-interactive contour plot
# and verify that it's providing the same information about the maximum

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

# code to create interactive log-likelihood surface plot 
# (this is a two-dimensional hypersurface embedded in 
# three-dimensional Euclidean space)

interactive.log.likelihood.perspective.plot <- 
  plot_ly( x = mu.grid, y = sigma.grid, z = ll.grid.t ) %>% 
  add_surface( contours = list( z = list( show = T, usecolormap = T,
  highlightcolor = '#ff0000', project = list( z = T ) ) ) ) %>%
  layout( title = 'Perspective plot of the log-likelihood function',
  scene = list( xaxis = list( title = 'mu' ), 
  yaxis = list( title = 'sigma' ), zaxis = list( title = 'loglikehood' ),
  camera = list( eye = list( x = 1.5, y = 0.5, z = -0.5 ) ) ) ) %>% 
  colorbar( title = 'Log likelihood' )

# '#ff0000' is a color called 'hex red', also known as 'rgb( 255, 0, 0 )'

# you can also play with this color choice

# the next command will show the log likelihood perspective plot 
# in your browser; please allow the blocked content if it's blocked

interactive.log.likelihood.perspective.plot

# single-left-click and move the mouse arrow to rotate the plot

# to check this result, we can also make a static perspective plot

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

# code used to create interactive likelihood contour plot

interactive.likelihood.contour.plot <- 
  plot_ly( x = mu.grid, y = sigma.grid, z = exp( ll.grid.t ), 
  type = 'contour', contours = list( showlabels = T ), 
  colorscale = 'Jet' ) %>% layout( xaxis = list( title = 'mu' ),
  yaxis = list( title = 'sigma' ), 
  title = 'Contour plot of the likelihood function' ) %>% 
  colorbar( title = 'LogLikelihood' )  

# the next command will show the likelihood contour plot in your browser;
# please allow the blocked content if it's blocked

interactive.likelihood.contour.plot

# to check this plot, we can make a static contour plot

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

# code to create likelihood perspective plot

interactive.likelihood.perspective.plot <- 
  plot_ly( x = mu.grid, y = sigma.grid, z = exp( ll.grid.t ) ) %>% 
  add_surface( contours = list( z = list( show = T, usecolormap = T,
  highlightcolor = '#ff0000', project = list( z = T ) ) ) ) %>%
  layout( title = 'Perspective plot of the likelihood function',
  scene = list( xaxis = list( title = 'mu' ), 
  yaxis = list( title = 'sigma' ),
  zaxis = list( title = 'likehood' ), 
  camera = list( eye = list( x = 1.5, y = 0.5, z = -0.5 ) ) ) ) %>% 
  colorbar( title = 'Likelihood' )

# the next command will show the likelihood perspective plot 
# in your browser; please allow the blocked content if it's blocked

interactive.likelihood.perspective.plot

# to check this, we can make a static perspective plot

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