Mathematical art (part 2): Unweaving matrices


In my previous blog post, I showed how you can use SAS to program a "weaving" algorithm that takes an image, cuts it into strips, and weaves the strips together to create mathematical art. I used matrices and heat maps for the computations and visualization.

At the end of the post, I presented an image of woven strips (shown at left) and challenged the reader to "deduce the original image from the weave." This article shows that you can "unweave" the image, reassemble the strips, and use interpolation to obtain an approximation to the original image. You can download the SAS/IML program that computes and visualizes the preimage.

Unweaving a weave

The first step in reconstructing the original image is to "pull apart" the horizontal and vertical strips, rotate the horizontal ones, put them next to the vertical strips, and close all the gaps. This requires writing the inverse of the weave algorithm. When you do this, you obtain the following image:


Notice that there are gaps in our knowledge about the original image. Due to the weaving process, there are places where a vertical strip was covered by a horizontal strip, or the other way around. These missing values are shown in white. Mathematically, the weaving algorithm represents a many-to-one function, so there is no inverse function that can restore the preimage from the image.

Nevertheless, the unweaving reveals the main features of the original image. And we can do more. If we know that the original image was created from a continuous mathematical function, you can use interpolation to try to fill in the gaps. For this image you can use the information from nearby cells (north, south, east, and west) to estimate a value for the missing cells. The following uses a simple estimation scheme, which is to assign each missing cell the average value of its neighbors. You can proceed in five steps:

  1. Approximate the corners. For the lower left corner, use the average of the cell above and the cell to the right. For the lower right corner, use the average of the cell above and the cell to the left.
  2. Approximate the left edge. For each missing cell, average the cells above, below, and to the right.
  3. Approximate the right edge. For each missing cell, average the cells above, below, and to the left.
  4. Approximate the bottom edge. For each missing cell, average the cells above, to the left, and to the right.
  5. Approximate the interior cells. For each missing cell, average the cells above, below, to the left, and to the right.

After implementing this process, you obtain the following image, which approximates the original image:


How well does this image approximate the original image? Very well! The mean square error is very small, and the maximum error values are within 20% of the true values. For comparison, here is the original image. You can see that it is very similar to the reconstructed image.


The approximation could be improved by using a more sophisticated interpolation scheme, such as bivariate linear interpolation.

Feel free to download the SAS/IML program and play with other weaves. How well can you recover an original image from its weave?
Post a Comment

Mathematical art: Weaving matrices


An artist friend of mine recently created a beautiful abstract image and described the process on her blog. She says that "after painting my initial square, I cut it into strips and split them down the middle, then wove them together.... I had no idea when I started piecing these together again what they would look like."

I marveled at the beauty of the finished result, but I also started thinking about the process of constructing the artwork. You cut an initial picture into n vertical strips. (For simplicity, assume that n is even.) Then you take the first n/2 strips and space them out so that there are gaps between them. You take the last n/2 strips and turn them sideways, weaving them in and out of the vertical strips.

As my friend said in a private communication, the final image "depends entirely on the arrangement of the elements.... There are so many versions of the final outcome." By this she means that the process is deterministic, although the artist can make choices regarding the number of strips, the orientation of the strips, and the weaving pattern. Yes, it is difficult to predict the final image while staring at the initial image, but this is the sort of algorithm that you could ask a computer to carry out.

So, of course, I had to learn to make images like this with SAS!

Weaving matrices

To make sure that I understood the algorithm, I started with a 10 x 10 numeric matrix X that contains the values 1–100 in row-major order. That is, the first row contains 1–10, the second row contains 11–20, and so forth. I then created a 10 x 10 matrix W (for "weave") that contained all missing values. I copied the first five columns of X into the even columns of W. I then copied the last five columns of X into the even rows of W. I then went back and "wove" the columns and rows by re-copying certain elements from the first five columns of X. You can view this construction process for an example matrix.

The process can be encapsulated into a SAS/IML function, which is shown below:

proc iml;
start weave(x);
   n = nrow(x);             /* n must be even */
   y = j(n, n, .);
   j = 1:n/2;               /* vertical strips */
   y[, 2*j] = x[, j];       /* copy to every other column */
   i = (n/2+1):n;           /* horizontal strips */
   y[2*i-n,] = x[, i]`;     /* copy to every other row */
   k = do(1, n/2, 2);       /* weaving: raise certain rows */
   y[2*k, 2*k] = x[2*k, k]; /* every other for even rows */
   k = do(2, n/2, 2);
   y[2*k, 2*k] = x[2*k,k];  /* every other for odd rows */
   return( y );
x = shape(1:100, 10);       /* original 10x10 matrix */
w = weave(x);
print w;

Weaving heat maps

Now the fun can begin. Recall that a heat map is essentially a visualization of a matrix of numbers. A heat map can be used to visualize mathematical functions of two variables, and you can use the HEATMAPCONT function in SAS/IML 13.1 to easily create a heat map.

You can use the EXPANDGRID function to define a grid of evenly spaced points in the square [-1, 1] x [-1, 1]. You can evaluate a mathematical function on the grid and reshape the results into a matrix. You can then create a heat map of the original function and the "weave" of the function. For details, you can download the complete SAS/IML program that I used to construct the images in this article.

For example, in the following image, the function on the left is f(x,y) = x. The heat map that visualizes the function shows that the function is constant on vertical strips. You can see that each strip is a constant color, and the image on the right has taken those strips and woven them together.


The woven image is beautiful, but you can make more interesting images by choosing more interesting functions to use as the original image. The left portion of the following image shows a heat map for the quadratic function f(x,y) = x2 + y2. The image on the right is formed by weaving the vertical strips from the left image.


You get the idea: You can visualize any mathematical function and then form "mathematical art" by cutting the function along strips and weaving the strips back together.

If you are adventurous, you can play with different color palettes for the images. I chose a muted rainbow palette, but you can experiment with other color ramps. The program for these examples include other color ramps that change the aesthetics of the images.


I'll leave you with one more image. It is the weave for a common function from multivariable calculus. Can you deduce the original image from the weave? If not, look in the program that creates these examples..

These images show flexibility of the SAS/IML language for creating a custom algorithm. These simple images lack the depth and complexity of "real" art, but they show how matrices can be used to illustrate the main ideas behind a simple artistic process.

Try it yourself! What images can you make? If you create an interesting picture, leave a comment that includes the mathematical formula.

Post a Comment

Compute the number of digits in an integer

The title of this blog post might seem strange, but I occasionally need to compute the number of digits in a number, usually because I am trying to stuff an integer value into a string. Each time, I have to derive the formula from scratch, so I am writing this article so that I and others can look up the formula next time we need it.

The problem: You have a positive integer, n, which has k digits. How can you efficiently find k as a function of n?

The solution is to use logarithms. Because n has k digits,
10k-1n < 10k, or
10k-1 < n+1 ≤ 10k.

Applying the base-10 logarithm produces
k-1 < log10(n+1) ≤ k.

The CEIL function in Base SAS computes the integer value greater or equal to its argument, which proves that
k = ceil( log10(n+1) ).

A canonical example, is n=1234, which has k=4 digits. How can you find k from n? The following program computes k:

data _null_;
n = 1234;      /* positive integer */
k = ceil(log10(n+1)); 
put "The integer " n= "has " k= "digits.";
The integer n=1234 has k=4 digits.

If you have a negative integer m ≤ -1 and the absolute value has k digits, then m requires a string with k+1 elements, because one element is needed for the negative sign.

Incidentally, there is nothing special about base 10. The proof generalizes to compute the number of digits required to represent a number in any base. For example, the number n=1234 written in base 2 requires k = ceil(log2(n+1)) = 11 digits. You can check that the 11-digit binary number 10011010010 equals the decimal value 1234.

You can also use these mathematical ideas to compute the next power of 2 (or power of 10) that is larger than a given number.

Post a Comment

Monitor convergence during simulation studies in SAS

Ugh! Your favorite regression procedure just printed a warning to the SAS log. Something is wrong, and your attempt to fit a model to the data has not succeeded. A typical message is "WARNING: The validity of the model fit is questionable," perhaps followed by some additional diagnostic messages about "quasi-separation" or "lack of convergence."

If your modeling toolkit includes procedures such as LOGISTIC, GENMOD, MIXED, NLIN, or PHREG, you might have experienced convergence problems. A small sample size or a misspecified model are among the reasons for lack of convergence. There are many papers that discuss how to handle convergence issues. Paul Allison (2008) wrote a paper on some reasons that a logistic model might fail to converge, including an explanation of quasi-complete separation. The documentation for the MIXED procedure includes a long list of potential reasons that a mixed model might fail to converge. Several of my friends from SAS Technical Support give advice on convergence in their SAS Global Forum paper (Kiernan, Tao, and Gibbs, 2012).

Although it can be frustrating to deal with convergence issues during a data analysis, lack-of-convergence during a simulation study can be maddening. Recall that the efficient way to implement a simulation study in SAS is to use BY-group processing to analyze thousands of independent samples. A simulation study is designed to run in batch mode without any human intervention. How, then, can the programmer deal with lack of convergence during a simulation?

Some SAS users (such as Chen and Dong, 2009) have suggested parsing the SAS log for notes and warning messages, but that approach is cumbersome. Furthermore, if you have turned off SAS notes during the simulation, then there are no notes in the log to parse!

The _STATUS_ variable in the OUTEST= data set

Fortunately, SAS procedures that perform nonlinear optimization provide diagnostic variables as part of their output. These variables are informally known as "status variables." You can monitor the status variables to determine the BY groups for which the optimization converged.

The easiest way to generate a status variable is to use the OUTEST= option to generate an output data set that contains parameter estimates. Not every procedure supports the OUTEST= option, but many do. Let's see how it works. In my article "Simulate many samples from a logistic regression model," I showed how to generate 100 samples of data that follow a logistic regression model. If you make the sample size very small (like N=20), PROC LOGISTIC will report convergence problems for some of the random samples, as follows:

%let N = 20;                                     /* N = sample size */
%let NumSamples = 100;                           /* number of samples */
/* Generate logistic data: See */
options nonotes;    /* turn of notes; use OUTEST= option */
proc logistic data=LogisticData noprint outest=PE;
   by SampleId;
   model y(Event='1') = x1 x2;
options notes;
proc freq data=PE;
tables _STATUS_ / nocum;

The output from PROC FREQ shows that 10% of the models did not converge. The OUTEST= data set has a variable named _STATUS_ that indicates whether the logistic regression algorithm converged. You can use the _STATUS_ variable to analyze only those parameter estimates from converged optimizations. For example, the following call to PROC MEANS produces summary statistics for the three parameter estimates in the model, but only for the converged models:

proc means data=PE nolabels;
where _STATUS_ = "0 Converged"     /* analyze estimates for converged models */
var Intercept x1 x2;

The ConvergenceStatus table

Notice that the _STATUS_ variable does not contain details about why the algorithm failed to converge. For that information, you can use the ConvergenceStatus ODS table, which is produced by all SAS/STAT regression procedures that involve nonlinear optimization. To save the ODS table to an output data set, you cannot use the NOPRINT option. Instead, use the %ODSOff and %ODSOn macros to suppress ODS output. Use the ODS OUTPUT statement to write the ConvergenceStatus tables to a SAS data set, as follows:

/* define %ODSOff and %ODSOn macros */
proc logistic data=LogisticData;
   by SampleId;
   model y(Event='1') = x1 x2;
   ods output ParameterEstimates=PE2 ConvergenceStatus=CS;
proc freq data=CS;
tables Status Reason / nocum;

The ConvergenceStatus table has a numerical variable named Status, which has the value 0 if the algorithm converged. The character variable Reason contains a short description of the reason that the algorithm terminated. The results of the PROC FREQ call shows that the logistic algorithm failed eight times because of "complete separation," converged 90 times, and failed twice due to "quasi-complete separation."

You can use the DATA step and the MERGE statement to merge the ConvergenceStatus information with the results of other ODS tables from the same analysis. For example, the following statements merge the convergence information with the ParameterEstimates table. PROC MEANS can be used to analyze the parameter estimates for the models that converged. The summary statistics are identical to the previous results and are not shown.

data All;
merge PE2 CS;
by SampleID;
proc means data=All;
where Status = 0;
class Variable;
var Estimate;

In conclusion, the ConvergenceStatus table provides information about the convergence for each BY group analysis. When used as part of a simulation study, you can use the ConvergenceStatus table to manage the analysis of the simulated data. You can count the number of samples that did not converge, you can tabulate the reasons for nonconvergence, and you can exclude nonconverged estimates from your Monte Carlo summaries.

Post a Comment

Video: Ten tips for simulating data with SAS

One of my presentations at SAS Global Forum 2015 was titled "Ten Tips for Simulating Data with SAS". The paper was published in the conference proceedings several months ago, but I recently recorded a short video that gives an overview of the 10 tips:

If your browser does not support embedded video, you can go directly to the video on YouTube.

The tips are based on material in my book Simulating Data with SAS, which contains hundreds of tips and techniques for simulating data, and several tips have been published as articles in my blog.

Post a Comment

She wants to be an airborne ranger

I wanna be an airborne ranger,
Live the life of guts and danger.*

If you are an 80's movie buff, you might remember the scene in The Breakfast Club where Bender, the juvenile delinquent played by Judd Nelson, distracts the principal by running through the school singing this song. Recently, the US Army Ranger School has been in the news because, for the first time, two women passed the rigorous training regime.

It is extremely difficult for anyone, regardless of gender, to pass the course. Only 94 out of 381 males passed this year's course, and only two of the original 19 females graduated.

Those numbers imply that that approximately 20% of males graduated, whereas only 10% of females graduated. However, the initial female group was so small that I wondered whether this difference was statistically significant. Do these data indicate that there is a higher probability of passing the ranger course if you are male?

It turns out that the difference is not statistically significant. You can use PROC FREQ in Base SAS to analyze the data, as follows:

data Rangers;
input Gender $ Finished $ Count;
M  Yes  94
M  No  381
F  Yes   2
F  No   19
proc Freq data=Rangers;
   weight Count;
   tables Gender * Finished / nocol nopercent 
          riskdiff(cl=score column=2 norisks);

The first row of the crosstabulation shows the counts and percentages of women who failed or passed the ranger course. The second row presents the counts and percentages for men. The first column shows that about 90% of the women failed the course, as compared with 80% of the men.

The options on the TABLES statement request several statistical tests. The CHISQ option requests a chi-square test for association between gender and finishing the course. Because one cell has fewer than 5 counts, it is best to look at the exact Fisher test.


The p-value is not small, so there is insufficient evidence to reject the hypothesis that finishing the course is independent of gender. The other statistical tests (not shown) show similar results. These data do not suggest that gender significantly affects the odds of failing or the risk of failing. The confidence interval for the odds ratio includes 1; the confidence interval for the risk difference includes 0.

To obtain statistical significance, a larger cadre of women would need to be in the study. With the US Army re-evaluating the role of women in combat positions, maybe that will occur in the near future.

For now, congratulations to Capt. Kristen Griest and 1st Lt. Shaye Haver, the hard-working women who made history by passing the ranger course, and to all the women who attempted. The data indicate that the pass rate between men and women is not statistically significant. I think that is quite an achievement.

* The second line of this cadence has many variations, including some that are not safe for work.

Post a Comment

Correlations between groups of variables

Typically a correlation analysis reports the correlations between all pairs of variables, including the variables with themselves. The resulting correlation matrix is square, symmetric, and has 1s on the main diagonal. But suppose you are interested in only specific combinations of variables. Perhaps you want the pairwise correlations between one group that contains n1 variables and another group of n2 variables.

No worries. The CORR procedure in Base SAS supports the WITH statement which enables you to compute and display a subset of the correlations between variables.

Correlations of certain variables WITH other variables

Often variables naturally fall into groups. For example, in the Sashelp.Cars data, the variables are related to four groups: price, power of the engine, fuel economy, and size of vehicle. In a medical study you might have certain variables related to physical characteristics that will not be affected by the study (age, height, gender) and other variables that might be affected (blood pressure, cholesterol, weight).

In many cases, you are interested in the correlations between the groups. For example, the following PROC CORR analysis uses the VAR statement to specify four variables in one group and the WITH statement to specify three variables in another group. The procedure computes the 12 correlations between the two sets of variables:

/* default handling of missing values is pairwise */
proc corr data=Sashelp.Heart noprob nosimple; 
var AgeCHDDiag Height Weight MRW;
with Diastolic Systolic Cholesterol;

Several variables in the Sashelp.Heart data set contain missing values. The table shows the correlations (top of each row) and the number of nonmissing values (bottom of each row) that are used to compute each correlation. The default handling of missing values in PROC CORR is pairwise exclusion, which is why different cells in the table report different numbers of nonmissing values.

The same computation in SAS/IML

This article was motivated by a question that a SAS customer asked on a discussion forum. The question was "How can you do this computation in SAS/IML?" The customer noted that the CORR function in SAS/IML accepts only a single data matrix.

Nevertheless, you can use the CORR function for this computation because the pairwise correlations that you want are a submatrix of the full correlation matrix. That means that you can compute the full correlation matrix (with pairwise exclusion of missing values) and use array indices to extract the relevant submatrix, as follows:

proc iml;
varNames = {"AgeCHDDiag" "Height" "Weight" "MRW"};
withNames = {"Diastolic" "Systolic" "Cholesterol"};
names = varNames || withNames;
use Sashelp.Heart;        
read all var names into X;               /* read VAR and WITH variables */
corr = corr(X, "pearson", "pairwise");   /* compute for correlations */
vIdx = loc( element(names, varNames) );  /* columns of 1st set of variables */
wIdx = loc( element(names, withNames) ); /* columns of 2nd set of variables */
withCorr = corr[wIdx, vIdx];             /* extract submatrix */
print withcorr[r=withNames c=varNames];

Notice that the column indices for the VAR and WITH variables are computed by using the LOC-ELEMENT technique, although for this example the columns are simply 1:4 and 5:7. The correlations are the same as were produced by PROC CORR.

However, you might be concerned that we are doing more work than is necessary. The CORR function computes the entire 7x7 matrix of correlations, which requires computing 28 pairwise correlations, when all we really need is 12. The computation is even less efficient if the VAR and WITH lists are vastly different in size, such as one variable in the first list and six in the second.

To be more efficient, you might want to compute only the correlations that you really need. The following module computes only the pairwise correlations that you specify:

/* Compute correlations between groups of variables. Use pairwise exclusion 
   of missing values. Specify the data matrix X and indices for the 
   columns that correspond to the VAR and WITH variables. */
start CorrWith(X, varIdx, withIdx, method="Pearson");
   withCorr = j(ncol(withIdx), ncol(varIdx));
   do i = 1 to ncol(withIdx);
      w = withIdx[i];               /* column for WITH var */
      do j = 1 to ncol(varIdx);
         v = varIdx[j];             /* column for VAR var */
         withCorr[i,j] = corr(X[,v||w], method, "pairwise")[1,2];
withCorr = CorrWith(X, vIdx, wIdx);
print withcorr[r=withNames c=varNames];

The output is identical to the previous output, and is not shown.

Post a Comment

Create heat maps with PROC SGPLOT

When SAS 9.4m3 was released last month (including SAS/STAT and SAS/IML 14.1), I was happy to see that a HEATMAP statement had been added to the SGPLOT procedure. Although heat maps in the SAS/IML language have been available for several releases, you previously had to use the Graph Template Language (GTL) to create a customized heat map in Base SAS.

To illustrate using the HEATMAP statement, consider the following counts for patients in a medical study, classified by ordinal variables that represent the patients' weight and smoking status. The data and the techniques used to create the table are discussed in the article "How to order categories in a two-way table with PROC FREQ." You can download the complete SAS program that is used to create the table and graphs in this analysis.


The categories and counts are contained in a data set called FreqOut, which was created by the FREQ procedure. The data set is arranged in "list form" (or "long form") as follows:


A basic heat map of a frequency table

The following statements create a basic heat map of the frequencies for the 15 cells of the table:

proc sgplot data=FreqOut;
   heatmap x=Weight_Cat y=Smoking_Cat / freq=Count 
         discretex discretey
         colormodel=TwoColorRamp outline;

On the HEATMAP statement, the DISCRETEX and DISCRETEY options are used because the underlying variables are numeric; a user-defined format is used to render the numbers as "Underweight," "Normal," "Overweight," and so forth. A two-color gradient ramp is good for showing counts. The default color ramp has three colors (blue, white, and red), which is good for visualizing deviations from a reference value.

The heat map makes it visually clear that most patients in the study are overweight non-smokers. Other categories with many patients include the overweight heavy-smokers, the normal-weight non-smokers, and the normal-weight heavy-smokers. Underweight patients (regardless of smoking status) are relatively rare.

Overlaying text on a heat map

Although this same heat map can be created in SAS/IML by using the HEATMAPCONT subroutine, using PROC SGPLOT makes it easy to overlay other plot types and change properties of the graph. For example, you might want to overlay the counts of each cell to create a color-coded table. You can easily add the TEXT statement to the previous SGPLOT statements. If you add the counts, you don't need to include the gradient color ramp. Furthermore, to make the heat map look more like the table, you can add the REVERSE option to the YAXIS statement. Lastly, the category labels make it clear that the variables are related to smoking and weight, so you can suppress the variable names. The following statements implement these revisions and create a new heat map:

proc sgplot data=FreqOut noautolegend;
   heatmap x=Weight_Cat y=Smoking_Cat / freq=Count 
         discretex discretey
         colormodel=TwoColorRamp outline;
   text x=Weight_Cat y=Smoking_Cat text=Count / textattrs=(size=16pt);
   yaxis display=(nolabel) reverse;
   xaxis display=(nolabel);

I'm excited by the HEATMAP statement in PROC SGPLOT. I think it is a great addition to one of my favorite SAS procedures.

Post a Comment

The drunkard's walk in 2-D

Last month I wrote about how to simulate a drunkard's walk in SAS for a drunkard who can move only left or right in one direction. A reader asked whether the problem could be generalized to two dimensions. Yes! This article shows how to simulate a 2-D drunkard's walk.

In fact, it is possible to simulate the drunkard's walk in d dimensions. The drunkard starts at the origin of the coordinate system. He chooses a random coordinate direction and takes a unit step in either the positive or negative direction for that coordinate. If the drunkard takes N steps, you can determine the probability that the drunkard is in a particular location.

The computational methods used to simulate the drunkard's walk in higher dimensions are similar to the 1-D walk, so see my previous article for the background information.

Visualizing a drunkard's walk

In two dimensions, you can use a series plot to visualize the path of the drunkard as he stumbles to the north, south, east, and west. The following PROC IML program uses the Table distribution to draw 10 random values from the set {1, 2}, which are the two coordinate directions. The SAMPLE function makes random draws from the set {-1, 1} and determines whether the drunkard steps forward or backward in the chosen coordinate direction:

proc iml;
call randseed(54321);
dim = 2;                   /* number of directions */
NumSteps = 10;             /* number of steps */
coord = j(NumSteps,1,.);   /* allocate vector for random direction */
call randgen(coord, "Table", j(1,dim,1)/dim); /* random directions */
incr = sample({-1 1}, 1 // NumSteps);         /* random step +/-1 */
step = T(1:NumSteps);
print step coord incr;

The value coord=1 means that the drunkard takes a step in the x direction. The value coord=2 means that the drunkard steps in the y direction. The table shows that the drunkard's first step is a backward step in the x direction (west). The next step is a negative step in the y direction (south), and so forth. The complete trajectory is west, south, south, east, west, north, north, east, south, and east.

These two random vectors can be combined to form an N x d matrix where each column is a coordinate direction, each row is a time step, and the values 1, -1, and 0 indicate that the drunkard moves forward, backward, or does not move, respectively. This matrix contains the same information as the previous table. You can use the SUB2NDX function to convert subscripts into indices, and therefore create the matrix, as follows:

x = j(NumSteps, dim, 0);      /* (N x d) array */
idx = sub2ndx(dimension(x), step||coord);
x[ idx ] = incr;
print x;

You can encapsulate the preceding statements into a module that implements the Drunkard's Walk algorithm. You can download the complete SAS program that generates the graphs and computations in this article.

To visualize the path, you can use the CUSUM function on each column of the x matrix to accumulate the steps. The accumulation creates a new matrix that contains the position of the drunkard after each step. The following series plot attempts to show the path of the drunkard, with locations labeled by time. For most random walks, there is considerable overplotting. The graph show that the drunkard started at the origin and sequentially moved west, south, south, east, west, and so forth. The final position of the drunkard is at (1, -1). This path is typical in that the drunkard staggers near the origin. Although it is possible for the drunkard to always choose the same direction, it is not probable.


Visualizing the probability distribution of the final position

Just as in the one-dimensional case, you can use simulation to visualize the probability distribution of the final position of the drunkard. If the drunkard takes an even number of steps, he must end up an even distance away from the origin, where the distance is calculated by using the Minkowski metric, also known as the "Manhattan" or "taxicab" metric. The probability is highest that the drunkard returns to the origin, followed by the points that are two units away, four units away, and so forth. The following SAS/IML statements call the DrunkardWalk function 10,000 times and plots the counts for each final position by using a heat map:

/* Assume DrunkardWalk function is defined (see link to program).
   Compute final positions of 10,000 walks */
NumSteps = 6;
NumWalks = 10000;
finalPos = j(NumWalks, dim);
do i = 1 to NumWalks;
   x = DrunkardWalk(dim, NumSteps);
   finalPos[i,] = x[+,];
create walk from finalPos[c={"x" "y"}];   /* write to SAS data set */
append from finalPos;
   proc freq data=walk noprint;           /* PROC FREQ computes counts */
   tables x*y / list sparse out=freqOut;
/* read results back into SAS/IML vectors */
use freqOut;   read all var {x y count};   close;
/* reshape counts into rectangular matrix and visualize with heat map */
nX = range(x) + 1;
nY = range(y) + 1;
W = shape(count, nX, nY);
call HeatmapCont(W) xValues = unique(x) yValues=unique(y);


The heat map estimates the probability distribution of the final position of the drunkard after taking six steps. Notice that the drunkard must land on a spot that is an even number of units (in the Minkowski metric) from the origin. About 10% of the time the drunkard returned to the origin. About 7% of the time, the drunkard ends at (1,1), (-1,1), (-1,-1), or (1,-1). About 5% of the time the drunkard ends at (2,0), (0,2), (-2,0), or (0,-2). Positions that are four units away are less probable. Final positions that are six units away occurred during the simulation, but were so rare that they are barely visible in the heat map. For example, only 14 walks ended at the location (0, 6).

The drunkard's walk in two or more dimensions is a fun simulation of a discrete distribution. SAS/IML provides the necessary functionality to simulate many walks and to analyze the distribution of the results.

Post a Comment

Those tricky PERCENT formats

When using SAS to format a number as a percentage, there is a little trick that you need to remember: the width of the formatted value must include room for the decimal point, the percent sign, and the possibility of two parentheses that indicate negative values. The field width must include room for the parentheses even if the number that you are formatting is strictly positive!

The PERCENTw.d format

The PERCENTw.d format takes a proportion and displays it as a percentage. For example, 0.123 gets displayed as 12.3%. Negative values are displayed by enclosing them in parentheses, so -0.123 is displayed as (12.3%). Thus for proportions in the range (-1, 1), you need to allow at least six spaces: one each for the decimal point and percent sign, and two each for the number and the parentheses. This leads to the following heuristic:

For the PERCENTw.d format, the field width should be at least 6 more than the number of decimal places that you want to display.

In other words, for the PERCENTw.d format, a good rule of thumb is to set w = d + 6. Smaller field widths will either steal space from the decimal width or display the value as "*%", which means "unable to display this value by using the current field width."

My favorite formats are PERCENT7.1 and PERCENT8.2. An exception to the above heuristic is that I sometimes use PERCENT5.0 for variables that appear along the axis of a graph; you can use 5 for the field width because there is no decimal point to display. The following SAS DATA step displays several positive and negative proportions in the PERCENT7.1 and PERCENT8.2 formats:

data Percent;
format x BEST8. 
       p1 PERCENT7.1
       p2 PERCENT8.2;
label x = "Raw value"
      p1 = "PERCENT7.1"
      p2 = "PERCENT8.2";
input x @@;
p1 = x; 
p2 = x; 
0.7543  0.0123 -0.0123 -0.7543
proc print data=Percent noobs label; run;

If you have values that are greater than 1, you will need to increase the field width, since these will be displayed as 150% (1.5) or 2500% (25).

The PERCENTNw.d format

Not every SAS programmer knows about the PERCENTNw.d format. Notice the extra 'N' in the name, which stands for 'negative.' This format is similar to the PERCENTw.d format except that negative values are displayed by using a negative sign, rather than parentheses.

For many applications, I prefer the PERCENTNw.d format because of the way that it handles negative values. The heuristic for specifying the field width is the same as for the PERCENTw.d format. This might be surprise you. Since the PERCENTNw.d format uses a negative sign instead of parentheses, shouldn't you be able to specify a smaller field width?

No. For ease of converting between the PERCENTw.d and PERCENTNw.d formats, the rules are the same. The documentation for the PERCENTNw.d format states that "the width of the output field must account for the minus sign (–), the percent sign (%), and a trailing blank [emphasis added], whether the number is negative or positive."

You can modify the previous DATA step to display the same raw values by using the PERCENTNw.d format. The results follow:


In summary, remember to pad the field width for these formats by including space for the decimal point, the percent sign, and the possible negative sign or parentheses. You must do this even if your values are positive.

This blog post was inspired by an exchange on the SAS-L discussion forum. If you are new to Base SAS programming, The "L" is a good place to lurk-and-learn or to ask questions. I also recommend the SAS Support Communities for learning SAS software and related topics such as statistics, data mining, and matrix computations.

Post a Comment