Showing posts with label proc iml. Show all posts
Showing posts with label proc iml. Show all posts

Thursday, March 15, 2012

Make a frequency function in SAS/IML


Aggregation is probably the most popular operation in the data world. R comes with a handy table() function. Usually in SAS, the FREQ procedure would deal with this job. It will be great if SAS/IML has an equivalent function. I just created a user-defined function or module for such a purpose. Since it contains a DO loop, the efficiency is not very ideal -- always 10 times slower than PROC FREQ for a simulated data set of one million records.

/* 1 - Use IML for simulation and aggregation */
proc iml;
start freq(invec);
x = t(unique(invec));
y = repeat(x, 1, 2);
do i = 1 to nrow(x);
y[i, 2] = ncol(loc(invec=y[i, 1]));
end;
return(y);
finish;
store module = freq;
quit;

proc iml;
load module = freq;
/* Simulate a vector with 1 million values */
test = abs(floor(rannor(1:1e6)*100));
t0 = time();
result = freq(test);
timer = time() - t0;
print timer;
/* Output the result matrix as SAS data set */
create a var{"level" "frequency"};
append from result;
close a;
quit;

proc sgplot data = a;
series x = level y = frequency;
run;

/* 2 - Use PROC FREQ for simulation and aggregation */
data test;
do i = 1 to 1e6;
test = abs(floor(rannor(1234)*100));
output;
end;
drop i;
run;

options fullstimer;
proc freq data = test noprint;
table test / out = b;
run;

Monday, January 23, 2012

A test for memory management of SAS/IML


Programming always involves the considerations for the efficiency and the memory usage. For efficient programming in SAS/IML, my shortcut is to look at the tip sheet from Rick Wicklin and search ways to simplify the codes. As for the memory management mechanism of SAS/IML, I only found one page of SAS/IML 9.2 User’s Guide on the Internet.



To see the performance of SAS/IML ‘s memory management, I designed a simple test, since the SHOW SPACE statement would indicate the memory details of SAS/IML. A simulated 200 rows * 300 columns matrix occupies about 400k memory. I just requested 1MB memory by specifying the WORKSIZE at the beginning, which means that three matrices of such size would blow the tiny work space away. First I assigned multiple references toward this matrix, and then changed one value in this matrix. Finally I cleared all of the objects in the memory and generated a matrix with the same size. The memory change is shown in the plot above.

/* Only request 1MB memory for this test*/
proc iml worksize=1048576;
/* 0 -- State 0*/
show space;
/* 1 -- State 1 (1 reference on 1 object)*/
x = j(200, 300, 0);
call randgen(x, "Normal");
show space;
/* 2 -- State 2 (2 references on 1 object)*/
y = x;
show space;
/* 3 -- State 3 (3 references on 1 object)*/
z = x;
show space;
/* 4 -- State 4 (2 refernces on 1 object and 2 refernces on another object)*/
x[1, 1] = 5;
w = x;
show space;
/* 5 -- State 5 (0 references on 2 objects) */
free x y z w;
show space;
/* 6 -- State 1 again (1 reference on 1 object)*/
x = j(200, 300, 0);
call randgen(x, "Normal");
show space;
quit;

Conclusions:
1. Interestingly, an object with two references costs more memory than the same object with three references. Apparently two identical copies of the same matrix exist at State 2, while three references all point to a single matrix at State 3.

SAS/IML seems to have a unique memory management system. The memory manager of IML first checks the size available of the workspace (memory allocated to IML). If the memory usage doesn’t pass the alarming line, it will generously make copies. Otherwise, it will become very aggressive and start to compress the memory frenetically.

2. Two similar matrices stay in the memory at State 4, although the difference between them is just one cell. Each of the matrices has 2 references.

SAS/IML follows the rule of copy-on-change. Once values in a matrix are changed, the memory manager will make a copy for the changed object.

3. Running the FREE Statement doesn't immediate clear the memory at State 5.

SAS/IML’s FREE statement only kills the references. The unreferenced objects will be discarded during the ensuing operations.

Thursday, January 19, 2012

6 ways to count odd numbers in SAS/IML


SAS/IML has a number of vector-wise subscripts/operators/functions available, which can make many things easy. A cheat sheet about them can be found at Rick Wicklin’s blog.

To try out those wonderful features( and their combinations?), I designed a test to use them for the total of the odd numbers in a random numeric sequence. A typical solution is always placing a cumulative counter in a looping structure for many programming languages. Therefore I just used such basic DO loop in SAS/IML as a benchmark.

At the beginning, the MOD function computes the modulo. Then the simplest ways is the SUM function to aggregate all odd number. Subscript, like [, +], serves the same purpose. The method of the CHOOSE + SUM functions quite assembles the DO loop and is the generalized form of the SUM only method(also brings overhead on resources). The LOC function can subset or index a vector and is combined with two other functions.

Observations:
1. All vector-wise methods beat the DO loop, especially with big dataset. The simpler the method is, the faster the result is.
2. The robust CHOOSE function. For binary conditions, the CHOOSE function can replace the if-else-then statements plus a DO loop and is much more efficient. Rick Wicklin has an article about this function.
3. The LOC function is very handy and plays a role like the WHERE statement in SAS’s DATA step or the which()/subset() functions in R.
4. The SUM function seems slightly faster than the subscript [, +]  for a vector.


proc iml;
a = t(do(1e6, 2e7, 1e6));
timer = j(nrow(a), 6);
do p = 1 to nrow(a);
n = a[p];
/* Simulate a numeric sequence */
x = ceil(ranuni(1:n)*100000);

/* 1 -- SUM function*/
t0 = time();
r1 = sum(mod(x, 2));
timer[p, 1] = time() - t0;

/* 2 -- Subscript + */
t0 = time();
r2 = mod(x, 2)[ , +];
timer[p, 2] = time() - t0;

/* 3 -- SUM + CHOOSE functions*/
t0 = time();
r3 = sum(choose(mod(x, 2), 1, 0));
timer[p, 3] = time() - t0;

/* 4 -- NCOL + LOC functions */
t0 = time();
r4 = ncol(loc(mod(x, 2) = 1));
timer[p, 4] = time() - t0;

/* 5 -- DO loop */
t0 = time();
r5 = 0;
do i = 1 to ncol(x);
if mod(x[i], 2) = 1 then r5 = r5 + 1;
end;
timer[p, 5] = time() - t0;

/* 6 -- COUNTMISS + LOC functions */
t0 = time();
x[loc(mod(x, 2) = 1)] = .;
r6 = countmiss(x);
timer[p, 6] = time() - t0;

/* Validate all results */
print r1 r2 r3 r4 r5 r6;
end;
t = a||timer;
create _1 from t;
append from t;
close _1;
quit;

data _2;
set _1;
length test $100.;
label col1 = "Number of observations"
time = "Time by seconds to count odd numbers";
test = "SUM function"; time = col2; output;
test = "Subscript + "; time = col3; output;
test = "SUM + CHOOSE functions"; time = col4; output;
test = "NCOL + LOC functions"; time = col5; output;
test = "DO loop"; time = col6; output;
test = "COUNTMISS + LOC functions"; time = col7; output;
keep test time col1;
run;

proc sgplot data = _2;
series x = col1 y = time / curvelabel group = test;
yaxis grid;
run;

Thursday, January 12, 2012

Random seeds

A footnote toward Rick Wilkin’s recent article “How to Lie with a Simulation”. (Sit in front of a laptop w/o SAS; have to port all SAS/IML codes into R)

Generated 10 seeds randomly to run Stochastic simulation of Buffon's needle experiment by Rick's method. Hardly converge any of them ….

# Replicate Rick Wicklin's SAS/IML codes for Buffon's needle experiment
simupi <- function(N, seed) {
set.seed(seed)
z <- matrix(runif(N*2, 0, 1), N, 2)
theta <- pi*z[, 1]
y <- z[, 2] / 2
P <- sum(y < sin(theta)/2) / N
piEst <- 2/P
Trials <- 1:N
Hits <- (y < sin(theta)/2)
Pr <- cumsum(Hits)/Trials
Est <- 2/Pr
cbind(Est, Trials, seed)
}

# Generated 10 seeds randomly
seed_vector <- floor(runif(1:10)*10000)

# Each simulation with 50000 iterations
N <- 50000

# Run these 10 simulations
rt <- list()
for (i in 1:length(seed_vector)) {
rt[[i]] <- simupi(N, seed_vector[i])
}
results <- as.data.frame(do.call("rbind", rt))
results$seed <- as.factor(results$seed)

# Visualize individual simulation results
library(ggplot2)
p <- qplot(x = Trials, y = Est, data = results, geom = "line",
color = seed, ylim = c(2.9, 3.5))
p + geom_line(aes(x = Trials, y = pi), color = "Black")
ggsave("c:/plot1.png")



Averaging all results of the 10 simulations out. Then the curve converges easily. The application of this Monte Carlo simulation in Buffon's needle experiment is explained here by Rick Wicklin.

# Visuazlie the average result
rtmean <- aggregate(Est ~ Trials, data = results, mean)
p <- qplot(x = Trials, y = Est, data = rtmean, geom = "line")
p + geom_line(aes(x = Trials, y = pi), color = "Red")
ggsave("c:/plot2.png")

Tuesday, January 10, 2012

Benchmarking of the CUSUM function in SAS/IML


Cumulative sums can be obtained in SAS’s DATA step by the RETAIN statement. As the codes below, a new variable of the cumulative values will be returned by DATA step’s implicit DO loop.

data one;
do i = 1 to 1e6;
z = ranuni(0);
output;
end;
drop i;
run;

data two;
set one;
retain y 0;
y + z;
run;

The same logic can be realized by an explicit DO loop in SAS/IML. However, SAS/IML has a function CUSUM which is specially made for cumulative sums. To compare the efficiency of the two methods, I tried an experiment: calculating the cumulative sums for a random vector with incremental sizes from 1 million to 20 million by 1 million. The result shows that the CUSUM function is always 100 times faster than a raw DO loop. With the increase of the vector size, the gap gets wider. When the number of elements in the vector reaches 20 million, a Do loop would ask almost 40 seconds, while the CUSUM function requires only 0.35 seconds. Therefore, the vectorized CUSUM function is going to be a good candidate to do such jobs in SAS/IML.

proc iml;
a = t(do(1e6, 2e7, 1e6));
timer = j(nrow(a), 2);
do p = 1 to nrow(a);
n = a[p];
z = t(ranuni(1:n));

t0 = time();
x = cusum(z);
timer[p, 1] = time() - t0;

y = j(n, 1, .);
t0 = time();
do i = 1 to n;
if i = 1 then y[i] = z[i];
else y[i] = y[i-1] + z[i];
end;
timer[p, 2] = time() - t0;
end;

t = a||timer;
create _1 from t;
append from t;
close _1;
quit;

data _2;
set _1;
length test $100.;
label col1 = "Number of observations"
time = "Time by seconds for cumulative summation";
test = "The CUSUM function"; time = col2; output;
test = "DO loop"; time = col3; output;
keep test time col1;
run;

proc sgplot data = _2;
series x = col1 y = time / curvelabel group = test;
yaxis grid;
run;
P.S.
Thanks to Rick’s correction, I included his improved DO loop. Comparing with my previous DO loop with a conditional statement, this new DO loop is obviously more efficient. The conclusion is still the same: the CUSUM function outperforms the DO loops in this competition. And with more data involved, the contrast gets even more striking.


proc iml;
a = t(do(1e6, 2e7, 1e6));
timer = j(nrow(a), 3);
do p = 1 to nrow(a);
n = a[p];
z = t(ranuni(1:n));
/* 1 -- The CUSUM function */
t0 = time();
x = cusum(z);
timer[p, 1] = time() - t0;
/* 2 -- Bad DO loop*/
y = j(n, 1, .);
t0 = time();
do i = 1 to n;
if i = 1 then y[i] = z[i];
else y[i] = y[i-1] + z[i];
end;
timer[p, 2] = time() - t0;
/* 3 -- Good DO loop */
w = z;
t0 = time();
do i = 2 to n;
w[i] = w[i-1] + z[i];
end;
timer[p, 3] = time() - t0;
end;

t = a||timer;
create _1 from t;
append from t;
close _1;
quit;

data _2;
set _1;
length test $100.;
label col1 = "Number of observations"
time = "Time by seconds for cumulative summation";
test = "The CUSUM function"; time = col2; output;
test = "bad DO loop"; time = col3; output;
test = "good DO loop"; time = col4; output;
keep test time col1;
run;

proc sgplot data = _2;
series x = col1 y = time / curvelabel group = test;
yaxis grid;
run;

Saturday, January 7, 2012

DO loop vs. vectorization in SAS/IML


Vectorization is an important skill for many matrix languages. From Rick Wiklin’s book about SAS/IML and his recent cheat sheet, I found a few new vector-wise functions since SAS 9.22. To compare the computation efficiency between the traditional do loop style and the vectorization style, I designed a simple test in SAS/IML: square a number sequence(from 1 to 10000) and calculate the time used.

Two modules were written according to these two coding styles. Each module was repeated 100 times, and system time consumed was recorded by SAS/IML’s time() function.

proc iml;
start module1; * Build the first module;
result1 = j(10000, 1, 1); * Preallocate memory to the testing vector;
do i = 1 to 10000; * Use a do-loop to square the sequence;
result1[i] = i**2;
end;
store result1; * Return the resulting object;
finish;
t1 = j(100, 1, 1); * Run the first test;
do m = 1 to 100;
t0 = time(); * Set a timer;
call module1;
t1[m] = time() - t0;
end;
store t1;
quit;

proc iml;
start module2; * Build the second module;
result2 = t(1:10000)##2; * Vectorise the sequence;
store result2; * Return the resulting object;
finish;
t2 = j(100, 1, 1); * Run the second test;
do m = 1 to 100;
t0 = time(); * Set a timer;
call module2;
t2[m] = time() - t0;
end;
store t2;
quit;

proc iml;
load result1 result2; * Validate the results;
print result1 result2;
quit;

Then the results were released to Base SAS and visualized by a box plot with the SG procedures. In this experiment, the winner is the vectorizing method: vectorization seems much faster than do loop in SAS/IML. Therefore, my conclusions are: (1) avoid the do loop if possible; (2)use those vector-wise functions/operators in SAS/IML; (3) always test the speed of modules/functions by SAS/IML’s time() function.

proc iml;
load t1 t2;
t = t1||t2;
create _1 from t;
append from t;
close _1;
print t;
quit;

data _2;
set _1;
length test $25.;
test = "do_loop"; time = col1; output;
test = "vectorization"; time = col2; output;
keep test time;
run;

proc sgplot data = _2;
vbox time / category = test;
yaxis grid;
run;

Thursday, June 2, 2011

Using Proc IML for credit risk validation


Validation step is crucial for a scorecard in credit risk industry. Gunter and Peter mentioned in their fantastic book [Ref. 1] that cumulative accuracy profile (CAP) and receiver operating characteristic (ROC) are two popular methods. Thus, the values of accuracy ratio from CAP (or I refer it as Gini coefficient) and area under curve(AUC) from ROC would be important metrics to evaluate the discriminatory power of the scorecard. And actually they can be derived from each other by their linear relationship.

In the latest post of his blog, Rick Wicklin introduced how to implement the trapezoidal rule or calculate trapezoid areas under curve by a function of Proc IML [Ref. 2]. Although the same methodology can be realized by Data Step array or Proc FCMP, the beauty of Rick’s method is that it avoids the loops through IML’s matrix operation and therefore is more efficient and scalable. In the example below, I built Rick’s function into a macro to calculate AUC and accuracy ratio for a tiny testing dataset. The ‘TrapIntegral’ function can be further applied for validation of large-scale credit risk records.

References:
1. Gunter Löeffler and Peter Posch. ‘Credit Risk Modeling using Excel and VBA’. The 2nd edition. Wiley.
2. Rick Wicklin. ‘The Trapezoidal Rule of Integration’. The Do Loop. 01JUN2011.

/*******************READ ME*********************************************
* - Using Proc IML for credit risk validation -
*
* SAS VERSION: 9.2.2
* EXCEL VERSION: 2007
* DATE: 02jun2011
* AUTHOR: hchao8@gmail.com
*
****************END OF READ ME******************************************/

****************(1) MODULE-BUILDING STEP********************************;
%macro auc(data =, path =, filename =);
/*****************************************************************
* MACRO: auc()
* GOAL: calcuate auc and accuracy ratio for default risk
* PARAMETERS: data = dataset to use
* path = output path
* filename = name for validation card
*****************************************************************/
options mprint mlogic;
ods listing close;
proc sql;
select sum(default) into :totaldef
from &data;
quit;

data _tmp01;
set &data nobs = nobs;
xratio = ifn(default = 0, 1, 0)/(nobs - &totaldef);
yratio = default/&totaldef;
run;
proc sort data = _tmp01 out = _tmp02;
by descending rating descending default;
run;
data _tmp03;
set _tmp02;
by descending rating;
retain x y;
if _n_ = 1 then do; x = 0; y = 0; end;
x + xratio;
y + yratio;
if last.rating;
run;
data _tmp04;
if _n_ = 1 then do; x = 0; y = 0; end; output;
set _tmp03(keep = x y);
run;

proc iml;
use _tmp04;
read all var{x y};
start TrapIntegral(x,y);
N = nrow(x);
dx = x[2:N] - x[1:N-1];
meanY = (y[2:N] + y[1:N-1])/2;
return( dx` * meanY );
finish;
area = TrapIntegral(x,y);
acuratio = 2*area - 1;
call symput('area', left(char(area)));
call symput('acuratio', left(char(acuratio)));
quit;

ods html file = "&path&filename..xls" gpath = "&path" style = harvest;
title; footnote;
proc print data = &data label noobs;
run;

proc sgplot data = _tmp04;
series x = x y = y ;
scatter x = x y = y;
band x =x upper = y lower = 0 / transparency=.5;
xaxis grid;
yaxis grid;
inset "AUC is: %sysfunc(putn(&area, 8.4));
Accuracy Ratio is: %sysfunc(putn(&acuratio, 8.4))"
/ position = bottomright border;
keylegend "scatter";
run;
ods html close;

proc datasets;
delete _:;
quit;
ods listing;
%mend auc;

****************(2) TESTING STEP****************************************;
data test;
input Observation Rating $ Default;
label rating = 'Rating(A is best)'
default = 'Default(1=default)';
cards;
1 A 0
2 A 0
3 A 0
4 B 1
5 B 0
6 B 0
7 C 1
8 C 1
9 C 1
10 C 0
;;;
run;

%auc(data = test, path = h:\, filename = valid);
****************END OF ALL CODING***************************************;

Saturday, December 13, 2008

Proc GLM v.s. IML in ANOVA

data thrust;
input speed100 speed250 speed400 speed550;
if _n_ le 15 then feed=0.015;
if _n_ le 10 then feed=0.010;
if _n_ le 5 then feed=0.005;
cards;
121 98 83 58
124 108 81 59
104 87 88 60
124 94 90 66
110 91 86 56
329 291 281 265
331 265 278 265
324 295 275 269
338 288 276 260
332 297 287 251
640 569 551 487
600 575 552 481
612 565 570 487
620 573 546 500
623 588 569 497
;
run;

data thrust2;
set thrust;
array speedarray(4) speed:;
do i = 1 to 4;
force = speedarray(i);
output;
end;
drop speed:;
run;

data thrust3;
set thrust2;
if i = 1 then speed = 100;
if i = 2 then speed = 250;
if i = 3 then speed = 400;
if i = 4 then speed = 550;
drop i;
run;

proc sort data=thrust3 out=thrust4;
by feed speed;
run;

proc glm data=thrust3;
title 'GLM solution';
class feed speed;
model force= feed speed feed*speed;
run;

proc iml;
title 'IML solution';
use thrust4;
read all var{force} into y;
Xzero=J(60,1);
Xfeed=I(3)@j(4,1)@j(5,1);
Xspeed=j(3,1)@I(4)@j(5,1);
Xint=I(3)@I(4)@j(5,1);
Xfull=Xzero||Xfeed||Xspeed||Xint;

SSE=y`*y-y`*Xfull*ginv(Xfull`*Xfull)*Xfull`*y;
SSTrtFeed=y`*((I(3)-(1/3)*J(3,3))@((1/4)*J(4,4))@((1/5)*J(5,5)))*y;
SSTrtSpeed=y`*(((1/3)*J(3,3))@(I(4)-(1/4)*J(4,4))@((1/5)*J(5,5)))*y;
SSInt=y`*((I(3)-(1/3)*J(3,3))@(I(4)-(1/4)*J(4,4))@((1/5)*J(5,5)))*y;
SSTotal=y`*(I(60)-(1/60)*J(60,60))*y;

print SSTrtFeed SSTrtSpeed SSInt SSE SSTotal;

MSTrtFeed=SSTrtFeed/2; MSTrtSpeed=SSTrtSpeed/3; MSInt=SSInt/6; MSE=SSE/48;
print MSTrtFeed MSTrtSpeed MSInt MSE ;

F_TrtFeed=MSTrtFeed/MSE; F_TrtSpeed=MSTrtSpeed/MSE; F_Int=MSInt/MSE;
Print F_TrtFeed F_TrtSpeed F_Int;

P_TrtFeed=1-probf(F_TrtFeed,2,48);
P_TrtSpeed=1-probf(F_TrtSpeed,3,48);
P_Int=1-probf(F_Int,6,48);
Print P_TrtFeed P_TrtSpeed P_Int;
run;quit;