--------------- Maple ----------------------

ll := ( c, sigma2, n, s ) -> c - n * log( sigma2 ) / 2 - 
  s / ( 2 * sigma2 );

plot( ll( 0, sigma2, 672, 128889.6 ), sigma2 = 160 .. 240 );

solve( diff( ll( c, sigma2, n, s ), sigma2 ) = 0, sigma2 );

evalf( 128889.6 / 672 );

evalf( ll( 0, 191.8, 672, 128889.6 ) );

plot( ll( 2102.168262, sigma2, 672, 128889.6 ), sigma2 = 160 .. 240 );

l := ( c, sigma2, n, s ) -> exp( ll( c, sigma2, n, s ) );

plot( l( 2102.168262, sigma2, 672, 128889.6 ), sigma2 = 160 .. 240 );

sichisq := ( theta, nu, s2 ) -> exp( ( nu / 2 ) * log( nu / 2 ) + 
  ( nu / 2 ) * log( s2 ) - ( 1 + nu / 2 ) * log( theta ) - 
  nu * s2 / ( 2 * theta ) - lnGAMMA( nu / 2 ) );

plot( sichisq( sigma2, 670, 191.8 * 672 / 670 ), 
  sigma2 = 160 .. 240 );

evalf( sichisq( 191.8, 670, 191.8 * 672 / 670 ) );

evalf( sichisq( 160, 670, 191.8 * 672 / 670 ) );

----------------- switch over to R now --------------------

sichisq <- function( theta, nu, s2 ) {

  exp( ( nu / 2 ) * log( nu / 2 ) + ( nu / 2 ) * log( s2 ) - 
    ( 1 + nu / 2 ) * log( theta ) - nu * s2 / ( 2 * theta ) - 
    lgamma( nu / 2 ) )

}

sichisq( 191.8, 670, 191.8 * 672 / 670 )

sigma2.grid = seq( 160, 240, length = 500 )

plot( sigma2.grid, sichisq( sigma2.grid, 670, 191.8 * 672 / 670 ),
  type = 'l', xlab = 'sigma2', ylab = 'Density', lty = 2 )

lines( sigma2.grid, sichisq( sigma2.grid, 672, 191.8 ), lty = 1 )

lines( sigma2.grid, 2.5 / sigma2.grid, lty = 3 )