Newton-Raphson Algorithm in SAS
Find the Bayesian posterior MLE estimate in SAS
Introduction:
In numerical analysis,
Newton's method(also known as the Newton
-Raphson method),
named after
Isaac Newton and
Joseph Raphson,
is a method for finding successively better approximations to the roots (or zeroes) of a real-valued function. How can we apply that in SAS?
Where or when we can apply? Find the Bayesian posterior estimate in SAS might be a good place.
Example: In this tutorial, we want to apply the Newton Raphson algorithm to find the Bayesian posterior MLE estimate.
Based on
previous tutorial,
on Bayesian analysis, the conjugate prior for binormial
distribution is beta distribution. Assume the batting probability(p)) of a group of players
follows beta distribution with parameters (a,b), to find the MLE estimate of those parameters, we can apply the Newton-Raphson algorithm.
Given the batting results of the group of players, say, player i, hiting K_i out of total N_i battings; the corresponding batting probability for
player i is: K_i/N_i. We can simply
apply univariate procedure to get the beta estimators, just as in the previous tutorial.
In that case, 2/20 will have the same effect as 20/200 since they are both p=0.2. To make good use of those K_i and N_i, we have to write down the likelihood
function, find the MLE of (a,b).
%macro Newton_beta(datain,n,s);
** get starting values for (a) and (b) using the mean and variance *;
data Datain;
set &datain(where=(&s<=&n and &n>0));
rate = &s/&n;
run;
proc means data=Datain(where=(&n>0)) noprint;
var rate;
output out=Dataout mean=mean stddev=stddev;
run;
data _null_;
set Dataout;
if stddev**2>mean*(1-mean) then stddev=sqrt(0.9*mean*(1-mean));
a = mean*((mean/(stddev**2))*(1-mean)-1);
b = a*(1-mean)/mean;
call symput("a",a);
call symput("b",b);
run;
Acknowledgement: The original SAS code is from the work by Dr. Speights.
Related links:
Continue to next:
Find the Bayesian posterior MLE estimate in SAS
SAS tutorial home
Back to:
ODS output options: NameLen: NO truncation in the model output
Statistics tutorial home