OPTIONS LINESIZE=75 NOCENTER;

* Calculate the absolute value of r(lambda) for different values;
* of lambda from 0.1 to 10.00 in steps of 0.1. ;

DATA grid_search;
  DO i = 1 to 1000;
    lambda = i/100;
    abs_r_lambda = ABS(lambda - (168/20) * (1 - EXP(-lambda)));
    OUTPUT;
  END;
  DROP i;

PROC PRINT DATA = grid_search NOOBS;
  TITLE 'Grid Evaluation';

RUN;

* Sort the data from the smallest absolute value of r(lambda) to the largest. ;

PROC SORT DATA = grid_search OUT = grid_search_sorted;
  BY abs_r_lambda;

PROC PRINT DATA = grid_search_sorted NOOBS;
  TITLE 'Sorted Grid Values';

RUN;

* The approximate maximum likelihood estimate is the first row of the ;
* sorted dataset (the absolute value of r(lambda) closest to zero. ;

DATA grid_mle;
  SET grid_search_sorted;
  IF (_N_ EQ 1);

* Print out the approximate MLE.;

PROC PRINT DATA = grid_mle NOOBS;
  TITLE 'MLE of Lambda by Grid Search';

RUN;

DATA newton_simple;
  
  /* a good starting guess for lambda is the sample mean */
  lambda = (168/20);

  /* set the previous value of lambda to 0 */


  previous_lambda = 0;

  /* thes variables store the current value of r(lambda) and the
     derivative r’(lambda) */

  r_lambda = 0;
  r_primed_lambda = 0;
  
  /* set the format of the variables - number are displayed up to 12
     characters long, with 10 decimal places. */

  FORMAT lambda 12.10 previous_lambda 12.10 r_lambda 12.10 
         r_primed_lambda 12.10;

  /* now update the Newton-Raphson step 10 times */

  DO iteration = 1 to 10;
    r_lambda = (lambda - (168/20) * (1 - EXP(-lambda)));
    r_primed_lambda = (1 - (168/20) * EXP(-lambda));
    previous_lambda = lambda;
    lambda = lambda - r_lambda/r_primed_lambda;
    OUTPUT;
  END;

/* Print out the resulting dataset */

PROC PRINT DATA=newton_simple NOOBS;
  ID iteration previous_lambda;
  TITLE 'MLE of Lambda by Newton-Raphson - Simple Approach';

RUN;

* A better approach for the problem ;
  
%LET meanx = (168/20);
%LET maxiter = 10;
%LET converge = 1e-8;

DATA newton_better;

  /* Initialize Newton scheme */

  lambda = &meanx;
  previous_lambda = 0;
  r_lambda = 0;
  r_primed_lambda = 0;
  iteration = 1;
  
  /* Set formating */
  
  FORMAT lambda 12.10 previous_lambda 12.10 r_lambda 12.10 
         r_primed_lambda 12.10;

  DO WHILE ((ABS(lambda-previous_lambda) > &converge) & 
              (iteration < &maxiter));
    r_lambda = (lambda - &meanx * (1 - EXP(-lambda)));
    r_primed_lambda = (1 - &meanx * EXP(-lambda));
    previous_lambda = lambda;
    lambda = lambda - r_lambda/r_primed_lambda;
    OUTPUT;
    iteration = iteration + 1;
  END;
  
PROC PRINT DATA=newton_better NOOBS;
  ID iteration previous_lambda;
  TITLE 'MLE of Lambda by Newton-Raphson - Better Approach';
  
RUN;