scaled.inverse.chisq.density <- 
  function( nu.0, sigma.0.squared, sigma.squared ) {

  log.density <- ( nu.0 / 2 ) * log( nu.0 / 2 ) - lgamma( nu.0 / 2 ) +
    ( nu.0 / 2 ) * log( sigma.0.squared ) - ( 1 + ( nu.0 / 2 ) ) *
    log( sigma.squared ) - nu.0 * sigma.0.squared / 
    ( 2 * sigma.squared )

  return( exp( log.density ) )

}

sigma.squared.grid <- seq( 160, 240, length = 500 )

n <- 672

sigma.hat.squared <- 191.8

plot( sigma.squared.grid, scaled.inverse.chisq.density( n,
  sigma.hat.squared, sigma.squared.grid ), type = 'l', lwd = 2,
  xlab = 'sigma squared', ylab = 'Density', 
  main = 'Problem 2 Part (c)' )

lines( sigma.squared.grid, 2.5 / sigma.squared.grid, lwd = 2,
  lty = 2, col = 'red' )

lines( sigma.squared.grid, scaled.inverse.chisq.density( n - 2,
  n * sigma.hat.squared / ( n - 2 ), sigma.squared.grid ),
  lwd = 2, lty = 3, col = 'blue' )

text( 225, 0.015, 'prior', col = 'red', cex = 1.25 )

text( 170, 0.03, 'posterior', col = 'black', cex = 1.25 )

text( 215, 0.02745932, 'likelihood', col = 'blue', cex = 1.25 )

######################################################################

dpareto <- function( alpha, beta, theta ) {

  density <- 0

  if ( theta >= beta ) { density <- alpha * ( beta^alpha ) *
    theta^( - ( alpha + 1 ) ) }

  return( density )

}

n <- 11

m <- 5.1

dpareto( n - 1, m, 6 )

n.grid <- 500

theta.grid <- seq( 0, 10, length = n.grid )

likelihood.grid <- rep( NA, n.grid )

for ( i in 1:n.grid ) {

  likelihood.grid[ i ] <- dpareto( n - 1, m, theta.grid[ i ] )

}

plot( theta.grid, likelihood.grid, type = 'l', lwd = 2,
  col = 'black', xlab = 'theta', ylab = 'Density' )