In a previous post, I described ways to create SAS/IML vectors that contain uniformly spaced values. The methods did not involve writing any loops.

This post describes how to perform a similar operation: creating evenly spaced values on a two-dimensional grid. The DATA step solution is simple, but an efficient SAS/IML solution requires thinking about the problem in terms of vectors and matrices.

### The DATA Step Solution

A grid of values is useful in many contexts, such as scoring and visualizing a statistical model. In the DATA step, you can use loops to create a grid of values. For example, to generate a 3 x 2 grid of values on the rectangle [0,2] x [0,1] you can use the following program:

data grid; keep x y; /** set parameters for the general problem **/ XMin = 0; XMax = 2; Nx = 3; YMin = 0; YMax = 1; Ny = 2; dx = (XMax-XMin)/(Nx-1); dy = (YMax-YMin)/(Ny-1); do y = YMin to YMax by dy; do x = XMin to XMax by dx; output; end; end; run; proc print; var x y; run; |

This solution uses a double loop to compute the values on a uniform grid. There are `Nx` points in the X direction, evenly spaced on the interval [XMin, XMax], and
`Ny` points in the Y direction, evenly spaced on the interval [YMin, YMax].

As I've written before, the SAS/IML language is a vector-matrix language, and an efficient SAS/IML program is one that avoids writing unnecessary loops. How can you write a similar program in PROC IML that avoids one or both loops?

### A SAS/IML Language Solution

**UPDATE 2014: In SAS/IML 12.3** (released with SAS 9.4), you can use the EXPANDGRID function to generate a grid in multiple dimensions! The following program is **not needed** if you have a recent version of SAS.

There are three simple SAS/IML functions that you can use to create a module that computes a grid of evenly spaced values:

- The DO function for generating sequential values, which is described in a previous post.
- The REPEAT function, which repeats a pattern a specified number of times.
- The SHAPE function, which reshapes a vector into a matrix. I also wrote a previous post on reshaping matrices.

The first step is to use the parameters in the problem to determine the points in the X and Y directions that will form the grid:

proc iml; XMin = 0; XMax = 2; Nx = 3; YMin = 0; YMax = 1; Ny = 2; XDiv = do(XMin, XMax, (XMax-XMin)/(Nx-1)); YDiv = do(YMin, YMax, (YMax-YMin)/(Ny-1)); print XDiv, YDiv; |

These are the coordinates of the grid points, but you need to arrange them into a matrix of ordered pairs that contains
`Nx` * `Ny` rows and 2 columns.

The X coordinates are easy to construct. To form the X coordinates you can repeat the XDiv vector `Ny` times, and then reshape those values into a column vector:

X = repeat(XDiv, Ny); X = shape(X, Nx*Ny); /** or use COLVEC module in IMLMLIB **/ print X; |

The Y coordinates are a little more difficult to construct. Recall that SAS/IML matrices are stored row-wise. This means that if you use the REPEAT function to form a set of points from the `YDiv` vector, the points will be in the wrong order:

Y = repeat(YDiv, Nx); print Y; /** wrong order: elements stored row-wise **/ |

What you really want is to traverse the `Y` matrix down the columns instead of across the rows. No problem! Use the T function to transpose the values, and then reshape those transposed values into a column:

Y = shape(T(Y), Nx*Ny); /** or use COLVEC module in IMLMLIB **/ |

### A SAS/IML Module That Defines Ordered Pairs on a Grid

You can put these statements into a module, which enables you to conveniently generate a grid of points:

/** Module to define ordered pairs in a uniform grid of points. Input: XMin = minimum value of the X coordinate of the grid XMax = maximum value of the X coordinate of the grid Nx = number of grid lines along the X coordinate YMin = minimum value of the Y coordinate of the grid YMax = maximum value of the Y coordinate of the grid Ny = number of grid lines along the Y coordinate Output: Ordered pairs are returned in a matrix with (Nx x Ny) rows and two columns. **/ start Define2DGrid(XMin, XMax, Nx, YMin, YMax, Ny); XDiv = do(XMin, XMax, (XMax-XMin)/(Nx-1)); YDiv = do(YMin, YMax, (YMax-YMin)/(Ny-1)); X = repeat(XDiv, Ny); X = shape(X, Nx*Ny); Y = repeat(YDiv, Nx); Y = shape(T(Y), Nx*Ny); return ( X || Y ); finish; /** Example: Use the module to define a grid of points and then evaluate a function f(x,y) on the grid. **/ p = Define2DGrid(-2, 2, 23, /** 23 points in [-2,2] **/ -1, 1, 11); /** 11 points in [-1,1] **/ x = p[,1]; y = p[,2]; z = x##2 - 2*y##2; /** evaluate a function z = f(x,y) **/ /** the following SAS/IML Studio statement creates a rotating plot **/ RotatingPlot.Create("Quadratic Surface", x, y, z); |

## 3 Comments

Pingback: Random uniform versus uniformly spaced: Applying statistics to show choir - The DO Loop

Pingback: Compute the multivariate normal denstity in SAS - The DO Loop

Pingback: Using Newton's method to find the zero of a function in SAS