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;
 

Related links:
Continue to next: Mystery of Proc Sort by Nodup and Nodupkey   SAS tutorial home
Back to: Newton-Raphson Algorithm in SAS   Statistics tutorial home