/* SAS program to accompany the article 
   "Two ways to compute maximum likelihood estimates in SAS"
   by Rick Wicklin, published 14JUN2017 on The DO Loop blog
   http://blogs.sas.com/content/iml/2017/06/14/maximum-likelihood-estimates-in-sas.html

   First read the article
   http://blogs.sas.com/content/iml/2017/06/12/log-likelihood-function-in-sas.html
*/

/********************************************************/
/* Construct log-likelihood function for binomial data  */
/********************************************************/

/* Create example Binomial(p, 10) data for unknown p */
data Binomial;  
input x @@;
datalines;
6 2 6 4 5 7 5 7 6 4 6 4 7 4 7 7 8 5 8 4 
;

/********************************************************/
/* Construct log-likelihood function for lognormal data */
/********************************************************/
/* Create example Lognormal(2, 0.5) data of size N=200 */
%let N = 200;         
data LN;
call streaminit(98765);
sigma = 0.5;               /* shape or log-scale parameter    */
zeta = 2;                  /* scale or log-location parameter */
do i = 1 to &N;
   z = rand("Normal", zeta, sigma);  /* z ~ N(zeta, sigma)    */
   x = exp(z);                       /* x ~ LogN(zeta, sigma) */
   output;
end;
keep x;
run;


/********************************************************************/
/* Method 1: Compute maximum likelihood estimates by using PROC IML */
/********************************************************************/

/* Example 1: MLE for binomial data */
proc iml;
/* log-likelihood function for binomial data */
start Binom_LL(p) global(x, NTrials);
   LL = sum( logpdf("Binomial", x, p, NTrials) );
   return( LL );
finish;

NTrials = 10;    /* number of trials (fixed) */
use Binomial; read all var "x"; close;

/* set constraint matrix, options, and initial guess for optimization */
con = { 0,      /* lower bounds: 0 < p     */
        1};     /* upper bounds:     p < 1 */
opt = {1,       /* find maximum of function   */
       2};      /* print some output      */
p0  = 0.5;      /* initial guess for solution */
call nlpnra(rc, p_MLE, "Binom_LL", p0, opt, con);
print p_MLE;
quit;

/* Example 2: MLE for lognormal data */
proc iml;
/* log-likelihood function for lognormal data */
start LogNormal_LL(param) global(x);
   mu = param[1];
   sigma = param[2];
   /*    LL = sum( logpdf("Lognormal", x, mu, sigma) );  */
   twopi = 2*constant('pi');
   LL = -nrow(x)/2*log(twopi*sigma##2)     /* use '##' to vectorize */
        - sum( (log(x)-mu)##2 ) / (2*sigma##2)
        - sum(log(x));  /* this term is constant w/r/t (mu, sigma) */
   return( LL );
finish;
use LN; read all var "x"; close;

/* constraint matrix */
con = { . 0,       /* lower bounds: 0 < sigma   */
        . .};      /* upper bounds:     N/A     */
opt = {1,          /* find maximum of function   */
       1};         /* print some output      */
param0  = {1 1};   /* initial guess for solution */
call nlpnra(rc, param_MLE, "LogNormal_LL", param0, opt, con);
print param_MLE[c={"mu" "sigma"}];
QUIT;


/************************************************/
/* Method 2: Compute MLEs by using PROC NLMIXED */
/************************************************/
/* Built-in LL for Bernoulli, binomial, Poisson, neg binomial,
                   normal, and gamma
   Not only get MLE, but standard error and 95% CI */
proc nlmixed data=Binomial;
   parms p = 0.5;
   bounds 0 < p < 1;
   NTrials = 10;
   model x ~ binomial(NTrials, p);
run;

proc nlmixed data=LN;
parms mu 1 sigma 1;
bounds 0 < sigma;
/* LL = logpdf("Lognormal", x, mu, sigma); */
sqrt2pi = sqrt(2*constant('pi'));
LL = -log(sigma) 
     - log(sqrt2pi*x)               /* this term is constant w/r/t (mu, sigma) */
     - (log(x)-mu)**2  / (2*sigma**2);
model x ~ general(LL);
run;