Friday, December 17, 2010

Trapezoidal Rule

TAD (Time Average Difference) is area under the curve (AUC).

Trapezoidal Rule

Using Trapezoidal Rule for the Area Under a Curve Calculation

Estimate the Area Under a ROC Curve

SAS Calculations of AUC for Multiple Metabolic Readings



/************************************************************************
AREA.SAS

DISCLAIMER: THIS INFORMATION IS PROVIDED BY SAS INSTITUTE INC. AS A SERVICE TO ITS USERS. IT IS PROVIDED "AS IS". THERE ARE NO WARRANTIES, EXPRESSED OR IMPLIED, AS TO MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE REGARDING THE ACCURACY OF THE MATERIALS OR CODE CONTAINED HEREIN.

PURPOSE: This program uses PROC EXPAND to calculate the approximate area under the curve for some sample data. The sample data should consist of (x,y) pairs.

DETAILS: For this example, the sample data is generated from a high degree polynomial. PROC EXPAND is then used to compute the approximate area under the curve using each of the following methods:

a. Cubic Spline interpolation. b. Trapezoid rule.

The exact area, given by the definite integral, is calculated for the polynomial curve in order to assess the precision of the approximations.
************************************************************************/

%let lower=-2; %let upper=1; %let interval=0.2;

* generate some data according to a high order polynomial; data kvm; do x=&lower to &upper by &interval; y=15+(x-2)*(x-1.5)*(x-1)*(x-.5)*x*(x+.5)*(x+1)*(x+1.5)*(x+2); output; end;

proc sort; by x;

/* PROC EXPAND will include a contribution for the last interval. For an accurate approximation to the integral, we need to make sure that this last contribution is negligible. So we'll append an additional x value which is extremely close to the last x value. Of course, the two Y values will be identical. But the result is that the last interval is extremely short, so any contribution to the integral approximation is negligible. */

proc print data=kvm(obs=50); title 'First few observations of the original data'; run;

data one; set kvm end=eof; output; if eof then do; x=x+(1e-10); output; end; run;

proc print data=one(obs=50); title 'First few observations of the original data'; run;

proc gplot data=one; title 'original series'; plot y*x; run;

proc expand data=one out=three method=spline ; convert y=total/observed=(beginning,total) transformout=(sum); id x; run;

proc sort data=three; by descending total;

proc print data=three(obs=1); var total; title 'Approximate Integral Using Spline method'; run;

*************************************************;

proc expand data=one out=three method=join; convert y=total/observed=(beginning,total) transformout=(sum); id x; run;

proc sort data=three; by descending total;

proc print data=three(obs=1); title 'Approximate Integral Using Trapezoid Rule'; var total; run;

***********************************************;

/* Since this data was generate using a high order polynomial, it's easy to compute the definite integral value. Then we can compare our approximations and assess their precision. */

data four; set one; x3= x*x*x; x5= x*x*x*x*x; x7= x*x*x*x*x*x*x; x9= x*x*x*x*x*x*x*x*x; run;

proc autoreg data=four outest=five noprint; title 'Calculate polynomial coefficients'; model y=x x3 x5 x7 x9; run;

proc print data=five; run;

data six; set five; title 'brute force (hand) integration of polynomial'; val=&upper; eval=intercept*val+x*(val**2)/2+x3/4*(val**4)+x5/6*(val**6) +x7/8*(val**8)+x9/10*(val**10); temp=eval; val=&lower; eval=intercept*val+x*(val**2)/2+x3/4*(val**4)+x5/6*(val**6) +x7/8*(val**8)+x9/10*(val**10); integral=temp-eval; keep integral; run;

proc print; run; title;