dm "log; clear; odsresults; clear;"; options linesize=64 nonumber nodate; *********************************************************************; * AUTHOR: Chris Bilder *; * DATE: 5-18-16 *; * PURPOSE: Basic matrix algebra *; *********************************************************************; title1 "Chris Bilder, STAT 850"; proc iml; *******************; * Basic operations ; A = {1 2 3, 4 5 6}; print A; B = {-1 10 -1, 5 5 8}; print B; print A B; print A, B; temp = A[1,1]; print temp; temp = A[1,]; print temp; temp = A[,1:2]; print temp; temp = A[,1]; print temp; C = A + B; print C; D = A - B; print D; *******************; * Matrix multiplication; A = {1 2 3, 4 5 6}; B = {3 0, 1 2, 0 1}; print A, B; E = A*B; print E; F = B*A; print F; A3 = A*3; print A3; *Elementwise multiplication; C = A#A; print C; *******************; * Inverse; A = {1 2, 3 4}; print A; Ainv = inv(A); print Ainv; check = A*Ainv; print check; *******************; * More matrix operations; A = {1 2, 3 4}; print A; Adiag = diag(A); print Adiag; *Pulls out diagonal of matrix and creates a new matrix (see if can improve); *Determinant; Adet = det(A); *Tranpose; Atran = t(A); *A` also works; print Adet, Atran; *Create diagonal matrix; * Create 1x5 vector and use elementwise multiplication with a 5x5 identity matrix; DiagMat = {1 2 3 4 5}#I(5); print DiagMat; *******************; * Eigenvalues; A = {1 0.5, 0.5 1.25}; print A; lambda = eigval(A); lambdavec = eigvec(A); print lambda, lambdavec; check1left = A*lambdavec[,1]; check1right = lambda[1]*lambdavec[,1]; print check1left, check1right; check2left = A*lambdavec[,2]; check2right = lambda[2]*lambdavec[,2]; print check2left, check2right; quit; **********************************************************; * Example of how to get data into proc iml from; * a SAS data set; *Read in HS and College GPA data set; data set1; infile "C:\data\gpa.csv" firstobs=2 delimiter=","; input HS College ; run; proc iml; use set1; read all var{HS} into Xpart; read all var{College} into Y; *Number of observations; n = nrow(Xpart); *nx1 vector of 1's; vec1 = J(n, 1, 1); *Horizontal concatenate; X = vec1 || Xpart; print X, Y; betahat = inv(t(X)*X) * t(X) * Y; print betahat; *Shows how to read in all variables at same time; read all into XY; print XY; *Show how to read in data under particular conditions; * I put the { } part as an example of how to read in 2 or more variables at a time; * in a general manner (key point: no comma needed to separate variable names; read all var {HS College} where(HS > 3.0 & College > 3.0) into BorHigher; print BorHigher; quit; ******************************************************************; *Show how to use R within SAS; proc iml; *Simple example; submit / R; x <- 1 x endsubmit; *Use a SAS data set in R; run ExportDataSetToR("set1", "Rset1"); submit / R; head(Rset1) mod.fit <- lm(formula = College ~ HS, data = Rset1) names(mod.fit) mod.fit$coefficients #summary(mod.fit) #Does not work! plot(x = Rset1$HS, y = Rset1$College, xlab = "HS GPA", ylab = "College GPA", main = "College GPA vs. HS GPA", xlim = c(0,4.5), ylim = c(0,4.5), col = "red", pch = 1, cex = 1.0, lwd = 2.0, panel.first = grid()) curve(expr = predict(object = mod.fit, newdata = data.frame(HS = x)), col= "blue", add = TRUE, lwd = 2, xlim = c(min(Rset1$HS), max(Rset1$HS))) curve(expr = predict(object = mod.fit, newdata = data.frame(HS = x), interval = "confidence", level = 0.95)[,2], col= "blue", add = TRUE, lwd = 2, xlim = c(min(Rset1$HS), max(Rset1$HS)), lty = "dashed") curve(expr = predict(object = mod.fit, newdata = data.frame(HS = x), interval = "confidence", level = 0.95)[,3], col= "blue", add = TRUE, lwd = 2, xlim = c(min(Rset1$HS), max(Rset1$HS)), lty = "dashed") endsubmit; *Use R vector in proc iml; run ImportMatrixFromR(betahats, "mod.fit$coefficients"); * call ImportMatrixFromR(betahats, "mod.fit$coefficients"); *Also works; print betahats; *Use R data frame outside of proc iml; run ImportDataSetFromR("set2", "Rset1"); *Use a proc iml matrix in R; use set1; read all into XY; print XY; run ExportMatrixToR(XY, "Rset2"); submit / R; head(Rset2) endsubmit; quit; proc print data = set2(obs=5);; run; *Check of R versions and other information about R from SAS; proc iml; submit / R; print("System information") Sys.info() print("\n R version") #\n does not work properly here R.Version() #Looks like SAS used my most current version of R print("\n Session information") sessionInfo() endsubmit; quit; *Verifies in log window that R is available; proc options option=RLANG; run; **************************************************; *Other things you can do in IML:; * 1) Plots; * 2) Call procs; * 3) Open files directly