/* SAS Support Pages
http://support.sas.com/documentation/cdl/en/imlug/64248/HTML/default/viewer.htm#imlug_r_sect012.htm
http://support.sas.com/documentation/cdl/en/imlug/64248/HTML/default/viewer.htm#imlug_r_sect013.htm
Edit "myuserid" in the setwd function of the R-code
*/

proc iml;
q = {3.7, 7.1, 2, 4.2, 5.3, 6.4, 8, 5.7, 3.1, 6.1, 4.4, 5.4, 9.5, 11.2};
RVar = "rq";
call ExportMatrixToR( q, RVar );

submit RVar / R;
   library(KernSmooth)
   idx <-which(!is.na(&RVar))        # must exclude missing values (NA)
   p <- &RVar[idx]                   #    from KernSmooth functions
   h = dpik(p)                       # Sheather-Jones plug-in bandwidth
   est <- bkde(p, bandwidth=h)       # est has 2 columns
endsubmit;

call ImportMatrixFromR( m, "est" );
/* estimate the density for q >= 8 */
x = m[,1];                /* x values for density */  
idx = loc( x>=8 );        /* find values x >= 8 */
y = m[idx, 2];            /* extract corresponding density values */

/* Use the trapezoidal rule to estimate the area under the density curve.
   The area of a trapezoid with base w and heights h1 and h2 is 
   w*(h1+h2)/2. */
w = m[2,1] - m[1,1];      
h1 = y[1:nrow(y)-1];      
h2 = y[2:nrow(y)];
Area = w * sum(h1+h2) / 2;
print Area;

submit / R;
   setwd('//ahs-sas-appserv/sasusers$/myuserid/My SAS Files/9.4')
   png("RPlot.png", 400, 400)
   hist(p, freq=FALSE) # histogram
   lines(est)          # kde overlay
endsubmit;

submit / R;
   setwd('//ahs-sas-appserv/sasusers$/myuserid/My SAS Files/9.4')
   pdf("RPlot.pdf", 3, 3)
   hist(p, freq=FALSE) # histogram
   lines(est)          # kde overlay
endsubmit;