/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.5 ] */

/* [wxMaxima: comment start ]
   Eric Doviak
   original: 22 April 2012
   updated:  30 April 2012

   Math Methods 7025X
   shape of the likelihood function 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
   First, load the "distrib" package.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load(distrib)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Now, randomly draw 100 values from the standard normal.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
N:100$
x:random_normal(0,1,N)$

mnx:(1/N)*sum(x[i],i,1,N)$
sdx:sqrt((1/N)*sum((x[i]-mnx)^2,i,1,N))$

print("")$
print("mean of x: ",mnx)$
print("std. dev.: ",sdx)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Set up the log-likelihood function.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
loglik(mu,sigma):= -N*log(sigma) - (N/2)*log(2*%pi) - (1/2)*sum(((x[i]-mu)/sigma)^2,i,1,N);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Deriviatives must be taken with respect to sigma^2, so define:

   gamma == sigma^2

   and take derivatives with respect to gamma.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
loglik(mu,gamma):= -(N/2)*log(gamma) - (N/2)*log(2*%pi) - (1/(2*gamma))*sum((x[i]-mu)^2,i,1,N);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Maximize it with respect to mu and gamma. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sol:lbfgs(-loglik(mu,gamma),'[mu,gamma],[0.01,0.99],0.0001,[-1,0])$
mu_max:subst(sol[1],mu)$
gamma_max:subst(sol[2],gamma)$

print("")$
print(mu," = ",mu_max)$
print(sigma^2," = ",gamma_max)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   To see the maximum, plot the log likelihood function in "mu-gamma" space.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
plot3d(loglik(mu,gamma),[mu,-0.5,0.5],[gamma,0.5,1.5])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Check to see if second-order conditions are satisfied.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/* set up the Hessian matrix */
dxx(mu,gamma):=''(diff(diff(loglik(mu,gamma),mu),mu))$
dyy(mu,gamma):=''(diff(diff(loglik(mu,gamma),gamma),gamma))$
dxy(mu,gamma):=''(diff(diff(loglik(mu,gamma),mu),gamma))$
H:matrix(
    [dxx(mu_max,gamma_max),dxy(mu_max,gamma_max)],
    [dxy(mu_max,gamma_max),dyy(mu_max,gamma_max)])$

print("")$
print("own-partials must be negative:")$
print("")$
print("d^2 loglik(mu,gamma)"/"(d mu)^2"," = ",dxx(mu_max,gamma_max))$
print("")$
print("d^2 loglik(mu,gamma)"/"(d gamma)^2"," = ",dyy(mu_max,gamma_max))$
print("")$
print("")$
print("the cross-partial:")$
print("")$
print("d^2 loglik(mu,gamma)"/"d mu d gamma"," = ",dxy(mu_max,gamma_max))$
print("")$
print("")$
print("the Hessian matrix:")$
print("H = ",H)$
print("")$
print("determinant of Hessian must be positive")$
print("det(H) = ",determinant(H))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
info:-1*invert(H)$

print("")$
print("the information matrix:")$
print("-1*(H^-1) = ",info)$
print("")$
print(mu,": ",mu_max,"   se:",sqrt(info[1,1]))$
print(sigma^2,": ",gamma_max,"   se:",sqrt(info[2,2]))$
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$