proc iml;
use Sashelp.Iris;
read all var {"SepalWidth"} into x;
close Sashelp.Iris;

/* write the log-likelihood function for normal dist */
start LogLik(param) global (x);
   mu = param[1];
   sigma2 = param[2]##2;
   n = nrow(x);
   c = x - mu;
   return( -n/2*log(sigma2) - 1/2/sigma2*sum(c##2) );
finish;

con = { .   0,  /* lower bounds: -infty < mu; 0 < sigma */
        .   .}; /* upper bounds:  mu < infty; sigma < infty */
p = {35.5 4.95};/* initial guess for solution */
optn = {1,  /* find max of function       */
        3}; /* print a lot of output      */
call nlpnra(rc,xres,"LogLik", p, optn, con);

params = Define2DGrid( 27, 36, 21, 
                       3,   6, 21 );
f = j(nrow(params),1);
do i = 1 to nrow(params);
   f[i] = LogLik( params[i,] );
end;

VarName = "SepalWidth";
declare Histogram hist;
hist = Histogram.Create("Distrbution", x);
hist.ShowDensity();
hist.SetAxisLabel(XAXIS, VarName);

declare DataObject dobj;
dobj = DataObject.Create("LogLiklihood", {"mu" "sigma" "LogLik"}, params||f);
declare RotatingPlot rot;
rot = RotatingPlot.Create(dobj, "mu", "sigma", "LogLik");
rot.SetDrawingMode(SMOOTHCOLORMESH);
rot.ShowFrame( true );
rot.ShowReferenceLines();
rot.SetWindowPosition(50,50);
rot.SetTitleText("Log-Likelihood Function for Normal Distribution", true);
pause "NoDialog:";


@ods output grad=grad IterHist=iterhist;
call nlpnra(rc,xres,"LogLik", p, optn, con);
@ods output close;
pause "NoDialog:";

submit;
data a;
merge iterhist(rename=(OptCrit=LogLik DifCrit=DifLogLik))
grad(rename=(x1=Mu x2=Sigma));
run;
/*
proc print data=a noobs;
var Iter Mu Sigma LogLik DifLogLik MaxGrad;
run;
*/
endsubmit;

use grad;
read all var {x1 x2} into path;
close grad;
path = p // path;

declare ContourPlot cont;
cont = ContourPlot.Create(dobj, "mu", "sigma", "LogLik");
cont.DrawUseDataCoordinates();
cont.DrawSetPenWidth(2);
cont.SetWindowPosition(50,50);
hist.DrawUseDataCoordinates();
hist.DrawSetPenColor(GREY);
t = do(min(x), max(x), (max(x)-min(x))/100)`;
do i = 1 to nrow(path);
   mu = path[i,1];
   sigma = path[i,2];
   cont.DrawMarker( mu, sigma, MARKER_CIRCLE, 7 );
   cont.SetTitleText( "Iteration ="+char(i-1,2), true);
   if i>1 then 
      cont.DrawLine(path[i-1,1], path[i-1,2], mu, sigma);
   
   hist.DrawLine(t, pdf("Normal", t, mu, sigma));
   pause "NoDialog:";
end;
cont.DrawMarker( mu, sigma, MARKER_STAR, 8 );
hist.DrawSetPenAttributes(BLACK, SOLID, 3);
hist.DrawLine(t, pdf("Normal", t, mu, sigma));