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