EZ Study

Actuarial Biology Chemistry Economics Calculators Confucius Engineer
Physics
C.S.

Newton-Raphson Algorithm in SAS
Find the Bayesian posterior MLE estimate in SAS

(Continued) from the previous tutorial on the macro:
``````*iterate using Newton Raphson technique to maximize likelihood *;

%do i=1 %to 100;

data tester;
set Datain end=eof;

**compute the likelihood and 1st and 2nd derivatives of (a) and (b)*;

retain lik d_a d_b d_a2 d_b2 d_ab 0;
sum1_1=0;
sum1_2=0;
sum1_3=0;
do i=0 to &s-1;
sum1_1=sum1_1+log(&a+i);
sum1_2=sum1_2+1/(&a+i);
sum1_3=sum1_3+1/((&a+i)**2);
end;
sum2_1=0;
sum2_2=0;
sum2_3=0;
do i=0 to &n-&s-1;
sum2_1=sum2_1+log(&b+i);
sum2_2=sum2_2+1/(&b+i);
sum2_3=sum2_3+1/((&b+i)**2);
end;
sum3_1=0;
sum3_2=0;
sum3_3=0;
do i=0 to &n-1;
sum3_1=sum3_1+log(&a+&b+i);
sum3_2=sum3_2+1/(&a+&b+i);
sum3_3=sum3_3+1/((&a+&b+i)**2);
end;

lik=lik+sum1_1+sum2_1-sum3_1;
d_a=d_a + sum1_2 - sum3_2;
d_b=d_b + sum2_2 - sum3_2;
d_a2=d_a2 - sum1_3 + sum3_3;
d_b2=d_b2 - sum2_3 + sum3_3;
d_ab=d_ab + sum3_3;

if eof then do;

a=&a;
b=&b;

det = d_a2*d_b2 - d_ab*d_ab;
a_new = &a - (d_a*d_b2 - d_b*d_ab)/det;
b_new = &b - (d_b*d_a2 - d_a*d_ab)/det;

error = sqrt(d_a*d_a+d_b*d_b);

call symput("a",a_new);
call symput("b",b_new);
if error<.001 then call symput("i",101);
iteration=&i;
mean = a_new/(a_new+b_new);
std_dev = sqrt(mean*b_new/((a_new+b_new)*(a_new+b_new+1)));
output;
end;

run;
%if &i=1 %then %do;
data Output1;
set tester;
run;
%end;
%else %do;
proc append base=Output1 data=tester;
run;
%end;

%end;
proc print data=Output1;
var iteration error mean std_dev a_new b_new;
run;

%mend Newton_beta;
``````