/* 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;