/* 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;