/* SAS program to accompany the article 
   "Two simple ways to construct a log-likelihood function in SAS"
   by Rick Wicklin, published 12JUN2017 on The DO Loop blog
   http://blogs.sas.com/content/iml/2017/06/12/log-likelihood-function-in-sas.html

   For the theory for likelihood and log-likelihood, see
   https://onlinecourses.science.psu.edu/stat504/node/27
*/

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

proc iml;
/* Method 1: Use LOGPDF. This method works in DATA step as well */
start Binom_LL1(p) global(x, NTrials);
   LL = sum( logpdf("Binomial", x, p, NTrials) );
   return( LL );
finish;
 
/* visualize log-likelihood function, which is a function of p */
NTrials = 10;    /* number of trials (fixed) */
use Binomial; read all var "x"; close;
 
p = do(0.01, 0.99, 0.01);      /* vector of parameter values */
LL = j(1, ncol(p), .);
do i = 1 to ncol(LL);
   LL[i] = Binom_LL1( p[i] );  /* evaluate LL for a sequence of p */
end;

title "Graph of Log-Likelihood Function";
title2 "Binomial Distribution, NTrials=10";
call series(p, LL) grid={x y} xvalues=do(0,1,0.1)
                   label={"Probability of Sucess (p)", "Log Likelihood"};

/* Method 2: Manually compute log likelihood by using formula
   PDF(x; p, NTrials) = comb(NTrials,x) # p##x # (1-p)##(NTrials-x)
*/
start Binom_LL2(p) global(x, NTrials);
   LL = sum(lcomb(NTrials, x)) + log(p)*sum(x) + log(1-p)*sum(NTrials-x);
   return( LL );
finish;

LL2 = Binom_LL2(p);
print (max(abs(LL-LL2)));
quit;

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

proc iml;
/* Method 1: use LOGPDF */
start LogNormal_LL1(param) global(x);
   mu = param[1];
   sigma = param[2];
   LL = sum( logpdf("Lognormal", x, mu, sigma) );
   return( LL );
finish;

/* Method 2: Manually compute log likelihood by using formula
   PDF(x; p, NTrials) = comb(NTrials,x) # p##x # (1-p)##(NTrials-x)
*/
start LogNormal_LL2(param) global(x);
   mu = param[1];
   sigma = param[2];
   twopi = 2*constant('pi');
   LL = -nrow(x)/2*log(twopi*sigma##2) 
        - 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;

/* graph of log-likelihood function for lognormal data */
param = expandgrid( do(1,3,0.1), do(0.25, 0.6, 0.025) );
LL = j(nrow(param), 1, .);
do i = 1 to nrow(LL);
   LL[i] = LogNormal_LL2( param[i,] );  /* evaluate LL for a sequence of p */
end;

Z = param || LL;
create graph from Z[c={"mu" "sigma" "LogLik"}]; 
append from Z;
close;
quit;

proc template;
define statgraph ContourPlotParm;
dynamic _X _Y _Z _TITLE;
begingraph;
   entrytitle _TITLE;
   layout overlay;
      contourplotparm x=_X y=_Y z=_Z /
        contourtype=fill nhint=18  colormodel=twocolorramp name="Contour";
      continuouslegend "Contour" / title=_Z;
   endlayout;
endgraph;
end;
run;

proc sgrender data=graph template=ContourPlotParm;
dynamic _TITLE="Graph of Log-Likelihood Function for Lognormal Distribution"
        _X="mu" _Y="sigma" _Z="LogLik";
run;