y_bar := 404.59; s_u := 6.466846; n := 100; likelihood := ( mu, sigma ) -> sigma^( - n ) * exp( - ( ( n - 1 ) * s_u^2 + n * ( y_bar - mu )^2 ) / ( 2 * sigma^2 ) ); plot3d( likelihood, 402.5 .. 406.5, 5.0 .. 8.0 ); log_likelihood := ( mu, sigma ) -> - n * log( sigma ) - ( ( n - 1 ) * s_u^2 + n * ( y_bar - mu )^2 ) / ( 2 * sigma^2 ); plot3d( log_likelihood, 402.5 .. 406.5, 5.0 .. 8.0 ); ################# new session assume( s_MLE >= 0, n > 0 ); log_likelihood := ( mu, sigma ) -> - n * log( sigma ) - ( n * s_MLE^2 + n * ( y_bar - mu )^2 ) / ( 2 * sigma^2 ); simplify( diff( log_likelihood( mu, sigma ), mu ) ); simplify( diff( log_likelihood( mu, sigma ), sigma ) ); diff( log_likelihood( mu, sigma ), sigma ); simplify( diff( log_likelihood( mu, sigma ), mu$2 ) ); diff( log_likelihood( mu, sigma ), sigma$2 ); simplify( diff( diff( log_likelihood( mu, sigma ), mu ), sigma ) ); Hessian := simplify( Matrix( [ [ diff( log_likelihood( mu, sigma ), mu$2 ), diff( diff( log_likelihood( mu, sigma ), mu ), sigma ) ], [ diff( diff( log_likelihood( mu, sigma ), mu ), sigma ), diff( log_likelihood( mu, sigma ), sigma$2 ) ] ] ) ); Information_Matrix := simplify( eval( - Hessian, [ mu = y_bar, sigma = s_MLE ] ) ); Covariance_Matrix := LinearAlgebra[ MatrixInverse ]( Information_Matrix ); simplify( sqrt( Covariance_Matrix[ 1, 1 ] ) ); simplify( sqrt( Covariance_Matrix[ 2, 2 ] ) );