*************************************;
*		STATISTICAL COMPUTING		*;
*		SAS: USING PROC IML			*;
*			SPRING 2014				*;
*************************************;


/************************************/
/*		GENERATING MATRICES			*/
/************************************/

PROC IML;
A	= 1; 		/* CREATE A SCALAR*/
B	= {1 2 3}; /* CREATE A ROW VECTOR OF LENGTH 3*/
C	= {	4,
		5,
		6};		/* CREATE A COLUMN VECTOR OF LENGTH 3*/
D	={	1 2,
		3 4,
		5 .};	/* CREATE A 3 BY 2 MATRIX WHERE THE 3,2 ELEMENT IS MISSING*/
PRINT A B C D;	/* DISPLAY THE  MATRICES IN THE OUTPUT*/
QUIT;

/************************************/
/*		MANIPULATE MATRICES			*/
/************************************/

PROC IML;
	REPEAT_O1={[2]"YES" [2] "NO"}; /*USING THE REPEAT FUNCTION TO FILL THE MATRIX*/
	REPEAT_O2={"YES" "YES" "NO" "NO"}; /* REPEATING ELEMENTS MANUALLY*/
		PRINT REPEAT_O1 REPEAT_O2;
/*******************************************/
	A={1 2,
	   3 4};	/* DEFINE MATRIX*/
	A1=A[2,1];	/* SELECT THE ELEMENT IN THE 2ND ROW, FIRST COLUMN: A1 SOULD EQUAL 3 */
	A2=A[1,];	/* SELECT THE FIRST ROW, A2 SHOULD EQUAL A 2 X 1 VECTOR {1 2} */
	A3=A[,1];	/* SELECT THE FIRST COLUMN, A3 SHOULD EQUAL A 1 X 2 VECTOR {1,3} */
	B={1 2 0 0, 3 4 0 0}; /* DEFINE A MATRIX B, WITH TWO SUBMATRICES A AND A 2 X 2 NULL MATRIX*/
		A_NEW=B[1:2,1:2]; /* RECOVER THE A MATRIX FROM B */
		A_NEW2=B[,{1 2}]; /*RECOVER THE A MATRIX FROM B, ANOTHER WAY TO WRITE IT*/
/*******************************************/
	C_ROWNM={M F}; /* SET ROW NAMES FOR MATRIX C*/
	C_COLNM={TRUE FALSE}; /* SET COL NAMES FOR MATRIX C*/
	C={10 25,9 18};
	PRINT A A1 A2 A3 B A_NEW 
			C[ROWNAME=C_ROWNM COLNAME=C_COLNM FORMAT=6.1 LABEL="MY MATRIX"] /*MODIFYING PRINTED OUTPUT FOR MATRIX C*/;
/*******************************************/
	C_NEW=C; /* CREATING A DUPLICATE MATRIX*/
		MATTRIB C_NEW ROWNAME=C_ROWNM COLNAME=C_COLNM FORMAT=6.1 LABEL="MY MATRIX"; /* PERMANANTLY CHANGING OUTPUT FORMAT*/
		PRINT C C_NEW; /* COMPARING DIFFERENT APPROACHES*/
/*******************************************/
	D=A[LOC(A>1)];/* SELECTING ONLY ELEMENTS THAT MEET RULE, NOTE THAT MATRIX STRUCTURE NOT RETAINED*/
		PRINT A D;
/*******************************************/
	E_TEMP=A; /* CREATING A DUPLICATE MATRIX*/
		E_TEMP[1,1]=25 /* CHANGING A SINGLE ELEMENT*/
			PRINT E_TEMP;
		E_TEMP[LOC(E_TEMP>5)]=.; /* SETTING ALL ELEMENTS MEETING RULE TO MISSING*/
			PRINT E_TEMP;
QUIT;

/************************************/
/*		MATRIX ALGEBRA				*/
/************************************/
PROC IML;
	*MATRIX ADDITION;
	A={1 3, 2 5}; /*DEFINE MATRIX*/
	B={-5 2, 7 0}; /*DEFINE MATRIX*/
	C=A+B; /* ADD A AND B*/
	PRINT C;

	*MATRIX MULTIPLICATION;
	A={2 3,4 5}; /*DEFINE MATRIX*/
	B={1 6,2 0}; /*DEFINE MATRIX*/
	AB=A*B; /*MULTIPLY A BY B*/
	BA=B*A; /* MULTIPLY B BY A*/
	PRINT A B AB BA; /* NOTE THAT MULTIPLICATION IS NOT COMMUTATIVE, AB DOESN'T EQUAL BA*/
QUIT;

/************************************/
/*		SYSTEM OF EQUATIONS			*/
/************************************/
PROC IML;
	A={3 2 -4,
		5 -4 0,
		0 3 10};
	B={11,9,42};
OPT1=SOLVE(A,B);
OPT2=INV(A)*B;
PRINT OPT1 OPT2;
QUIT;

/************************************/
/*		READING SAS DATASETS		*/
/************************************/

*CREATING A SAS DATASET TO WORK WITH;
DATA MYDATA;
SET SASHELP.CARS;
RUN;

PROC IML;
	USE MYDATA VAR {MSRP MPG_CITY MPG_HIGHWAY} ; /* OPEN DATASET*/
	READ ALL VAR _ALL_ WHERE (MSRP<12000) INTO CAR_MAT; /* READ DATASET*/
		Z=NROW(CAR_MAT); /* FIGURE OUT HOW MANY ROWS*/
	PRINT Z CAR_MAT[COLNAME={MSRP CITY HWY}]; /* LOOK AT DATA*/

	MEAN_LEV=CAR_MAT[:,]; /*TAKE THE MEAN LEVEL FOR EACH VARIABLE*/
	PRINT MEAN_LEV;
QUIT;

/************************************/
/*		CONDITIONAL PROCESS			*/
/************************************/
PROC IML;
	A={12 22 33}; /* CREATE MATRIX*/
IF MAX(A)<20 /* CONDITIONAL ARGUMENT*/
	THEN P=1; /* ACTION IF TRUE*/
	ELSE P=0; /* ACTION IF FALSE*/
PRINT P;
QUIT;

PROC IML;
	Y=0;
	DO I=1 TO 3;
		Y=Y+1;
		PRINT Y;
	END;
QUIT;
PROC IML;
COUNT=1;
	DO WHILE(COUNT<3);
		COUNT=COUNT+1;
		PRINT "SUCCESS";
	END;

COUNT=1;
	DO UNTIL(COUNT>3);
		COUNT=COUNT+1;
		PRINT "SECOND ATTEMPT";
	END;
QUIT;

/************************************/
/*		PROC						*/
/************************************/
PROC IML;
Q={2 5 7 9};
	CREATE MYDATA VAR{Q};
	APPEND;
	CLOSE MYDATA;
*Table=“Moments”;
SUBMIT;
	*SUBMIT table;
	*ODS SELECT &table;
PROC UNIVARIATE DATA=MYDATA;
		VAR Q;
		ODS OUTPUT MOMENTS=MOMENTS;
RUN;
ENDSUBMIT;
USE MOMENTS;
READ ALL VAR{NVALUE1 LABEL1};
CLOSE MOMENTS;
LABL ="MY OUTPUT";
PRINT NVALUE1[ROWNAME=LABEL1 LABEL=LABL];
QUIT;

PROC OPTIONS OPTION=RLANG;run;

proc iml;
/* Comparison of matrix operations in IML and R */
print "---------- SAS/IML Results -----------------";
x = 1:3; /* vector of sequence 1,2,3 */
m = {1 2 3, 4 5 6, 7 8 9}; /* 3 x 3 matrix */
q = m * t(x); /* matrix multiplication */
print q;

print "------------- R Results --------------------";
submit / R;
rx <- matrix( 1:3, nrow=1) # vector of sequence 1,2,3
rm <- matrix( 1:9, nrow=3, byrow=TRUE) # 3 x 3 matrix
rq <- rm %*% t(rx) # matrix multiplication
print(rq);
endsubmit;
submit/R;
plot(rq,main="main")
endsubmit;
quit;

proc iml;
use Sashelp.Class;
read all var {Weight Height};
close Sashelp.Class;
/* send matrices to R */
call ExportMatrixToR(Weight, "w");
call ExportMatrixToR(Height, "h");
submit / R;
Model <- lm(w ~ h, na.action="na.exclude") # a
ParamEst <- coef(Model) # b
Pred <- fitted(Model)
Resid <- residuals(Model)
endsubmit;
call ImportMatrixFromR(pe, "ParamEst");
print pe[r={"Intercept" "Height"}];
ht = T( do(55, 70, 5) );
A = j(nrow(ht),1,1) || ht;
pred_wt = A * pe;
print ht pred_wt;

YVar = "Weight";
XVar = "Height";
submit XVar YVar / R;
Model <- lm(&YVar ~ &XVar, data=Class, na.action="na.exclude")
print (Model$call)
endsubmit;

quit;