A variation of the Pareto chart (Wilkinson, TAS, 2006) uses confidence intervals for multinomial counts to enhance the information on the chart. This article supplies the necessary background for how to estimate confidence intervals for multinomial counts. A subsequent article shows how to construct Wilkinson's Pareto chart in SAS.

The multinomial distribution

Suppose a categorical variable, X, can take on k different values. A statistical model of the variable is to assume that X can be modeled by a multinomial distribution. A multinomial distribution assumes that the probability of the i_th value is p_i, where Σ_i p_i = 1.

In a sample of size N, there will be N₁ items of the first type, N₂ items of the second type, and so forth, where Σ_i N_i = N. The k-tuple of counts (N₁, N₂, ..., N_k) is a single observation from the multinomial distribution.

The expected value of N_i is N*p_i, but in an actual sample the count could be higher or lower than that value. Confidence intervals enable you to answer the question, "what is the likely range of each count?"

A previous article describes how to construct simultaneous confidence intervals for the proportions. But the Wilkinson variation of the Pareto chart uses marginal confidence limits, which are different. Marginal confidence limits can be easily generated by running a simulation of the multinomial distribution and using percentiles of the results. In SAS, the easiest way to simulate from the multinomial distribution is to use the RandMultinomial function in SAS IML.

Simulating multinomial frequencies

Suppose a quality engineer samples 40 defective items. She finds that there are six causes for the defective items. In the sample, N₁=17 items are defective because of the first cause, N₂=10 are defective because of the second cause, and the remaining causes are responsible for 5, 3, 3, and 2 counts, respectively. The best estimate of the unknown probabilities is obtained by dividing the counts by N=40. Thus, the empirical probabilities are p = {0.425, 0.25, 0.125, 0.075, 0.075, 0.05}.

The following call to PROC IML specifies these parameters and generates five random samples. Each sample is a random draw from the multinomial distribution.

/* CI for marginal probability parameters in a multinomial distribution */
proc iml;
/* observed frequencies */
freq = {17, 10, 5, 3, 3, 2};
N = sum(freq);      /* sample size */
p = freq / N;
 
/* 5 random samples of size N from Multinom(p) */
call randseed(123);
X5 = randmultinomial(5, N, p);
print X5[c=('C1':'C6') r=('S1':'S5')];

Each row shows the count in a random sample. The i_th column shows the counts for the i_th defective cause across all five samples. In this small simulation, the counts for the first cause range from 11-18, the counts for the second cause range from 7-11, and so forth. Because the probabilities for the 4th-6th causes are relatively small, some samples did not contain any items for those causes.

The small output gives us the intuition behind confidence intervals. We'd like to know ranges that cover the counts for 95% of the samples. You can do this by generating a large number of samples and calculating percentiles.

Confidence intervals for the multinomial counts

The previous output shows the range of counts for each variable for five random samples. Now imagine generating many more random samples from the same multinomial distribution. If you compute the 2.5th and 97.5th percentiles for each column, you obtain a simulation-based estimate of a 95% confidence interval for each column's values. The following IML statements perform these calculations for a simulation that uses 10,000 random draws.

/* Use simulation to estimate (1-alpha)*100% confidence intervals 
   for counts for each category.
   - N = sample size
   - p = parameter vector for Multinomial(N, p) distribution
   - B is the number of Monte Carlo simulations
   - alpha = significance level (default 0.05) ==> 95% CL   
*/
start MultinomialCL(B, N, p, alpha=0.05);
   X = RandMultinomial(B, N, p); /* generate B random samples */
   /* Find the (alpha / 2) and (1 - alpha / 2) percentiles.
      For alpha = 0.05, this estimates the 2.5th and 97.5th percentiles. */
   w = (alpha/2) // (1-alpha/2);
   call qntl(CL, X, w);       /* QNTL computes the quantiles of each column */
   return( CL );
finish;
 
/* call function and print the results in a table */
limits = MultinomialCL(10000, N, p);
Category = 'C1':'C6';
print limits[r={'Lower95' 'Upper95'} c=Category L="Marginal Limits (N=40)"];

The table shows 95% confidence limits for the counts from a multinomial sample of size N=40 that has a vector of probability parameters p. For the most probable category (C1), you can expect the count to be between 11 and 23 for 95% of samples. For the least probable category (C6), you can expect the count to be between 0 and 5.

Although the table provides all relevant information, I like to graph the confidence intervals as bars in a chart that graphs the observed frequency and confidence intervals versus the categories, as follows:

/* graph the empirical frequencies and the marginal CLs */
Lower95 = limits[1,];
Upper95 = limits[2,];
create MultinomialLimits var {'Category' 'Freq' 'Lower95' 'Upper95'};
   append;
close;
QUIT;
 
title "95% Confidence Limits for Frequencies";
proc sgplot data=MultinomialLimits noautolegend;
   highlow x=Category low=Lower95 high=Upper95 / 
           type=bar barwidth=0.8;* primary=true;
   scatter x=Category y=Freq;
   yaxis offsetmin=0 grid;
   xaxis type=discrete discreteorder=data;
run;

The graph enables you to see at a glance the ranges of counts for a sample of size N from the multinomial(N, p) distribution.

Summary

The statistics of small random samples can be highly variable. This article looks at the counts in a multinomial distribution. By running a simple simulation study, you can construct 95% confidence intervals for the counts. This is helpful for quality engineers who need to know whether an increase or decrease in an observed count is associated with a change in the manufacturing process or is merely random variation. These confidence intervals are used in constructing Wilkinson's variation of the Pareto chart, which I will construct in a future article.

The post Confidence limits for counts in a multinomial distribution appeared first on The DO Loop.

However, there are instances in which a numerical matrix is merely a way to organize a bunch of numbers. For example, in a simulation study, you might construct a matrix whose i_th column contains the generated data that results from the i_th independent random sample. In this instance, sorting each column of the matrix is useful in analyzing the sampling distribution of the quantiles. For example, if you sort all columns, then the first row of the sorted matrix contains the distribution of the minimum statistic, and the last row contains the distribution of the maximum statistic.

This article shows how to sort the rows and columns of a matrix independently in the SAS IML matrix language.

A matrix from a simulation study

Suppose you want to run a simulation study to approximate the sampling distribution of quantiles in small, normally distributed, random samples. You decide to study samples of size N=25 and you want to graph the approximate distribution of the quantiles based on generating B samples from N(0,1).

In a matrix language like SAS IML, it is efficient to simulate all random samples by using a single call to the RANDFUN functions. You can fill an N x B matrix with N(0,1) variates and interpret each column as a random sample of size N. In a real simulation study, you would choose a large number for B (such as 10,000), but for this article I will choose only B=200 samples because I want to show some visualizations by using a heat map. The following SAS IML program generates the data for the simulation study:

proc iml;
/* simulation study: If X ~ N(0,1), 
   what is the distribution of the sample quantiles? */
N = 25;                      /* sample size */
B = 200;                     /* number of simulations of X ~ N(0,1) */
call randseed(1234);
M = randfun(N//B, "Normal"); /* NxB matrix; each col is random sample */

If you sort each column in the matrix, then it is easier to compute certain rank-based statistics and quantiles. For example, in the sorted matrix, the first row will contain the minimum values for each sample. Similarly, the last row contains the maximum values, and the 13th row contains the median values. (You could also use the QNTL function, which computes quantiles for each column in a matrix.)

Sorting each column of a matrix

If you use CALL SORT in IML (or the SORT procedure in Base SAS), the rows of the matrix are arranged in order according to one or more key columns. To sort the columns independently, you need to loop over the columns, sort each one, and then overwrite the column with the sorted version.

The following IML module, SortMat, can sort the columns or rows of the matrix in either ascending or descending order. The first argument should be either "row" or "col". The second argument is the matrix to sort. The third argument is optional. If specified, use 0 to sort the rows or columns in ascending order and use 1 to sort in descending order.

/* sort the columns of a matrix independently */
start SortCols(M, descend=0);
  S = M;
  if ^descend then do;
     do i = 1 to ncol(S);
        x = S[,i];
        call sort(x);        /* sort each column in ascending order */
        S[,i] = x;
     end;
  end;
  else do;
     do i = 1 to ncol(S);
        x = S[,i];
        call sort(x, 1, 1); /* sort each column in descending order */
        S[,i] = x;
     end;
  end;
  return( S );
finish;
 
/* sort the rows or columns of a matrix independently.
   SYNTAX:
   A_row = SortMat("row", A <,descend=0>)
   A_col = SortMat("col", A <,descend=0>)
*/
start SortMat(direction, M, descend=0);
  isRow = (ksubstr(upcase(direction), 1, 3) = "ROW");
  if isRow then 
     return( T(SortCols(M`, descend)) );
  else 
     return( SortCols(M, descend) );
finish;
store module=(SortCols SortMat);
 
/* sort matrix by cols and make a heat map to visualize the result */
t = "Matrix Sorted Ascending by Cols";
A = SortMat("col", M);
ods graphics / width=480px height=240px;
call heatmapcont(A) colorramp="ThreeColor" displayoutlines=0 xaxistop=1 title=t;
/* Note: To sort the columns in descending order, use 
   A = SortMat("col", M, 1);
*/

The heat map visualizes the column-wise sorting. For each column, the values in the sorted matrix are stored in increasing order. Thus, the first row contains the smallest value in each sample, the second row contains the second smallest value, and so on, until the last row, which contains the largest value in each sample. Equivalently, each column contains the order statistics for the sample.

If you want to analyze, for example, the distribution of the minimum values in the samples, you can extract the first row. The following statements extract the first row and create a histogram of the sampling distribution of the minimum statistic:

/* The smallest values are in the first row of the sorted matrix */
sample_min = A[1,];
title "Distribution of the Sample Min for N(0,1) Data";
call histogram(sample_min);

The histogram shows the distribution of the minimum statistic. This distribution is one of the extreme-value distributions, known as the Gumbel distribution. For a standardized normal sample of this size, it is common that the minimum value is near -2, but values less than -3 are possible.

Sorting each row of a matrix

In the same way, you can pass "row" as the first argument of the SortMat function to independently sort each row of a matrix:

/* Suppose you organize the simulation so that each row is a random sample. Then the
   simulated samples are the rows of M`. You can sort the rows by calling the SortMat("row",...) function. */
M = M`;
t = "Matrix Sorted Ascending by Rows";
B = SortMat("row", M);
ods graphics / width=280px height=480px;
call heatmapcont(B) colorramp="ThreeColor" displayoutlines=0 xaxistop=1 title=t;
/* To sort in descending order, use
B = SortMat("row", M, 1);  */

Instead of generating new data, I simply transposed the previous simulated matrix. For the transposed matrix, each row is a random sample. If you sort by rows, the columns contain the distributions of the quantiles. The i_th column contains the i_th smallest value in each sample.

Summary

This article creates a SAS IML function that can independently sort each column or each row of a matrix. One application is that you can easily analyze the distribution of order statistics, including minima, maxima, and other quantiles. Another application is to sort the columns of a lasagna plot, which can visualize changes to the distribution of a response variable over time.

The post Sort the rows or columns of a matrix independently appeared first on The DO Loop.

There are two common variations of the Pareto plot. The most familiar is a bar chart of frequencies that is overlaid with a cumulative frequency curve. This graph is shown to the right. The other is a "cascade chart" (sometimes called a waterfall chart) that displays the cumulative frequencies as a staircase of bars.

These charts are produced automatically by the PARETO procedure in SAS/QC software. The PARETO procedure has many options for displaying Pareto charts. However, not every SAS user has a license for SAS/QC. This article shows how to create basic Pareto charts by using Base SAS procedures and the SGPLOT procedure.

The data

To illustrate a Pareto chart, let's use a modification of data from the documentation of the PARETO procedure. A random sample of 40 defective parts are examined, and the cause of each defect is recorded, as follows:

/* Modification of PROC PARETO example */
data Failure;
INFILE DATALINES delimiter=',';
length Cause $ 16;
label Cause = 'Cause of Failure';
input Cause @@;
datalines;
Corrosion, Oxide Defect, Contamination, Oxide Defect
Oxide Defect, Oxide Defect, Contamination, Metallization
Oxide Defect, Contamination, Contamination, Oxide Defect
Contamination, Contamination, Contamination, Corrosion
Silicon Defect, Contamination, Contamination, Contamination
Contamination, Contamination, Doping, Oxide Defect
Oxide Defect, Metallization, Contamination, Corrosion
Silicon Defect, Contamination, Corrosion, Corrosion
Metallization, Oxide Defect, Contamination, Contamination
Oxide Defect, Doping, Doping, Contamination
;

If you have a license for SAS/QC software, you can create a basic Pareto chart with an overlaid cumulative curve by using the PARETO procedure. So that you can easily reuse the SAS code for your own data, I have defined macro variables for the name of the data set and the name of the categorical variable.

%let DSName = Failure;
%let VarName = Cause;
proc pareto data=&DSName;
   vbar &VarName;
run;

The graph is shown at the top of this article. It shows that "Contamination" is responsible for 42.5% of the failures. The first three categories account for 80% of the defects, so those are the issues that quality engineers should strive to improve. Whereas the classic "Pareto rule" is that 20% of the issues result in 80% of the defects, the cumulative curve shows that, for this process, 50% of the issues (the first three categories) are responsible for about 80% of the defects.

A Pareto chart in Base SAS

If you do not have a license for SAS/QC software, you can still create a basic Pareto chart by using PROC FREQ to compute the frequencies and PROC SGPLOT to display a sorted bar chart and a cumulative curve. There are a few useful options to know about:

1. The PROC FREQ statement supports the ORDER=FREQ option, which sorts the frequencies in descending order.
2. The TABLE statement supports the OUT= option, which you can use to create a data set that contains the data for the Pareto chart. Use the OUTCUM option to output the cumulative statistics.
3. The SGPLOT procedure supports the VBAR and VLINE statements, which enables you to overlay a bar chart and a cumulative curve. You can use the Y2AXIS option on the VLINE statement to display a separate axis for the cumulative scale.

/* create the Pareto charts by hand in SGPLOT */
proc freq data=&DSName order=Freq noprint;
   where not missing(&VarName);
   table &VarName / out=FreqOut outcum;
run;
 
title "Standard Pareto Chart";
proc sgplot data=FreqOut noautolegend;
   vbar &VarName / response=Percent;
   xaxis type=discrete discreteorder=data;
   yaxis grid min=0 max=100 offsetmin=0 label="Percent";
   /* overlay the cumulative percentage on the Y2 axis */
   vline &VarName /response=cum_Pct markers datalabel y2axis;
   y2axis min=0 max=100 offsetmin=0 label="Cumulative Percent";
   format cum_Pct best4.;
run;

The output is similar to the previous Pareto chart. I decided to add grid lines and to add markers to the cumulative curve. Otherwise, the information in this chart is the same as the one produced by PROC PARETO.

Problem with the classical Pareto chart

The classical Pareto chart suffers from three problems (Wilkinson, TAS, 2006):

1. Dual scale: Mathematically, the cumulative curve is the integral (or sum) of the bar heights. It can be confusing to plot both quantities in the same graph. The "dual scale" problem is even worse if you show the scale of the bars as counts instead of percentages, as is often done.
2. Range: When there are many categories (or when no category is responsible for most of the defects), the bars are short, whereas the cumulative curve will always range to 100%. In that situation, it is difficult to judge the height of the bars, which get squashed down to the bottom of the frame.
3. Interpolation: It is wrong to use line segments to connect the cumulative frequencies. It makes it seem like the cumulative probability is a continuous function that increases linearly between categories. It is not. For each category, it is a single number. The number increases discontinuously. The categories are discrete so there is nothing "between" them.
4. Reference distribution: Implicitly, calling something a "Pareto chart" implies that the data follows a Pareto distribution. In a Pareto distribution, the frequencies follow a power law. The assumption of an underlying power law justifies focusing on the issue that has the largest frequency. If the categories are equally likely to occur, you might as well focus on the issues that are cheapest or quickest to resolve.

A Pareto cascade chart

An alternative to the classical Pareto chart addresses the first three problems. I call the revised chart a cascade chart or a cumulative Pareto chart. A bar chart places the base of each bar along a common baseline at zero. In contrast, a cascade chart aligns the base of each bar with the top of the preceding bar. This results in a staircase-shaped graph. The height of each "step" is proportional to the relative frequency of each category. The tops of the steps show the cumulative distribution of frequencies.

If you have access to the PARETO procedure, you can use the CHARTTYPE=CUMULATIVE option to create a cascade chart, as follows:

proc pareto data=&DSName;
   vbar &VarName / charttype=cumulative;
run;

In this graph, the tops of the "steps" show the cumulative distribution of frequencies. The heights of the steps show each category's contributions.

You can also create the chart in Base SAS. First, run the PROC FREQ step that was shown earlier. This creates the FreqOut data set. The following DATA step uses the LAG function to compute the base level for each step. You can then use a high-low plot to display the staircase-shaped cascade chart:

/* add upper and lower variables for the HIGHLOW plot */
data CumPareto;
set FreqOut;
_Lower = lag(cum_Pct);
if _N_ = 1 then 
   _Lower = 0; 
run;
 
title "Cumulative Pareto Chart";
proc sgplot data=CumPareto;
   highlow x=&VarName low=_Lower high=cum_Pct / 
           type=bar barwidth=1 highlabel=Percent;
   yaxis offsetmin=0 grid;
   xaxis type=discrete discreteorder=data;
   format Percent best4.;
run;

Again, I have added gridlines and labels to the bars, but otherwise the information is the same as the cumulative chart that is created by PROC PARETO. This chart has only a single axis, which shows the cumulative probability. Short bars are no longer squashed at the bottom of the chart. By adding labels, you can see the size of the relative frequencies and the cumulative distribution on a common scale.

There is a third type of Pareto chart that addresses the problem of not being able to see a reference distribution. I will discuss that variation in a separate article.

Summary

In SAS, the PARETO procedure in SAS/QC software provides options to create many types of Pareto charts. If you do not have a license for SAS/QC software, this article shows how to use Base SAS to create two simple Pareto charts. You can use a Pareto chart to identify the most frequent categories. In practice, the most frequent categories are analyzed by quality engineers in hopes that addressing those issues will lead to a major improvement in the quality of a manufacturing process.

To simplify the creation of basic Pareto charts, I have encapsulated the techniques in this article into two SAS macros, which are shown in the Appendix.

Appendix: SAS macros that create basic Pareto charts

/* Macros to create a basic Pareto chart, written by Rick Wicklin. For details, see
   https://blogs.sas.com/content/iml/2026/06/22/pareto-charts-sas.html
*/
%macro StdPareto(DSName, VarName);
proc freq data=&DSName order=Freq noprint;
   where not missing(&VarName);
   table &VarName / out=_FreqOut outcum;
run;
 
proc sgplot data=_FreqOut noautolegend;
   vbar &VarName / response=Percent;
   xaxis type=discrete discreteorder=data;
   yaxis grid min=0 max=100 offsetmin=0 label="Percent";
   /* overlay the cumulative percentage on the Y2 axis */
   vline &VarName /response=cum_Pct markers datalabel y2axis;
   y2axis min=0 max=100 offsetmin=0 label="Cumulative Percent";
   format cum_Pct best4.;
run;
%mend;
 
%macro CumPareto(DSName, VarName);
   proc freq data=&DSName order=Freq noprint;
      where not missing(&VarName);
      table &VarName / out=_FreqOut outcum;
   run;
 
   /* add upper and lower variables for the HIGHLOW plot */
   data _CumPareto;
   set _FreqOut;
   _Lower = lag(cum_Pct);
   if _N_ = 1 then 
      _Lower = 0; 
   run;
 
   proc sgplot data=_CumPareto;
      highlow x=&VarName low=_Lower high=cum_Pct / 
              type=bar barwidth=1 highlabel=Percent;
      yaxis  offsetmin=0 grid;
      xaxis type=discrete discreteorder=data;
      format Percent best4.;
   run;
%mend;
 
/* show how to call each macro. The categorical variable can be 
   numeric or character */
title;
/* character variable */
%StdPareto(sashelp.cars, Type);
%CumPareto(sashelp.cars, Type);
/* numeric variable */
%StdPareto(sashelp.cars, Cylinders);
%CumPareto(sashelp.cars, Cylinders);

The post Create two types of Pareto charts in SAS appeared first on The DO Loop.

A previous article shows how to construct confidence bands for the CDF of continuous univariate data. The bands can be added to a plot of the empirical CDF (ECDF) for the data. One of the drawbacks of an ECDF plot is that the cumulative distribution is an S-shaped curve. If you display two S-shaped curves that represent cumulative distributions, it is difficult for the human eye to detect differences between them, which is why statisticians plot histograms (estimates of the probability density) more often than ECDF curves (estimates of the cumulative distribution).

If you want to compare the data distribution to a theoretical parametric distribution (such as normal or Weibull), the quantile-quantile (Q-Q) plot is a useful alternative to the CDF plot. A Q-Q plot is a visual representation of a statistical goodness-of-fit test: it helps you assess whether the data are a random sample from a theoretical distribution. Data that are sampled from the specified parametric distribution appear to be linear in a Q-Q plot.

More correctly, they appear to be linear for most random samples. As discussed in a previous article, a Q-Q plot might look nonlinear due to random variation, especially for a small sample. Consequently, it can be useful to add a confidence band to the Q-Q plot.

The relationship between the CDF plot and the Q-Q plot

It turns out that you can transform the ECDF and confidence bands from the previous article to "straighten out" the cumulative distributions. The result is a Q-Q plot with a confidence band.

Recall that if X is distributed according to some distribution, F, then F(X) is uniformly distributed. And if U is uniformly distributed, then F^-1 is distributed according to F. Details and an example are discussed in a previous article on the probability integral transformation. You can use this result to "straighten out" the ECDF and the confidence bands, provided that you choose a distribution to use for the transformation. I will transform the simple Kolmogorov band from the previous article, but you can apply this technique to any other confidence band.

The most familiar distribution is the standard normal distribution. The standard cumulative distribution function is usually denoted as F = Φ. In SAS, you can apply Φ by calling the function CDF("Normal", x). You can apply the inverse transformation, Φ^-1, by calling the SAS function QUANTILE("Normal", p). The inverse normal CDF transformation is shown to the right. Notice that the function is very steep when the probability is near 0 or 1. This means that a small change in probability results in a very large change in the associated quantiles.

Constructing the Q-Q Plot Confidence Bands in SAS IML

For data, let's analyze the breaking strength data of fiber-optic cords from previous posts.

A previous article shows how to construct a normal Q-Q plot in SAS by using the following steps:

1. Sort the data.
2. Compute n evenly spaced points in the interval (0,1), where n is the number of data points in your sample. SAS procedures often use Blom's formula, where the i_th point is v_i = (i - 0.375) / (n + 0.25).
3. Compute the quantiles (inverse CDF) of the evenly spaced points. This gives you quantiles of the standard normal distribution.
4. Create a scatter plot of the sorted data versus the standard quantiles computed in Step 3. It is traditional to put the data on the vertical axis and the theoretical quantiles on the horizontal axis. This is in contrast to histograms and ECDF plots, which place the data on the horizontal axis.
If you want to add confidence bands to the Q-Q plot, you can transform confidence bands for the CDF. The transformation maps probabilities in (0,1) into standardized quantiles.

The following SAS IML program uses the ECDF and ECDF_KSCL modules from earlier articles. For your convenience, you can download these functions from GitHub. The program compute Kolmogorov bands for the ECDF, transforms them into quantile scale, and writes the results to a data set for graphing.

/* Before running this program, STORE the ECDF and ECDF_KSCL modules and define the 
   Cord data set. See
   https://blogs.sas.com/content/iml/2026/05/26/create-ecdf.html
   https://blogs.sas.com/content/iml/2026/06/08/confidence-bands-ecdf.html
*/
proc iml;
/* In SAS Viya, the ECDF function is built into SAS IML and does not need to be loaded. */
load module=(ECDF ECDF_KSCL);  
 
/* Read and sort the data */
use Cord;  read all var "Strength" into x;  close;
call sort(x);
n = countn(x);
 
/* Compute theoretical plotting positions and map to standardized quantiles (Blom, 1958) */
v = ((1:n) - 0.375) / (n + 0.25);  
q = quantile("Normal", v);
 
/* Compute the ECDF and the Kolmogorov probability bands */
y = ECDF(x);
KS_band = ECDF_KSCL(y);   /* 95% CL */
 
/* Optionally, transform probability bands to quantile bounds.
   The domain of the QUANTILE function is the open interval (0,1),
   so clip the band values.
*/
F_Lower = choose(KS_band[,1] > 0, KS_band[,1], .);
F_Upper = choose(KS_band[,2] < 1, KS_band[,2], .);
/* For a normal Q-Q plot, use the quantile of the normal distribution.
   Use the quantile function for other distributions (e.g., "Exponential") to 
   construct confidence intervals for other Q-Q plots. */
Q_Lower = quantile("Normal", F_Lower);
Q_Upper = quantile("Normal", F_Upper);
 
/* Create a dataset for plotting */
create QQ_Bands var {"x" "q" "Q_Lower" "Q_Upper"};
append;
close;
QUIT;
 
/* Create the Q-Q plot. The data is plotted on the vertical axis.
   The theoretical quantiles are plotted on the horizontal axis. */
title "Normal Q-Q Plot with 95% Confidence Bands";
proc sgplot data=QQ_Bands noautolegend;
   /* Draw the confidence bands */
   series x=Q_Lower y=x / lineattrs=(color=gray);
   series x=Q_Upper y=x / lineattrs=(color=gray);
   /* Overlay a scatter plot of the quantiles of the data vs the standard normal quantiles. */
   scatter x=q y=x;
   xaxis label="Theoretical Normal Quantiles" grid;
   yaxis label="Sample Quantiles (Strength)" grid;
run;

For your convenience, I have also written an IML function (QQ_KSCL) that encapsulates the statements that generate the Q-Q plot and confidence bands.

Interpreting the Q-Q Plot Bands

Intuitively, the Q-Q plot and its confidence bands are a visual depiction of a hypothesis test. The Q-Q plot indicates whether the data might be a random sample from a normal distribution. If the scatter plot is approximately linear, then the assumption of normality seems reasonable.

But what does "approximately linear" mean? By how much can a Q-Q plot deviate from linearity without rejecting the assumption of normality? The confidence bands indicate that the scatter plot can bend a little bit in the center of the distribution and can bend more severely in the tails. It is said that George Box used to place a fat pencil on a Q-Q plot. If he could adjust the pencil so that it covered all points, then he assumed that the data were approximately normal. The confidence band provides a similar visual indicator. If you can draw a straight line that sits entirely within the confidence bands, the data might be approximately normal. You can run a more formal analysis to test that hypothesis.

The widening of the confidence bands in the extreme quantiles explains why points in the tails of a Q-Q plot often deviate from the linear reference line or exhibit a slight bend. Without confidence bands, you might mistakenly conclude that these deviations are evidence of skewness, heavy tails, or non-normality. The confidence bands indicate that deviations in the tail are expected due to natural sampling variability. If the scatter points are inside the bands, there is insufficient evidence to reject the null hypothesis that the data comes from a normal distribution.

A normal Q-Q plot for non-normal data

Before concluding, let's take a quick look at a similar graph for data that are not normal. I generated a random sample of 50 data from an exponential distribution and ran the same SAS code. The Q-Q plot for the exponential data looks like the following graph:

If you apply Box's "fat pencil" test, it is clear that you cannot cover the scatter points with a pencil. You should conclude that the normal distribution is not a good fit for these data.

The Kolmogorov bands in this plot are so wide that you might actually be able to draw a line that stays inside these confidence bands. However, the scatter points are not near that line. A weakness of the Kolmogorov-Smirnov test for normality (and the Kolmogorov confidence bands) is low statistical power for small samples. The Kolmogorov bounds themselves do not rule out the possibility of fitting a normal distribution to these data, but the nonlinear shape of the data points suggest that the fit is not good.

Summary

The standard normal quantile function is a transformation that maps probabilities onto quantiles. By transforming the non-parametric Kolmogorov bands from an ECDF, you can generate simultaneous confidence bands for a normal Q-Q plot. You could use a different quantile function to obtain Q-Q plots for other distributions.

For moderately sized samples, if you can draw a straight line that stays inside the bands, you might want to perform a formal test for the hypothesis that the data are normally distributed. The Kolmogorov bands have low power, so be careful using them for small samples.

You can download all functions from the ECDF and Q-Q plot articles. You can also download a program that creates all tables and graphs in this series of blog posts.

The post Create a confidence band for a normal Q-Q plot appeared first on The DO Loop.

In statistics, we assume that a sample is a random draw from an underlying distribution. Therefore, an ECDF has a sampling variability. As with other point estimates in statistics, you can construct a confidence band for the ECDF. The confidence band assumes that the samples are drawn randomly and independently from a population. If so, the true CDF of the population is entirely within the confidence bands for 95% of random samples.

This article shows how to use SAS to create a Kolmogorov confidence band for the ECDF. The Kolmogorov is nonparametric, which means it does not assume that the data is generated from some well-known parametric distribution such as Normal, exponential, beta, Weibull, and so forth.

The Kolmogorov confidence bands for an ECDF

One method for constructing confidence bands for an ECDF uses the Kolmogorov-Smirnov (K-S) test statistic, often denoted as D. As shown in a previous article, the K-S statistic measures the maximum absolute distance between the ECDF and the true CDF. By using the asymptotic distribution of the K-S statistic, you can calculate a 95% confidence envelope. The mathematics behind building the confidence band is presented in lecture notes from Pennsylvania State University and in lecture notes by Charles Geyer at U. MN.

To construct the 95% K-S confidence bands, you need the critical value of the Kolmogorov distribution. As noted in a previous article about critical values for the K-S statistic, the two-sided 95% critical value for the asymptotic K-S distribution is approximately 1.36. (More precisely, 1.358099.) The width of the confidence band at a point on the ECDF is the critical value divided by the square root of the sample size: D_crit/√n. Because probabilities cannot be less than 0 or greater than 1, the confidence bands are truncated at these values. (Note: for small samples and for other significance levels, you can use a SAS IML program to compute the exact critical values for any significance level.)

The following SAS IML program calculates a 95% confidence band for the ECDF of the Strength variable on the Cord data set. The program uses the ECDF module from the previous post, and the data from the Cord data set (N=50), which is defined in the same article. So, define the data set and STORE the ECDF function before running the following program.

/* Before running this program in SAS 9.4, define and STORE the ECDF function from 
   https://blogs.sas.com/content/iml/2026/05/26/create-ecdf.html 
   Also define the CORD data set from that blog post. */
proc iml;
load module=(ECDF);
 
/* The asymptotic 95th percentile of the K-S distribution is 1.358099.
   See https://blogs.sas.com/content/iml/2019/05/20/critical-values-kolmogorov-test.html
   or the lecture notes from Charles J. Geyer at U. MN. at
   https://www.stat.umn.edu/geyer/5601/examp/kolmogorov.html
*/
ks_crit = 1.358099;
 
/* read the data for the ECDF */
use Cord;  read all var "Strength" into x;  close;
 
/* evaluate the ECDF at all data points */
call sort(x);
ecdf = ECDF(x);
n = countn(x);   /* count the nonmissing values */
 
/* Compute the standard error for the K-S bands */
SE_KS = ks_crit / sqrt(n);
 
/* construct upper and lower 95% CI bands centered on the ECDF.
   The >< operator returns the minimum, ensuring the max value is 1.
   The <> operator returns the maximum, ensuring the min value is 0.
   See https://blogs.sas.com/content/iml/2026/02/04/clip-values.html */
CL95_Upper = (ecdf + SE_KS) >< 1;
CL95_Lower = (ecdf - SE_KS) <> 0;
 
create ECDF_bands var {"x" "ECDF" "CL95_Upper" "CL95_Lower" };
  append;
close;
QUIT;
 
/* Graph the step functions by using PROC SGPLOT */
title "Empirical CDF and 95% Confidence Bands";
proc sgplot data=ECDF_bands noautolegend;
   label x="Breaking Strength (psi)" ECDF="Cumulative Proportion";
   step x=X y=ECDF / lineattrs=(thickness=2);
   step x=X y=CL95_Lower / lineattrs=(color=gray);
   step x=X y=CL95_Upper / lineattrs=(color=gray);
   xaxis grid label="x";
   yaxis grid min=0 offsetmin=0.03;
run;

The SGPLOT procedure overlays three separate STEP statements to display the ECDF, the 95% lower band, and the 95% upper band. I did not use the BAND statement (which creates a filled region) because it is designed for continuous curves. You could use a POLYGON statement if you need to display a filled region. I leave that as an exercise.

How to interpret the confidence limits

The K-S bands are simultaneous confidence bands, not pointwise confidence intervals. There is a big difference!

If these were pointwise intervals, a 95% interval at a specific value like x = 7 would mean, "the true CDF value at x = 7 is within the vertical range for 95% of random samples." However, the interpretation is different for simultaneous bands. The bands indicate that, in 95% of random samples, the entire CDF curve for the population will lie between these upper and lower step functions for all values of x. This is a stronger statement, which is why the Kolmogorov bands are generally wide for small samples. The Cord data set has 50 observations, so wide bands.

Incidentally, it is straightforward to wrap the computation in this section into a SAS IML function that computes the confidence bands for any ECDF. Furthermore, you can compute the bands for any confidence level, not just for the 95% level. The Appendix defines a SAS IML function that computes Kolmogorov confidence bands for an ECDF for any confidence level.

Advantages and disadvantages of Kolmogorov bands

Some advantages of Kolmogorov bands include:

• Nonparametric: The Kolmogorov bands do not assume any specific distribution for the data. (For example, it is not assumed to be normal.) Thus, the bands are very general.
• Easy to Compute: The K-S bands are easy to compute. You merely need to add and subtract a term that depends on the sample size and the critical value of the Kolmogorov D statistic.

Kolmogorov bands have a few statistical shortcomings:

• Constant Width: Before the limits are truncated at 0 and 1, the K-S bands have the same vertical width at the center of the data as they do in the tails. In reality, the variance of an ECDF approaches 0 at the tails (because the CDF must exactly reach 0 and 1). Consequently, K-S bands are overly conservative (too wide) in the tails.
• Continuous Assumption: The critical values of the standard K-S test assume the underlying population distribution is perfectly continuous. If your data contains many tied values (perhaps due to rounding or limitations of the instruments used to measure the data), the bands will be conservative, meaning the true coverage probability might be higher than the nominal 95%.

Summary

Plotting an ECDF is one way to summarize univariate data. Like any statistic, the ECDF has uncertainty due to random sampling. You can use the asymptotic distribution of the Kolmogorov-Smirnov statistic to calculate simultaneous 95% confidence bands for the ECDF. Thus, you can overlay confidence bands on an ECDF. The K-S bands are somewhat conservative at the extremes of the data, but they are a classic, nonparametric, and easily interpretable way to visualize uncertainty in an ECDF.

The ECDF is a valuable tool for theoretical purposes and for ECDF-based goodness-of-fit tests. However, when comparing the data to a parametric distribution, many practitioners prefer Q-Q plots. Q-Q plots are preferred because it is easier to judge whether points fall on a straight line than it is to see the difference between two S-shaped curves. You can use the inverse-CDF transformation to map the K-S bands from the ECDF plot onto a Q-Q plot. I will demonstrate this transformation in a future article.

Appendix: A SAS IML function that computes Kolmogorov confidence bands for an ECDF

It is straightforward to wrap the computation in this article into a SAS IML function. The following program defines a function named ECDF_KSCL. It computes the K-S confidence bands for any confidence level in the range [0.75, 0.99], in increments of 0.01. If you pass in a level such as 0.975, it will be rounded to the nearest 0.01. The critical values for the K-S statistic were computed by using an IML program similar to the one in the lecture notes by Charles Geyer at U. MN.

proc iml;
/* Compute the Kolmogorov-Smirnov confidence bands for a given ECDF.
   The CL parameter specifies the confidence level, which defaults to 0.95.
   The CL parameter should be in the range [0.75, 0.99]. */
start ECDF_KSCL(_ECDF, CL=0.95);
   ECDF = colvec(_ECDF);      /* ensure column vector */
   n = countn(ECDF);          /* count the nonmissing values */
   KS_band = j(nrow(ECDF), 2, .);
   /* Enter a pre-tabulated set of critical values for confidence levels */
   dAlpha = 0.01;
   alpha = do(0.01, 0.25, dAlpha);
   conf_level = round(1 - alpha, dAlpha);
   ks_crit = { 1.6276236 1.517427 1.4490862 1.3985734 1.3580986 
               1.3241093 1.2946745 1.2686236 1.2451911 1.2238479 
               1.2042125 1.1860005 1.1689938 1.1530214 1.1379465 
               1.1236583 1.110065 1.0970904 1.0846701 1.0727492 
               1.0612805 1.0502232 1.0395418 1.0292048 1.0191847 };
   j = loc( conf_level = round(CL, dAlpha) );
   if ncol(j)=0 then do;
      print "ERROR in ECDF_KSCL: The confidence level must be between 0.75 and 0.99";
      return( KS_band );
   end;
   /* Compute the standard error for the K-S bands */
   SE_KS = ks_crit[j] / sqrt(n);
 
   /* construct upper and lower 95% CI bands centered on the ECDF.
      The >< operator returns the minimum, ensuring max value is 1.
      The <> operator returns the maximum, ensuring min value is 0.
      See https://blogs.sas.com/content/iml/2026/02/04/clip-values.html */
   CL_Upper = (ECDF + SE_KS) >< 1;
   CL_Lower = (ECDF - SE_KS) <> 0;
   KS_band = CL_Lower || CL_Upper;
   return KS_Band;
finish;
store module=(ECDF_KSCL);
QUIT;
 
/* Test the ECDF_KSCL function */
proc iml;
load module=(ECDF ECDF_KSCL);
 
/* read the data for the ECDF */
use Cord;  read all var "Strength" into x;  close;
/* evaluate the ECDF at all data points */
call sort(x);
ecdf = ECDF(x);
CL = ECDF_KSCL(ecdf, 0.90); /* get bands for any confidence levels in [0.75, 0.99] */
CL_Lower = CL[,1];
CL_Upper = CL[,2];
create ECDF_bands var {"x" "ECDF" "CL_Lower" "CL_Upper" };
  append;
close;
QUIT;
 
/* Graph the step functions by using PROC SGPLOT */
title "Empirical CDF and 90% Confidence Bands";
proc sgplot data=ECDF_bands noautolegend;
   label x="Breaking Strength (psi)" ECDF="Cumulative Proportion";
   step x=X y=ECDF / lineattrs=(thickness=2);
   step x=X y=CL_Lower / lineattrs=(color=gray);
   step x=X y=CL_Upper / lineattrs=(color=gray);
   xaxis grid label="x";
   yaxis grid min=0 offsetmin=0.03;
run;

The post How to add confidence bands on the ECDF plot in SAS appeared first on The DO Loop.

In statistics, we assume that a sample is a random draw from an underlying distribution. Each random sample has a different ECDF. Therefore, you can think of an ECDF as a "point estimate" of the underlying cumulative distribution. If you were to draw a different sample of the same size, you would obtain a different ECDF function. Assuming that the data came from the same data-generating process (that is, the same probability distribution) an interesting question is: How much variability should you expect in a random ECDF?

I like to answer questions like this by using a simulation approach. A simulation builds intuition with a minimal amount of math. After understanding the main ideas and concepts, the math usually becomes easier to understand.

This article creates a simulation in SAS that enables you to visualize the variability for the ECDF of a random normal sample.

How much should you expect an ECDF to vary?

This section uses simulation to generate many normal samples. For each sample, you can compute the ECDF plot and display them in a single graph. This visualizes the concept of variability of the ECDF for a random sample. To see variability in histograms and Q-Q plots, you can read a previous article.

A previous article showed how to compute the ECDF of a real data set containing the breaking strength of fiber-optic cords. That data appears to be approximately normally distributed with mean μ = 7 and standard deviation σ = 0.13.

Let's simulate 1,000 random samples from N(7, 0.13) and compute the ECDF for each. Let's then plot the ECDFs and overlay two curves: the true normal CDF and the ECDF for the real data. To do this, you can use the ECDF function in SAS IML. The ECDF function is natively supported in SAS Viya, or you can load the SAS 9.4 function that I created in the last article:

/* Before running this program in SAS 9.4, define and STORE the ECDF function from 
   https://blogs.sas.com/content/iml/2026/05/26/create-ecdf.html */
proc iml;
load module=(ECDF);
/* One way to visualize the variability in the ECDF estimate is 
   to run a simulation.
   Repeat B times:
     Generate a random sample of size 50 from N(mu, sigma)
     Compute the ECDF
   Then plot the B ECDFs from the simulated samples and overlay 
   the ECDF for the data. The following carries out that plan.
*/
mu = 7;
sigma = 0.13;
NSamples = 1000;
NObs = 50;
X = j(NObs, NSamples, .);     /* each column of X is a random sample */
call randseed(1234);
call randgen(X, "Normal", mu, sigma);
 
ECDF = j(NObs, NSamples, .);
do i = 1 to NSamples;
  ECDF[,i] = ecdf( X[,i] );   /* the i_th column is an ECDF of X[,i] */
end;
 
/* the CREATE/APPEND statement outputs data row-wise, so transpose.
   Now the i_th row X[i,] is a sample and ECDF[i,] is the ECDF. */
X = X`;
ECDF = ECDF`;
create ECDF_Sim var {"X" "ECDF"};  append;  close;
 
/* create some overlays */
/* First, overlay the parent probability distribution */
t = do(mu - 3.5*sigma, mu + 3.5*sigma, sigma/100);
CDF = cdf("Normal", t, mu, sigma);
create CDF var {"t" "CDF"};  append;  close;
 
/* Next, read data, compute ECDF, write to data set */
use Cord;  read all var "Strength";  close;
call sort(Strength);
ecdf_data = ECDF(Strength);
create ECDF_Data var {"Strength" "ECDF_Data"};  append;  close;
QUIT;
 
/* merge the ECDF, the simulated ECDFs, and the CDF from N(mu, sigma) */ 
data ECDF_All;
set ECDF_Sim CDF ECDF_Data;
run;
 
/* overlay the 1,000 simulated ECDFs and the real ECDF */ 
title "1,000 Random Normal ECDFs";
proc sgplot data=ECDF_All noautolegend;
   label x="Breaking Strength (psi)" ECDF="Cumulative Proportion";
   scatter x=X y=ECDF / transparency=0.98;
   series x=t y=CDF / lineattrs=(color=black thickness=2);
   step x=Strength y=ECDF_Data / lineattrs=(color=DarkRed thickness=3);
   xaxis grid label="x";
   yaxis grid;
run;

Let's understand this graph. The blue point cloud contains points for 1,000 ECDFs from random samples of size N=50 from a normal distribution. (Technically, I should have overlaid 1,000 step functions, but plotting the points with high transparency creates a cleaner "density cloud" effect.) On top of the cloud, I overlaid a black curve, which is the true CDF of the population distribution. I also overlaid a dark red step function, which is the ECDF of the Strength variable in the Cord data set. This ECDF seems to fit inside the blue cloud, which indicates that you cannot rule out that the Strength data came from an N(7, 0.13) distribution.

You can use this visual cloud to answer two related questions about sampling variability:

• Given a specific value X, what is its likely percentile? You can estimate the answer by drawing a vertical line on the graph. The intersection of that line with the vertical spread of the blue cloud indicates the variability of the quantiles. For example, if you draw a vertical line at the mean X=7, the quantile range spans approximately [0.35, 0.65]. This indicates that the value 7 in any given random sample of size 50 could reasonably fall anywhere between the 35th and 65th percentiles.
• Given a specific percentile, what is the likely range of X values? You can estimate this by drawing a horizontal line. The intersection of the horizontal line and the horizontal width of the blue cloud indicates the variability of the X values for that exact quantile. For example, if you draw a line for the median (q=0.5), the intersection with the cloud appears to span X values roughly in the range [6.95, 7.05]. This indicates where the sample median is likely to fall.

Notice that the dark red ECDF from the Cord data fits almost entirely inside the densest part of the blue cloud. This suggests that it is reasonable to hypothesize that the Strength data came from an N(7, 0.13) distribution. You might want to run a formal statistical test for this hypothesis. You can use PROC UNIVARIATE to run several ECDF-based goodness-of-fit tests, as follows:

proc univariate data=Cord normal;
  cdfplot Strength / normal(mu=7 sigma=0.13);
  ods select TestsForNormality CDFPlot;
run;

The output includes the TestsForNormality table, which shows tests statistics and p-values for several tests for normality. All tests fail to reject the null hypothesis, which tells us that the data are consistent with N(7, 0.13).

Summary

Statistics on small random samples have large variability. Graphical plots of the data (such as histograms, Q-Q plots, and ECDF plots) will naturally vary from sample to sample. The simulation approach in this article visualizes a region where an ECDF is likely to be for a random sample. If your observed data falls within this region, it is plausible the data came from the hypothetical distribution. This intuitive idea motivates looking at adding confidence intervals to Q-Q plots and ECDF plots, which is a topic for a future article.

The post How to visualize the variability in an ECDF appeared first on The DO Loop.

This article shows how to construct and plot an ECDF in SAS. The easy way is to use PROC UNIVARIATE. However, for certain data analysis applications, it is useful to be able to call a function in the SAS IML matrix language that manually evaluates the ECDF. The ECDF function was added to the SAS IML language in the 2025.01 release of SAS Viya. For customers who are still on SAS 9.4, this article provides a similar ECDF function that you can use.

What is an ECDF?

Given a set of univariate data x₁ ≤ x₂ ≤ ... ≤ x_n, the empirical distribution function (ECDF) is the step function defined by:
      F(t) = (number of data values ≤ t) / n

The ECDF is a piecewise-constant function that is flat on each half-open interval [x_i, x_i+1). At each unique data value, it increases by an amount proportional to the frequency of data at that value.

The ECDF in PROC UNIVARIATE

SAS provides the CDFPLOT statement in PROC UNIVARIATE, which computes and displays the ECDF and optionally overlays the CDF of a specified probability distribution. For example, the PROC UNIVARIATE documentation includes data about the breaking strength (in PSI) for 50 fiber-optic cords:

data Cord;
   label Strength="Breaking Strength (psi)";
   input Strength @@;
datalines;
6.94 6.97 7.11 6.95 7.12 6.70 7.13 7.34 6.90 6.83
7.06 6.89 7.28 6.93 7.05 7.00 7.04 7.21 7.08 7.01
7.05 7.11 7.03 6.98 7.04 7.08 6.87 6.81 7.11 6.74
6.95 7.05 6.98 6.94 7.06 7.12 7.19 7.12 7.01 6.84
6.91 6.89 7.23 6.98 6.93 6.83 6.99 7.00 6.97 7.01
;
 
title 'Cumulative Distribution Function of Breaking Strength';
proc univariate data=Cord noprint;
   cdfplot Strength / odstitle = title;
run;

From the plot, you can see that there is about a 50% chance that a random fiber-optic cord will break if subjected to 7.0 PSI of stress. If the stress is kept below 6.9 PSI, there is only a 20% chance of random cord will break. The graph can be improved to answer questions like these by using the STATREF= option on the CDFPLOT statement to add vertical reference lines at the data value that correspond to arbitrary percentiles. For example, the option statref=P 10 Q1 Q2 Q3 P 90 add vertical lines at the 10th, 25th, 50th, 75th, and 90th percentiles of the data.

In addition to displaying the ECDF, the CDFPLOT statement enables you to overlay the theoretical CDF of common probability distributions. This is useful for visualizing ECDF-based hypothesis tests. For example, the Kolmogorov-Smirnov (K-S) test for normality uses the largest vertical gap between the ECDF and the hypothesized CDF to construct a test statistic. If you add the NORMAL option to the CDFPLOT statement, you can visualize the K-S test and other ECDF-based goodness-of-fit tests for normality.

How to construct an ECDF manually

When conducting a simulation study, a bootstrap analysis, or exploring ECDF-based statistics, it is convenient to compute ECDF values in SAS IML. To construct the ECDF manually, do the following:

• Find the unique values in the data and the number of observations for each unique value. In SAS IML, you can use the TABULATE function for that.
• Since the ECDF is a piecewise constant function that steps up at these locations, use these unique values as the endpoints of bins. You can use the BIN function for that.
• For each bin, associate the value of the ECDF function at the left endpoint of the bin. Since the ECDF is a piecewise constant function that increases at these locations, this tells you the values of the ECDF for any point in any bin.
• For points that are outside of any bin, the ECDF is either 0 or 1. If t is less than the smallest data value, ECDF(t) = 0. If t is greater than or equal to the largest data value, then ECDF(t) = 1.

In SAS Viya, the ECDF function is supported natively by the SAS IML language. For customers who are still on SAS 9.4, the following function manually implements the basic ECDF algorithm:

proc iml;
/* Compute the empirical distribution function (ECDF) if a column vector of data.
   Optionally, evaluate the ECDF at values in a column vector, t. 
   If t is not specified, use t=x. 
   SYNTAX:
     y = ECDF(x); 
     t = T( do( min(x), max(x), (max(x)-min(x))/100 ) );
     y = ECDF(x, t); 
*/
start ECDF(_x, _t=);
    x = colvec(_x);      /* ensure a column vector */
    /* Tabulate the unique values of x.
       levels = unique sorted elements of x
       freq = number of duplicates for each level */
    call tabulate(levels, freq, x);
    levels = colvec(levels);
    freq = colvec(freq);
    y = cusum(freq) / sum(freq);
 
    if isSkipped(_t) then
       t = x;
    else
       t = colvec(_t);      /* ensure column vector */
    /* build a step function by assigning t to different bins, using levels as endpoints */
    bins = bin(t, levels);
 
    /* assign ECDF[i] for nonmissing bins */
    ecdf = j(nrow(t), 1, .);    /* initialize ECDF to missing */
    idx = loc(bins ^= .);
    if ncol(idx) > 0 then
        ecdf[idx] = y[ bins[idx] ];
 
    /* Three reasons bins[i] could be a missing value:
       1. if t[i] is missing. No need to change ECDF[i] b/c it was initialized to missing.
       2. if t[i] < min(x). Set ECDF(t[i]) = 0.
       3. if t[i] > max(x). Set ECDF(t[i]) = 1.
    */
    LT_idx = loc(t^=. & t<min(x));
    if ncol(LT_idx) > 0 then 
       ecdf[LT_idx] = 0;
    GT_idx = loc(t^=. & t>=max(x));
    if ncol(GT_idx) > 0 then 
       ecdf[GT_idx] = 1;
    return ecdf;
finish;
store module=(ECDF);
QUIT;

Let's call this function on the same data that was used by PROC UNIVARIATE earlier in this article. Suppose you want to estimate the probability that the breaking strength of a fiber optic cord is less than certain specified test values. You can read the Cord data set and call the ECDF function, as follows:

proc iml;
load module=(ECDF);
 
use Cord;
read all var "Strength" into x;
close;
 
/* estimate the probability that the breaking strength is less than 
   6.8, 6.9, 7.0, and 7.1 psi. Note that some of these values 
   are not observed data values.
*/
t = {6.8, 6.9, 7.0, 7.1};
ecdf_t = ECDF(x, t);
print t ecdf_t;

The output estimates the probability that a random cord will break when subjected to the specified stresses. About 4% will break if the stress is 6.8 PSI, whereas about 76% will break is the stress is as large as 7.1 PSI.

Notice that some of the specified locations are not data values. For each specified value, the ECDF returns an estimate of the probability that a random fiber-optic cord will break when subjected to a stress of that magnitude.

Graph the ECDF

Typically, an ECDF plot evaluates the ECDF function at the observed data values rather than arbitrary points. As discussed in a previous article, it is best to use the STEP statement in PROC SGPLOT to visualize an ECDF function.

The following IML statements sort the data, evaluate the ECDF at the data locations, and write the resulting calculations to a SAS data set. The SGPLOT procedure then visualizes the ECDF.

/* typically, an ECDF plot evaluates the ECDF at all data points */
call sort(x);
ecdf = ECDF(x);
 
create ECDF var {"x" "ECDF"}; append; close;
QUIT;
 
/* Graph a step function. See
   https://blogs.sas.com/content/iml/2016/09/06/graph-step-function-sas.html */
title "Empirical CDF";
proc sgplot data=ECDF noautolegend;
   label x="Breaking Strength (psi)" ECDF="Cumulative Proportion";
   step x=x y=ecdf;
   fringe x;
   xaxis grid label="x" offsetmin=0.05 offsetmax=0.05;
   yaxis grid min=0 offsetmin=0.03;
run;

The resulting graph is similar to the one produced by PROC UNIVARIATE earlier. However, by using PROC SGPLOT, you have more freedom to enhance the plot with features such as a fringe plot, custom reference lines, grid lines, and so on. The grid lines enable you to quickly estimate that about 50% of random cords will creak if subjected to a stress of 7.0 PSI.

Summary

The empirical CDF (ECDF) is a fundamental nonparametric tool for exploring data and comparing univariate distributions. In Base SAS, PROC UNIVARIATE provides a quick and easy way to generate these plots automatically. In addition, it can be useful to be able to evaluate the ECDF function programmatically. The ECDF function was added to SAS IML in SAS Viya in Release 2025.01. For SAS 9.4 customers, this article provides a similar function.

The post Create and graph the ECDF function in SAS appeared first on The DO Loop.

You might hope that higher-dimensional Childs-Sun formulas are similarly simple. Alas, they are not! For the higher-dimensional formulas, the probabilities are defined recursively. As the dimension increases, the number of terms explodes combinatorially. This article shows a SAS IML module that implements the Childs-Sun formula for the orthant probability of a 4-D standardized normal distribution. After implementing the 4-D case, I have no desire to attempt higher dimensions!

The Childs-Sun formula for dimension 4

Appendix 1 shows a snippet from p. 12 of Genz and Bretz (2009), which describes the Childs-Sun formula for an MVN orthant probability. The formula relies on the Childs-Sun reduction method (Childs, 1967; Sun, 1988). Similar to the lower-dimensional formulas, the 4-D formula contains a sum of arcsines of the correlation coefficients. However, the formula adds additional terms that involve an integral, I₄. The integrands are complicated expressions of the off-diagonal elements of a 4x4 correlation matrix.

In 4-D, there are 16 orthants. As discussed in a previous article, if you can find the probability in the "first orthant" (the one where all components have positive signs), you can use a clever "reflection trick" to find the probability in the other 15 octants.

In Genz and Bretz (2009), the Childs-Sun formula for the orthant probability, which they call P₄, has an unfortunate typo. The correct formula is shown in Appendix 1, where the red "delete" symbol removes an erroneous "plus sign." It took me many hours to decipher the notation. By consulting other references (Gehrlein and Saniga (1977); Sun and Asano (1990)), and with some assistance from the Google Gemini AI, I was able to decipher the formulas and create a SAS IML function that evaluates the orthant probability. The IML function is presented in this article.

Conditional probability and the geometry of the formula

I usually explain the mathematical details of my SAS programs, but I'm not going to do that today. The formula is too complicated to describe in a short blog post. Some details are in the original Childs (1967) paper, which is unfortunately behind a paywall.

The general formula (see Appendix 1) is very complicated. But if you restrict to k=4 dimensions (hence, n=2), the formula looks something like this:
$P_4 = \frac{1}{16} + \frac{1}{8\pi} \sum_{i < j}^4 \sin^{-1}(R_{ij}) + \frac{1}{4\pi^2} I_4(R)$
where I₄ is a complicated sum of one-variable integrals. Each integral depends on a lower-dimensional integration, called I₂, which is why this is a recursive formula. The notation is explained further in Appendix 2.

The one-variable integrals are the key. In general, to evaluate the probability of a 4-D distribution, you must solve a 4-D volume integral. The Childs-Sun method uses the geometry of the multivariate normal distribution to reduce this volume to a sum of iterated integrals that accumulate lower-dimensional slices.

This formula is an application of conditional probability. The dummy integration variable, x, is a continuous variable on [0,1]. For a fixed value of x, the term R¹ⁱ(x) evaluates to a conditional correlation matrix. The variables inside the integrand (e.g., the quadratic denominator w = 1 - R_1i² x²) are conditional variances and partial correlations.

The I₂ function evaluates a 2-D orthant probability (which has an exact arcsine solution) for the remaining two variables, given that the other variables are fixed at a slice determined by x. You can use the QUAD subroutine in SAS IML to integrate over x in [0,1], which aggregates these 2-D conditional probabilities to sweep out the exact probability mass of the 4-D volume.

Implementing the P4 formula in SAS IML

The main SAS IML function is called P4_childs, which computes the 4-D orthant probability. It calls the QUAD subroutine to compute I₄. The integrand is ProbInt4, which builds terms from the off-diagonal elements of the 4x4 correlation matrix. The DO loop evaluates three different conditional correlation matrices.

proc iml;
/* The base case: P2. For a scalar correlation, rho, returning the orthant probability, P2.
   See https://blogs.sas.com/content/iml/2026/05/11/mvn-orthant-probability.html */
start P2(rho);
  return( 0.25 + arsin(rho) / (2*constant("pi")) );
finish;
 
/* The sum of the I4 integrands. Instead of evaluating the three integrals separately,
   evaluate the sum of the integrals and perform one numerical integration.   
   Evaluates the Sun (1988) residual partial correlation matrix.
   Whereas Genz/Bretz Eqn 2.8-2.9 uses a general covariance matrix, this function assumes 
   R is a correlation matrix, so R[i,i]=1  */
start ProbInt4(x) global(R_global);
  R = R_global;
  pi = constant("pi");
  x2 = x # x;
  sum = 0;
 
  /* The sum of three terms: i=2, 3, and 4 */
  do i = 2 to 4;
    /* Determine the indices (j, k) of the two remaining variables */
    if i=2 then do;      j=3; k=4; end;
    else if i=3 then do; j=2; k=4; end;
    else do;             j=2; k=3; end;
 
    /* R_ii term: The quadratic denominator */
    w = 1 - R[1,i]##2 * x2;          /* R[i,i]=1 for all i */
    if w< 1E-12 then w = 1E-12;      /* Numerical safeguard */
 
    /* Compute the conditional variances and covariance. These are a_11, a_12, etc, in Sun and Asano */
    a_jj = 1      -R[1,j]##2 * x2   -((R[i,j] - R[1,i]*R[1,j]*x2)##2) / w;
    a_kk = 1      -R[1,k]##2 * x2   -((R[i,k] - R[1,i]*R[1,k]*x2)##2) / w;
    a_jk = R[j,k] -R[1,j]*R[1,k]*x2 -(R[i,j] - R[1,i]*R[1,j]*x2)*(R[i,k] - R[1,i]*R[1,k]*x2) / w;
 
    /* rho is the off-diagonal element of the conditional correlation matrix */
    rho = a_jk / sqrt(a_jj * a_kk);
 
    /* ensure rho is in [-1,1]. See
       https://blogs.sas.com/content/iml/2026/02/04/clip-values.html */
    rho = ( -1 <> (rho >< 1) ); /* clamp to [-1,1] */
 
    /* Call P2 function to get the I2 term */
    I2_val = 2 * pi * ( P2(rho) - 0.25 );     /* this actually simplifies to ARSIN(rho) :-) */
 
    /* Add the i-th term to the sum */
    sum = sum + (R[1,i] / sqrt(w)) * I2_val;
  end;
  return(sum);
finish;
 
/* The main P4 function, which evaluates the Childs-Sun formula for 4 dimensions */
start P4_Childs(R) global(R_global);
  R_global = R; /* Set the global correlation matrix for QUAD */
 
  /* Integrate ProbInt4 over the interval [0, 1] */
  call quad(I4, "ProbInt4", {0 1});
 
  /* Apply the Childs (1967) formula */
  pi = constant("pi");
  sum_arsin = sum( arsin(R[{2 3 4 7 8 12}]) );  /* sum of arcsine for all 6 off-diagonal elements */
  prob = 1/16 + sum_arsin / (8*pi) + I4 / (4 * pi**2);  /* fixes the typo in Genz and Bretz, Eqn 2.8 */
  return(prob);
finish;
 
/* --- Simple Test --- */
R = {1.0  0.5  0.3  0.2,
     0.5  1.0  0.4  0.3,
     0.3  0.4  1.0  0.5,
     0.2  0.3  0.5  1.0};
 
ans = p4_Childs(R);
print ans;

When you run this code, the exact Childs-Sun formula produces an orthant probability of 0.1611218 for the given matrix. You can verify that value by using a Monte Carlo simulation. The following Monte Carlo simulation uses one million random points from the MVN(0, R) distribution. In addition to the Monte Carlo estimate, I generate a 95% confidence interval for the probability.

/* Use Monte Carlo simulation to estimate the probability that a 
   MVN random variable is less than a specified value in each coordinate.
   Sigma is a kxk covariance matrix; mu is an option row vector with k elements.
   The row vector b specifies the upper limits of integration. 
   The function returns an estimate of 
   P(X1<b[1] & X2<b[2] & ... & Xk<b[k] | X~MVN(mu, Sigma))
   by simulating N random variates from MVN(mu, Sigma) and returning the proportion
   that are in the specified region.
*/
start MC_CDFMVN(N, b, Sigma, mu=j(1,ncol(Sigma),0));
   X = randnormal(N, mu, Sigma);
   inRegion = j(N,1,1);
   do i = 1 to ncol(b);
      v = (X[,i] < b[i]);
      inRegion = inRegion & v;
   end;
   return mean(inRegion);
finish;
 
N = 1E6;
prob_MC = MC_CDFMVN(N, {0 0 0 0}, R);
SE_MC = sqrt( prob_MC * (1-prob_MC)/N );
Lower95 = prob_MC - 1.96*SE_MC;
Upper95 = prob_MC + 1.96*SE_MC;
print prob_MC Lower95 Upper95;

The output from the Monte Carlo simulation shows that the Monte Carlo estimate is close to the value returned by the Childs-Sun formula. In addition, the 95% confidence interval includes the exact value, which gives me "confidence" (pun intended!) that I implemented the Childs-Sun formula correctly. You can repeat the computations for other 4x4 correlation matrices and obtain similar results.

Summary

As my dear mother used to say when faced with a challenging situation, "this is not for the faint of heart." The Childs-Sun formula for the orthant probability of a multivariate normal distribution is a complicated equation. In 2-D and 3-D, it simplifies into a sum of arcsines. In higher dimensions, you must compute conditional correlation matrices and integrate them to evaluate the formula. I have shown a SAS IML function that implements the formula in 4-D.

Despite the complexity, this is an amazing mathematical and computational result! The Childs-Sun formula is an exact analytical solution to a 4-dimensional multivariate normal volume integral over an orthant. And, as shown previously, you can use the P4_Childs function to compute the MVN probability for any 4x4 correlation matrix in any orthant.

Appendix 1: The general Childs-Sun formula for orthant probabilities

The following image is from p. 12 of Genz and Bretz (2009), Computation of Multivariate Normal and t Probabilities. The last term of Eqn 2.8 has a typo in Genz and Bretz. The equation contains a plus sign that should not be there. This typo caused me hours of confusion! I finally cross-referenced the equation with Sun and Asano (1990) and the original article by Childs (1967), which both display the correct equation.

Note the somewhat confusing notation in the Childs-Sun formula. The superscript on the argument to the I_2j function describes the submatrix to use in computing the conditional probabilities. This is not explained further in G&B. For the details of the implementation, see the SAS IML function in this article. The superscript on the argument to the I₄ function describes the submatrix to use in computing the conditional probabilities. This is not explained further in G&B. For the details of the implementation, see the SAS IML function in this article.

Appendix 2: The sum notation

The Genz and Bretz notation for the sum of integrals is very dense. Here is what it means in the 4-D case. For this case, n=2, which means that the outer sum reduces to j=2. The inner sum looks like
$\sum_{i_1 < \ldots < i_4}^4 I_4(R^{i_1,\ldots,i_4})$.

The notation is a combination sum. It means, "Sum over all possible combinations of the distinct indices chosen from the set {1, 2, 3, 4} subject to the constraint that the combinations are in increasing order." Well, for four indices, there in only one combination that satisfies this constraint, and it is the tuple (1, 2, 3, 4). So the inner summation also vanishes and becomes I₄(R). You can drop the superscript because you are sending in all four columns of the 4x4 matrix, R.

The I₄ term is defined recursively in terms of I₂. Now the sum matters. When j=1, the summation is
$\sum_{i_1 < i_2}^4 I_2(R^{i_1,i_2})$. There are six ways to choose two distinct ordered indices from the set {1,2,3,4}. Since I₂(R^{i1, i2}) = arcsine(R[i1, i2]), we get the sum-of-arcsine terms $\sum_{i < j} \sin^{-1}( R_{ij} )$.

The post Compute the orthant probability for the 4-D multivariate normal distribution appeared first on The DO Loop.

The probability depends on the density (PDF) of the distribution as well as the region of integration. One special region is the orthant. An orthant is a generalization of a quadrant. In two dimensions, the Cartesian coordinate system divides the plane into four quadrants. In three dimensions, the Cartesian coordinate system divides the plane into eight octants. In k dimensions, there are 2^k orthants.

It is rare to find an exact formula for the probability of a multivariate distribution. Amazingly, there are formulas for the probability of a standardized multivariate normal (MVN) distribution in an orthant (Genz and Bretz (2009) p. 11-13, who cite Sheppard (1898), D. R. Childs (1967), and H. J. Sun (1988)).

This article shows the formulas for dimensions 1, 2, and 3. In a subsequent article, I will discuss a SAS IML function that implements a formula for dimension k=4. The general case was solved by Childs (1967) and extended by Sun (1988), so I refer to these formulas as the Childs-Sun formulas.

As shown in a previous article, it suffices to compute a probability for the standardized MVN distribution because a probability for a general MVN distribution can always be found by standardizing the problem. To save typing, this article deals exclusively with the standardized MVN distribution and often omits the word "standardized".

Orthant probabilities

For a vector X, the signs of X are constant within each orthant. For example, a 2-D vector has signs (+1, +1) when it is in the first quadrant, and signs (-1, -1) when it is in the third quadrant. The term "orthant probability" traditionally refers to the first orthant, where all signs are positive.

The MVN distribution is symmetric under the transformation X → -X, so the probability that X~MVN(0,R) is in the "all-signs-positive" octant equals the probability that X is in the "all-signs-negative" octant. The "all-signs-negative" octant is important in probability and statistics. It represents the "CDF region," which is a left-tail probability. For example, for a 2-D random vector X=(X1,X2), the CDF evaluated at (0,0) is Φ₂(0,0) = P(X1<0 and X2<0), which is the probability that X is in the third quadrant. In terms of an integral,
$\Phi_2(0,0; \rho) = \int_{-\infty}^{0} \int_{-\infty}^{0} \phi_2(x_1, x_2; \rho)\, dx_2 \, dx_1$
Here, ρ is the Pearson correlation between X1 and X2, φ₂ is the bivariate standard normal density, and Φ₂ is the bivariate standard normal CDF.

This concept generalizes to higher-dimensional distributions: the "CDF orthant" is the "all-signs-negative" orthant for which every coordinate of X has a negative sign. The CDF of the standardized normal distribution at the origin is the integral of the distribution's density over that orthant. This article shows some Childs-Sun formulas for these probabilities for a general correlation matrix.

Orthant probabilities in dimensions 1 and 2

For the one-dimensional problem, the 1-D orthant is the set {X < 0}. The CDF at zero is P₁ = Φ₁(0) = 0.5. You can obtain this value by solving an integration problem, but it is also evident from the fact that the normal distribution is symmetric about 0. For this article, the important fact is that there is a formula for the 1-D orthant probability: P₁ = 1/2.

For 2-D, the correlation matrix has one parameter, ρ = R[1,2], which is the correlation between the bivariate normal variates X1 and X2. The 2-D (bivariate) formula (originally found by Sheppard, 1898) is
$P_2 = \Phi_2(0,0; \rho) = \frac{1}{4} + \frac{1}{2\pi} \sin^{-1}(\rho)$

Let's run a simple DATA step to verify this result. In SAS, you can use the PROBBNRM function to compute a bivariate normal probability:

/* 2-D orthant probability for various rho values */
data MVN2D;
input rho @@;
pi = constant('pi');
P2_formula = 1/4 + 1/(2*pi) * arsin(rho);
P2_integral = probbnrm(0, 0, rho);
diff = P2_formula - P2_integral;
drop pi;
datalines;
-0.8 -0.5 0 0.5 0.8 
;
 
proc print data=MVN2D noobs; 
run;

As shown in the output, the formula agrees with the PROBBNRM function to within machine precision.

Orthant probabilities in dimension 3

In 3-D the standard multivariate normal distribution has three parameters: the off-diagonal elements of the 3x3 correlation matrix, R. The formula for the orthant probability in dimension 3 is given by
$P_3 = \frac{1}{8} + \frac{1}{4\pi} \left( \sin^{-1}(R_{12}) + \sin^{-1}(R_{13}) + \sin^{-1}(R_{23}) \right)$

Notice that you can write P3 by using a summation formula in which you sum over the indices of the upper triangular elements of the correlation matrix, R.

To validate this formula, let's use a Monte Carlo estimation routine, which was developed in a previous article. The following SAS IML program evaluates the formula for a few correlation matrices, and compares the formula and the Monte Carlo estimate.

/* 3-D orthant probability for various correlation matrices, R */
proc iml;
/* Use Monte Carlo simulation to estimate the probability that a 
   MVN random variable is less than a specified value in each coordinate.
   See https://blogs.sas.com/content/iml/2026/05/04/standardize-mvn-probabilities.html
   The function estimates CDF(b).
*/
start MC_CDFMVN(N, b, Sigma, mu=j(1,ncol(Sigma),0));
   X = randnormal(N, mu, Sigma);
   inRegion = j(N, 1, 1);
   do i = 1 to ncol(b);
      v = (X[,i] < b[i]);       /* binary variable for the i_th coordinate */
      inRegion = inRegion & v;  /* find the intersection of all conditions */
   end;
   return mean(inRegion);       /* proportion of random points that satisfy criterion */
finish;
 
/* Evaluate the MVN CDF at the origin (="orthant probability") by using the exact P3 formula */
start Orthant_P3(R);
  pi = constant("pi");
  upper_tri_corrs = R[ {2 3 6} ];
  sum_arsin = sum( arsin( upper_tri_corrs ) );
  prob = 1/8 + sum_arsin / (4 * pi);
  return(prob);
finish;
 
N = 1E6;     /* size of Monte Carlo sample */
b = {0 0 0}; /* the origin */
call randseed(1234);
 
/* Example 1: AR(1) correlation matrix */
R = {1    0.5 0.25,
     0.5  1   0.5,
     0.25 0.5 1    };
P3_formula = Orthant_P3(R);
P3_MC_est  = MC_CDFMVN(N, b, R);
diff       = P3_formula - P3_MC_est;
print P3_formula P3_MC_est diff;
 
/* Example 2: positive and negative correlations */
R = {1   -0.5 -0.1,
    -0.5  1    0.8,
    -0.1  0.8  1    };
P3_formula = Orthant_P3(R);
P3_MC_est  = MC_CDFMVN(N, b, R);
diff       = P3_formula - P3_MC_est;
print P3_formula P3_MC_est diff;

Results are shown for two different correlation matrices. The formula in the Orthant_P3 function is exact. The Monte Carlo estimate is close to the exact value.

The computation for other 2-D quadrants

The formula is more powerful than it initially seems because of reflective symmetries of the multivariate normal distribution. The Childs-Sun formula enables you to compute the probability in all orthants, not just the "CDF orthant" where every component of the random variable has a negative sign. This is because you can change the signs of the components of the MVN random variable to obtain the other quadrants. Geometrically, "changing the sign" means considering the transformation X_i → -X_i for some coordinate. This transformation changes the sign of the correlation between X_i and the other variables. Thus, you can obtain the correlation matrix for any transformation that reflects one or more coordinates.

Let's see how this works in two dimensions. Let ρ be the off-diagonal (R[1,2]) correlation. Here is how to get the probability in all four quadrants:

1. Quadrant I: This is the orthant probability, P2 = 1/4 + arsin(ρ) / (2π).
2. Quadrant II: To find the probability in the region Q2 = {X1 < 0, X2 > 0}, define Y1 = -X1 and Y2 = X2. Y is bivariate normal with correlation -ρ. Because arsin(-ρ) = -arsin(ρ), the probability in Q2 is 1/4 - arsin(ρ) / (2π).
3. Quadrant III: To find the probability in the region Q3 = {X1 < 0, X2 < 0}, define Y1 = -X1 and Y2 = -X2. Y is bivariate normal with correlation ρ, so the probability over Q3 is P2.
4. Quadrant IV: To find the probability in the region Q4 = {X1 > 0, X2 < 0}, define Y1 = X1 and Y2 = -X2. Y is bivariate normal with correlation -ρ, therefore the probability in Q4 is 1/4 - arsin(ρ) / (2π).

The computation for other k-D orthants

The same idea applies for any dimension. The Childs-Sun formula gives the orthant probability in the first quadrant where all coordinates are positive. Every orthant is identified by a sign vector s=(s1, s2, ..., sk), where each component is either +1 or -1. To compute the probability that X~MVN(0, R) is in the orthant with the sign vector s, do the following:

1. Define a new MVN variable, Y = (s1*X1, s2*X2, ..., sk*Xk).
2. Compute the correlation matrix for Y. This requires changing some of the signs of the off-diagonal elements in R. See the Appendix for details.
3. Use the Childs-Sun formula to find the probability that Y is in the first orthant. This is also the probability that X is in the orthant with sign vector s.

Summary

In summary, this article shows three Childs-Sunformulas for computing the orthant probability for a standardized multivariate normal distribution with correlation matrix R in dimension k:

• If k=1, the univariate probability is $P_1 = \frac{1}{2}$.
• If k=2, the bivariate probability is $P_2 = \frac{1}{4} + \frac{1}{2\pi} \sin^{-1}(R_{12})$.
• If k=3, the trivariate probability is $P_3 = \frac{1}{8} + \frac{1}{4\pi} \left( \sin^{-1}(R_{12}) + \sin^{-1}(R_{13}) + \sin^{-1}(R_{23}) \right)$

You can apply these formulas as written to compute the probability in the first orthant (all signs positive) or in the "CDF orthant" (all signs negative). You can use the reflective symmetry of the MVN distribution to find the probabilities in every other orthant.

Appendix: How the correlation matrix changes under coordinate reflections

This main article indicates that you can apply reflection transformations for coordinates of X~MVN(0,R) and find a new correlation matrix, C, such that Childs-Sun formula applies to C gives the probability that X is in some orthant. This appendix shows the linear algebra steps.

Suppose X~MVN(0,R), and you want to find the probability that X is in an orthant defined by a sign vector s=(s1, s2, ..., sk). In this vector, each s_i is either +1 or -1. We can use s to define a new random vector, Y, such that Y is in the first orthant exactly when X is in the target orthant.

Let S be the k x k diagonal matrix whose diagonal elements are the components of the sign vector, s. Define the linear transformation Y = S X. Because S is a diagonal matrix of signs, the i_th component of Y is simply Y_i = s_i * X_i. By definition, if X is in the target orthant, then the sign of X_i matches s_i. Therefore, the product s_i * X_i is always positive, which means Y is in the first orthant.

Because Y is a linear transformation of MVN variables, it is also multivariate normal. The expected value of Y is E[Y] = S E[X] = 0. The covariance matrix for Y (call it C) is given by the standard formula for the variance of a linear transformation:
$C = \text{Cov}(Y) = S \text{Cov}(X) S^T = S R S$

If you multiply the matrices, you will find that the elements of C are given by C_ij = s_i s_j R_ij. Because s_i² = 1, the diagonal elements of C are 1, which means C is a valid correlation matrix. The off-diagonal elements of C simply flip the sign of the correlation whenever exactly one of the two corresponding variables is reflected.

In summary, Y~MVN(0, C). The Childs-Sun formula with the correlation matrix C gives the probability that Y is in the first orthant. This is exactly the probability that X is in the orthant with the given sign vector, s.

The post Compute multivariate normal orthant probabilities appeared first on The DO Loop.

The concept of a "Z-score" also applies to multivariate normal distributions. In multidimensional statistics, you can standardize a multivariate normal (MVN) probability calculation by transforming the covariance matrix into the corresponding correlation matrix. This article shows how to shift and scale high-dimensional MVN distributions, both mathematically and by using the SAS IML language.

The 1-D Z-score transformation

Let's start by reviewing the Z-score transformation for the 1-D normal distribution. Let X be a normally distributed random variable with mean μ and standard deviation σ. In notation, X ∼ N(μ, σ). (Some authors use N(μ, σ²), but the SAS probability functions (CDF, PDF, QUANTILE, etc.) use the standard deviation as a parameter.) If you want to compute a left-tail probability, such as P(X < b), you can compute the probability on the original (data) scale, or you can apply a transformation and compute the probability on the standardized scale. This is useful in numerical computations because it reduces the problem to two steps: first, standardize the problem, then compute a probability on the standardized problem. Each of these steps can be developed, debugged, and tested independently.

Let's apply the Z-score transformation. Define a new random variable Z = (X - μ)/σ. Because a linear transformation of a normal random variable is also normal, Z follows a standard normal distribution: Z ∼ N(0, 1). By applying the same transformation to the integration limit b, we define a standardized limit c = (b - μ)/σ. Consequently, computing the probability on the original data scale is mathematically equivalent to computing the probability on the standardized scale:
      P(X < b | X~N(μ, σ)) = P(Z < c | Z~N(0, 1))

In SAS, you can compute these probabilities directly using the CDF function. The following SAS IML program demonstrates this equivalence. The goal is to find the probability that a random variable from N(μ, σ) is less than x. You can compute this directly on the original scale, or you can convert the limit to a Z-score and evaluate it using the standard normal distribution:

data Prob;
input mu sigma @@;
do x = -1 to 1;
   /* 1. Compute probability on original scale */
   prob_x = CDF("Normal", x, mu, sigma);
   /* 2. Compute probability on standardized scale */
   z = (x - mu)/sigma;
   prob_z = CDF("Normal", z);  /* N(0,1) */
   output;
end;
datalines;
 2   3
 1.5 2
-0.5   0.5
run;
 
proc print noobs; run;

The program uses several examples of x, μ and σ. For all examples, the PROB_X and the PROB_Z columns are identical. Geometrically, the X variable represents a data value whereas the Z variable measures the number of standard deviations from the mean. As shown in the next section, you can generalize this concept to multivariate normal distributions.

Standardizing the multivariate normal distribution

Suppose you want to compute the joint cumulative probability for a multivariate normal distribution. Let X be a k-dimensional random variable such that X ∼ MVN(μ, Σ), where μ = (μ₁, μ₂, ..., μ_k) is the mean vector and Σ is a k x k covariance matrix. The CDF or left-tail probability is P(X < b), where b = (b₁, b₂, ..., b_k) is a vector of upper limits, and the inequality is evaluated component by component.

To standardize this problem, you can standardize each original variable. Instead of doing this one variable at a time, you can use matrix-vector notation and transform all variables simultaneously. Let D be the diagonal matrix of standard deviations. That is, the diagonal elements of D are the square roots of the diagonal elements of Σ. Then you can express the covariance matrix as Σ = D*R*D, where the matrix R has 1s on the diagonal. The matrix R is the correlation matrix for the original variables (and the covariance matrix for the new standardized variables).

The standardized random variables are defined by Z = D^-1(X - μ). Z is multivariate normal. To find its parameters, you can compute its expected value and variance. The expected value is simply a vector of zeros. The variance of Z is R. Thus, the transformed variable Z follows a standardized multivariate normal distribution with mean zero and correlation matrix R. In notation, Z ~ MVN(0, R).

To complete the probability transformation, apply the same linear transformation to the integration limits, b. The standardized limits are c = D^-1(b - μ). A probability for X on the data scale (covariance scale) is equal to the probability for Z on the correlation scale:
      P(X < b | X ~ MVN(μ, Σ)) = P(Z < c | Z ~ MVN(0, R))

This is why, for example, the PROBBNRM function in SAS computes probabilities for the standardized variable Z~BVN(0, R), where BVN is the bivariate normal distribution, and R is a 2x2 correlation matrix. The developer who implemented the function assumed that users can transform any problem to the standardized scale.

Converting from covariance to correlation in SAS IML

The SAS IML language supports the COV2CORR function, which returns the correlation matrix that is associated with a covariance matrix. By using the math in the previous section, you can write a SAS IML function that will transform the limits of integration. The following IML function computes c = D^-1(b - μ). The SAS IML language enables you to perform this operation efficiently by using vectorized division, rather than multiplying by the inverse of a diagonal matrix:

proc iml;
/* Sigma is a kxk covariance matrix and b is a row vector with k elements.
   Convert b to inv(D)*(b-mu), where D is the diagonal matrix of
   standard deviations: sqrt(vecdiag(Sigma)).
   The mu vector is optional and defaults to the zero vector.
*/
start Xform_Limits_Cov2Corr( b, Sigma, mu=j(1,ncol(Sigma),0) );
   D = sqrt(vecdiag(Sigma));      /* extract diagonal (variances) for X */
   c = (b - mu) / rowvec(D);      /* new limits for Z */
   return c;
finish;

As we did for the 1-D case, let's show the equivalency between the CDF on the data scale and the CDF for the standardized problem. The exact computation of multivariate normal probabilities in high dimensions is difficult. Instead, I will use a simple Monte Carlo simulation to estimate the probabilities. This technique is easy to code and will suffice to demonstrate the equivalency of the probabilities.

/* Use Monte Carlo simulation to estimate the probability that a 
   MVN random variable is less than a specified value in each coordinate.
   Sigma is a kxk covariance matrix; mu is an optional row vector with k elements.
   The row vector b specifies the upper limits of integration. 
   The function returns an estimate of 
   P(X1<b[1] & X2<b[2] & ... & Xk<b[k] | X~MVN(mu, Sigma))
   by simulating N random variates from MVN(mu, Sigma) and returning the proportion
   that are in the specified region.
*/
start MC_CDFMVN(N, b, Sigma, mu=j(1,ncol(Sigma),0));
   X = randnormal(N, mu, Sigma);
   inRegion = j(N, 1, 1);
   do i = 1 to ncol(b);
      v = (X[,i] < b[i]);       /* binary variable for the i_th coordinate */
      inRegion = inRegion & v;  /* find the intersection of all conditions */
   end;
   return mean(inRegion);       /* proportion of random points that satisfy criterion */
finish;

An Example: The 5-D Min Matrix

You can use the built-in COV2CORR function and the newly defined Xform_Limits_Cov2Corr function to transform a problem on the covariance (data) scale into an equivalent standardized problem. We will use a 5-dimensional "Min Matrix" as our covariance matrix. The Min Matrix is a structured matrix where the (i, j)-th element is equal to min(i, j). The mean vector, μ, and the integration limits, b, are arbitrary vectors.

The following program estimates the probability twice, both times using one million random variates in a Monte Carlo simulation. First, it uses the MC_CDFMVN function for the 5-D Min Matrix. Then it transforms the limits, creates the associated correlation matrix, and estimates the probability on the standardized scale.

So that the estimates are the same, and not just close to each other, you can use the RANDSEED subroutine to reset the random number stream prior to each simulation. This ensures that both computations use the exact same stream of pseudo-random numbers.

/* Example: use a 5-D "Min Matrix". 
   See: https://blogs.sas.com/content/iml/2026/02/16/min-matrix.html
*/
 
/* Define the parameters for MVN(mu, Sigma) */
Sigma = {1 1 1 1 1,
         1 2 2 2 2,
         1 2 3 3 3,
         1 2 3 4 4,
         1 2 3 4 5 };
mu = {1 2 3 4 5};
N = 1E6;
 
/* estimate the CDF probability at b for MVN(mu, Sigma)*/
b = {-1 2 2 3 7};
call randseed(1234, 1);    /* reset the RNG */
prob_cov = MC_CDFMVN(N, b, Sigma, mu);
 
/* standardize the parameters to the correlation scale */
c = Xform_Limits_Cov2Corr(b, Sigma, mu);
R = cov2corr(Sigma);
 
/* estimate the CDF probability at c for MVN(0, R)*/
call randseed(1234, 1);    /* reset the RNG */
prob_corr = MC_CDFMVN(N, c, R);
 
print prob_cov prob_corr;

The output shows that the probability estimate for the original problem is identical to the estimate for the standardized problem.

Summary

You can use techniques from linear algebra to extend the concept of a Z-score to multivariate normal distributions. You do not have to physically standardize each column of the data. Rather, you can convert the covariance matrix to the associated correlation matrix, apply a linear transformation to the limits of integration, and—voilà!—you can compute all multivariate normal probabilities in a standardized way.

This transformation is not purely academic. In computational statistics, algorithms that evaluate multivariate distributions are often implemented for the standardized correlation scale to ensure stable numerical computations. By scaling the problem and then computing answers on the correlation scale, the numerical routines tend to be more accurate and reliable.

The post On standardizing multivariate normal probabilities appeared first on The DO Loop.