Unfortunately, truncation has a known statistical problem: the density is biased near the boundary x=0 (Silverman, 1986). This article demonstrates an alternative: the log-KDE method. The log-KDE method guarantees that the estimated density is strictly positive and integrates to unity. Furthermore, the resulting density estimate is not biased near the boundary.

The Log-KDE method

Let X be a strictly positive random variable. We want to compute a nonparametric density estimate of the distribution of X. The log-KDE method consists of three steps:

1. Log-transform the strictly positive data: Y = log(X).
2. Fit a standard kernel density estimate to the transformed data, Y.
3. Back-transform the density to the original scale.

Step 3 is a little tricky. Let f_Y(y) be the KDE on the log scale. When you back-transform the density to the data scale, you must apply a change-of-variables formula, also known as the chain rule or the Jacobian of the transformation. This is familiar to students of integral calculus: when you change variables inside an integral, you also must change the differential.

For a general change-of-variables formula, consider the integral
∫ f_Y(y) dy
If you represent this integral in terms of x by making the substitution y = log(x), then dy = dy/dx dx = (1/x) dx. Thought of as measures, the back-transformed density is therefore
f_X(x) = f_Y(log(x)) (1/x)

Implement the method in SAS

Let's implement the log-KDE method for the same `Bulkhead` dataset used in the previous article. The data represent the position of a rivet relative to a bulkhead, which must be a positive quantity. The following DATA step defines the data.

data Bulkhead;
input Position @@;
datalines;
 2.44  6.84 16.84  3.76  5.58 14.53 12.12  8.00  8.66  6.63  9.51  6.37 11.78
 6.14  9.04  2.88 21.28  8.43 18.91  4.45  4.59 20.42  9.03  6.55  7.15  4.10
 8.72  5.17  2.41 10.93 14.08  5.44  2.81  3.99  8.40 11.28  3.97 18.99 10.51
 9.52  2.01  1.5 
;
%let DSName  = Bulkhead;
%let varName = Position;

First, use the DATA step to log-transform the data, Then, use PROC KDE to estimate the density on the unbounded log scale. You can use the OUT= option to save the evaluation points and density values to a dataset, as follows:

/* 1. Log-transform the data */
data LogData;
  set &DSName;
  Log_&varName = log(&varName);
run;
 
/* 2. Fit a density estimate to the transformed data */
proc kde data=LogData;
   univar Log_&varName / unistats 
          out=KDEOUT_Log(rename=(value=log_value density=KDE_LogDensity));
run;

The 'KDEOUT_Log' data set contains the density (KDE_LogDensity) of the log-transformed variable (log_value). Next, use the DATA step to apply the change-of-variables formula, which maps the log-density back to the data scale:

/* 3. Back-transform the density using f_X(x) = f_Y(log(x)) * (1/x), where 1/x is the Jacobian */
data KDEOUT;
   set KDEOUT_Log;
   x = exp(log_value);
   KDE_Density = KDE_LogDensity / x;  
   drop var log_value KDE_LogDensity;
run;

Visualize the back-transformed KDE

You can overlay the newly calculated density curve (in the KDE_Density variable), onto a histogram of the original data. As shown in a previous article, you can overlay custom density curves on a histogram by using the %EmulateHistogram macro.

/* Create macro variables for histogram bins and observation count.
   Download the macro from https://blogs.sas.com/content/iml/2025/10/13/high-low-emulate-histogram.html */
%EmulateHistogram(dsIn=&DSName, varIn=&varName);
 
/* Scale the density to the percentage (or count) scale of the histogram */
data All;
   set &DSName _HistBins KDEOut;
   KDE_Percent = 100    * &_binWidth * KDE_Density;
   KDE_Count   = &_NOBS * &_binWidth * KDE_Density;
run;
 
/* 4) Overlay on a histogram */
title "Back-Transform of Fit of log(Position)";
proc sgplot data=All noautolegend;
   highlow x=_midpt_ low=_zero_ high=_pct_ / type=bar barwidth=1;
   series x=x y=KDE_Percent;
   fringe &varName;
   yaxis min=0 offsetmin=0 grid;
   xaxis values=(0 to &_binEnd by &_binWidth);
run;

Notice that the KDE is zero when the Position variable is negative. Furthermore, this KDE integrates to unity. The estimate of the KDE near zero is data dependent. For these data, the KDE drops smoothly to zero at Positon=0. Contrast that with the estimate in the previous article, which predicted a nonzero probability when Position ≈ 0.

Peculiarities of the log-KDE method

It is important to understand that this log-KDE estimate is different from the standard KDE on the data scale. Most KDEs use a constant bandwidth for the kernels. Intuitively, you can visualize the KDE as the sum of a bunch of little bump functions (the kernels) centered at each data point. The most common bump functions a Gaussian in shape. But that is not true for the log-KDE method.

For the log-KDE method, Gaussian kernels are used to fit the density on the log scale. But when you back-transform the density, the Gaussian kernels are transformed into lognormal kernels. (Recall that the definition of a lognormal variable is one whose logarithm is normal.) Furthermore, the bandwidth of the lognormal kernels is not uniform.

To see this, consider a Gaussian kernel centered on the point y_i = log(x_i) in the log-scale coordinate system. It has a bandwidth, h, which determines the scale of the kernel. The back-transformation maps the Gaussian kernel to a lognormal kernel that has location parameter log(x_i) and scale parameter h. The standard deviation of that lognormal distribution is therefore
exp(log(x_i) + h²/2) sqrt(exp(h²) - 1)
or
x_i exp(h²/2) sqrt(exp(h²) - 1)

The interpretation of this formula is that the width of the (lognormal) kernels on the data scale are proportional to the data value, x_i. For data values much less than 1, the kernels on the data scale are tall and narrow whereas for values much greater than 1, the kernels on the data scale are short and broad.

For example, suppose you have the data values {0.2, 0.5, 1, 2}. These get log-transformed and the KDE is computed on the log scale by placing a Gaussian density at log(x_i) for each i. Now consider what these bell-shaped curves will look like when they are back-transformed. Each becomes a lognormal curve, but the width of those curves is not constant. This is illustrated in the following graph, which shows the kernels on the data scale for the data values {0.2, 0.5, 1, 2}.

As you can see, the kernel based at x=0.2 is tall and narrow, whereas the kernel based at x=2 is broad and flat.

Because of this, you might notice spikes in the log-KDE for data that are near 0. For example, if you were to add a new data point to the Bulkhead data set for Position=0.5, the log-KDE will become bimodal and exhibit a new wiggle near that location.

Summary

This article shows how to use the log-KDE method to estimate the density of strictly positive data. The method is guaranteed to estimate a density who support is strictly positive. Furthermore, the resulting density estimate is not biased near the boundary. However, geometrically, this is achieved by using lognormal kernels with nonconstant bandwidths. The standard deviation of a lognormal kernel is proportional to the data value at which the kernel is based.

Further Reading

Jones, A. T., Nguyen, H. D., & McLachlan, G. J. (2018). "Kernel density estimation on positive data via the logKDE package for R." The logKDE package is available on CRAN.

The post The log-KDE method: Density estimates for positive data appeared first on The DO Loop.

However, the KDE's flexibility comes at a cost. If the underlying distribution is known to have certain properties, the KDE does not necessarily preserve those properties. A property that arises often in practice is positivity. Many real-world measurements are strictly positive, such as mass, length, and many health-related clinical lab values. If you have a distribution of lengths, and you fit a KDE to the data, it is disconcerting (and wrong!) to display a density estimate that predicts positive density for negative lengths!

There are two ways to deal with this problem. This article shows how to use SAS to truncate the density estimate by setting the density to zero when the variable is negative. A subsequent article shows how to log-transform the data, fit a KDE on the log scale, and back-transform the results in SAS.

Fit a KDE to positive data

The following data set shows the position (in millimeters) of a rivet relative to a bulkhead. The distance is a positive quantity. For convenience, I assign the name of the data sets and the variables to macro variables. The remainder of the examples use only the macro variables, which makes it easy for you to run the code on your own data.

data Bulkhead;
input Position @@;
datalines;
 2.44  6.84 16.84  3.76  5.58 14.53 12.12  8.00  8.66  6.63  9.51  6.37 11.78
 6.14  9.04  2.88 21.28  8.43 18.91  4.45  4.59 20.42  9.03  6.55  7.15  4.10
 8.72  5.17  2.41 10.93 14.08  5.44  2.81  3.99  8.40 11.28  3.97 18.99 10.51
 9.52  2.01  1.5
;
 
%let DSName  = Bulkhead;
%let varName = Position;

You can use either PROC UNIVARIATE or PROC KDE to fit a kernel density estimate. I will use PROC KDE because, as we will see in the next section, it supports options to truncate the KDE. The following call fits the KDE to the data and visualizes the result. The UNISTATS option displays descriptive statistics for the data and shows the bandwidth chosen for the fit:

proc kde data=&DSName;
   univar &varName / unistats;
run;

I have highlighted a few entries in the output tables. Although the range of the data is [1.5, 21.28], the Controls table shows that the range of the KDE is [-4.99, 27.77]. The graph of the KDE curve confirms this fact. The KDE extends about three bandwidths past the extremes of the data, and the kernel bandwidth for these data is 2.13. However, it does not make sense to estimate the density to be nonzero for Position < 0 because Position must be a positive quantity.

PROC KDE provides options that enable you to truncate the density estimate in two ways:

• Truncate the KDE to the range of the data by using the TRUNCATE option.
• Truncate the KDE to an arbitrary range by using the GRIDL= option to specify a lower limit and/or the GRIDU= option to specify an upper limit.

It is somewhat unfortunate that PROC KDE does not have an option to display the histogram and KDE on the percentage scale. Most graphs in this article are shown on the frequency scale. You can rescale the graphs to display them on the percent or density scale. The last graph in this article is shown on the percentage scale.

Truncate the density into the data range

The following call uses the TRUNCATE option to restrict the KDE to the range of the data:

/* truncate KDE to data range */
proc kde data=&DSName;
   univar &varName / truncate out=KDE_truncData; /* save KDE to data set for later analysis */
run;

Notice that the KDE is truncated on the interval [1.5, 21.28], which is the range of the data.

Truncate the density onto the positive half-line

The Position variable cannot be negative, but there is no intrinsic reason why it cannot be smaller than 1.5. There is also no obvious reason to bound the density from above. An alternative way to truncate the density is to use the GRIDL=0 option to truncate the density to strictly positive values.

proc kde data=&DSName;
   univar &varName / gridl=0 out=KDE_trunc0; /* save KDE to data set for later analysis */
run;

This KDE is restricted to positive values. It is unbounded above, but the density is exactly zero for negative values of the Position variable. If you have an upper limit, you can set it by using the GRIDU= option. For example, if you are plotting the density of student scores on an exam that has 100 as a perfect score, you could use GRIDL=0 GRIDU=100 to bound the KDE to the range [0, 100].

Limitations of the truncation method

There are two problems with truncating a KDE:

• On a probability scale, the area under a true density must integrate to 1. If you tell PROC KDE to truncate the KDE, it does so by throwing away some area. Thus, area is less than 1 for the restricted KDE.
• Even if the estimate is rescaled, the KDE underestimates the expected value of distribution near 0 (Silverman, 1986, Section 2.10).

You can rescale to the density to correct for the first problem. In the previous sections, I wrote the values of the KDE to the data sets KDE_truncData and KDE_trunc0, respectively. The following macro computes the area under these curves by using a trapezoidal approximation to the area.

/* the area under a density estimate should be 1. See
   https://blogs.sas.com/content/iml/2011/07/08/the-area-under-a-density-estimate-curve-nonparametric-estimates.html
   The following macro computes the area under a KDE estimate that is created by PROC KDE and written
   to the 'DSName' data set by using the OUT= option on the UNIVAR statement.
   In the data set the x variable is named 'value', and the 'density' variable is f(x).
   */
%macro AreaUnderKDE(DS);
data _null_;
   set &DS end=EOF;
   dx = dif(value);
   meanY = (density + lag(density)) / 2;
   Area + dx*meanY;
   if EOF then do;
      call symputx('Area', Area, "G"); /* create global macro */
   end;
run;
%mend;
 
%AreaUnderKDE(KDE_truncData);
%put Trunc Data &=Area;
%AreaUnderKDE(KDE_trunc0);
%put Trunc to 0 &=Area;

Trunc Data AREA=0.9031115838
Trunc to 0 AREA=0.97522813

The log shows that the area is less than 1 for both KDE curves. This can be a problem if you intend to use the KDE to estimate probabilities. However, it is simple enough to rescale the truncated density by dividing it by its area:

/* If 'density' is a variable that does not integrates to 1,
   you can renormalize by dividing by the Area under the density curve. */
%macro Renormalize(DS, Area);
   data &DS;
      set &DS;
      density = density / &Area;
   run;
%mend;
 
%Renormalize(KDE_trunc0, &Area);

The following graph overlays the renormalized truncated KDE on a histogram and fringe plot of the data. The KDE is truncated at 0 because the measured quantity must be positive. The graph was created by using the %EmulateHistogram macro and the process described in the article, "Overlay multiple custom density curves on a histogram in SAS."

The second issue (underestimating the KDE near the boundary) does not have a simple solution. Some researchers suggest reflecting the data across x=0 before fitting the KDE. I don't want to say more about that method here, but I have used it previously to create smoothers of periodic data.

Summary

If you naively compute a kernel density estimate for data that are inherently positive, the KDE might estimate density for negative values of the variable. PROC KDE in SAS provides two methods for correcting this problem:

• The TRUNCATE option on the UNIVAR statement truncates the KDE to the range of the data.
• The GRIDL=L and GRIDU=U options truncate the KDE to the range [L, U]. For example, you can use GRIDL=0 to truncate the KDE onto the positive axis.

Technically, if you truncate the KDE, you should renormalize it so that the area under the curve is unity. This article includes two macros (%AreaUnderKDE and %Renormalize) that you can use for that purpose. I do not deal with the fact that the truncated KDE underestimates the expected density near the boundary. In a subsequent article, I present an alternative method: log-transforming the data.

The post Kernel density estimates for positive data appeared first on The DO Loop.

This article describes an alternative to the standard Pareto chart and shows how to create the alternative version in SAS.

A problem with the standard Pareto chart

Not every situation obeys the 80-20 rule. For your process, you might find that 80% of the defects are produced by 33% of the causes. Or, maybe 75% of the defects are produced by 10% of the causes. If your process is modeled by a Pareto distribution (a specific power-law distribution), you can use SAS to fit your data to a Pareto distribution. However, most practitioners do not do that. They just display a bar chart of the causes where the bars are displayed in descending order and then try to address the most frequent causes. For example, if a quality engineer looks at the Pareto chart to the right, she might decide to address the "contamination" cause first, followed by the "oxide defect" issue. Together, these two categories account for 67.5% of the defects.

But there is a fundamental problem with this strategy. In small samples, the empirical frequencies in the data might be much different than the underlying probabilities for the process. Thus, the "most frequent" category in the data might not be the most probable cause of the defects. It would be nice to know whether the tall bars to the left just happen to be tall because of random variation or whether they are actually the most probable categories and warrant your attention.

In a short three-page paper, Wilkinson (2006, TAS, "Revising the Pareto Chart") proposes an alternative to the Pareto chart that uses the idea of statistical significance to give additional insight into the tallest bars. The Wilkinson Pareto chart provides acceptance intervals that indicate whether each category appears more often than would be expected if all categories are equally probable. This enables you to assess whether the tallest bars are "sufficiently tall" and are, in fact, more probable than the other bars.

A Pareto chart with acceptance intervals

If you have a license for SAS/QC software, you can create Wilkinson's alternative Pareto chart by using the CHARTTYPE=INTERVALS option on the VBAR statement in PROC PARETO. A previous article showed how to create two standard Pareto charts in SAS. Let's analyze the same data set, which contains the causes of 40 randomly selected defective components in a manufacturing process. The following DATA step defines the data. The call to PROC PARETO creates Wilkinson's acceptance-interval chart:

/* Modification of PROC PARETO example data */
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
;
 
%let DSName = Failure;
%let VarName = Cause;
 
/* Wilkinson variant of Pareto chart. Which categories are 
   significantly different from a uniform distribution of 
   causes? */
proc pareto data=&DSName;
   vbar &VarName / charttype=intervals;
run;

What exactly are we looking at? The documentation for PROC PARETO states, "the most frequently occurring problem, Contamination, occurs more frequently than the first-ranked cause from a random sample of ... uniformly distributed causes. This result indicates that addressing contamination problems should be given a high priority."

The documentation does not explain the term "first-ranked cause" or the connection to "uniformly distributed causes."

Understanding the Pareto chart with acceptance intervals

The chart shows both the observed percentages and the expected percentages under the null hypothesis that all categories are equally likely to occur. However, the fact that a Pareto chart is sorted, introduces a complication into what would otherwise be a simple process.

This example has 6 categories. If the categories are equally likely, we would expect each bar in a bar chart to contain about 16.67% of the defects. In the observed sample, the most frequent category (Contamination) contains 42.5% of the defects. Intuitively, we suspect that 42.5% is so much greater than 16.67% that this number is unlikely to be observed by chance if the categories are equally likely. On the other hand, the second most frequent category (Oxide Defect) contains 25% of the defects, which is not too much larger than 16.67%. Perhaps that category is no more likely than any of the other remaining categories?

Wilkinson's chart visualizes the comparison between the observed percentages and the null hypothesis where every category is equally likely. For each category, the graph shows two features:

• The dots show the observed percentage of defects in each category. Thus, the first category (Contamination) shows a dot at 42.5%, the second shows a dot at 25%, the third shows a dot at 12.5%, and so forth.
• The bars represent 95% confidence intervals under the null hypothesis that all six causes are equally probable and under the assumption that the bars are being displayed in decreasing order of size. Wilkinson calls these the acceptance regions. If the observed percentage is outside and above the acceptance region (as it is for the Contamination category), it indicates that the most-frequent category is observed significantly more often than would be expected for the most-frequent category under the null scenario. If the observed percentage is inside the acceptance region (for example, Oxide Defect, Metallization, and Silicon Defect), it indicates that the observed frequency is not different from what would be expected under the null scenario. If the observed percentage is below the acceptance region (Corrosion and Doping), the observed frequency is less than what would be expected under the null scenario.

When I first saw this plot, I was a little confused because I expected the acceptance regions to be centered around 16.7%, which is the expected value for six equally probable categories. However, they are not! The reason for this discrepancy is that the acceptance regions represent the range of 95% of the ranked categories under the hypothesis that all categories are equally probable.

If all categories were equally probable, random variation causes some categories to have more observations than others. When you sort the categories, the first-ranked bar is almost always taller than 16.7%. Similarly, the last-ranked bar is almost always less than 16.7%. Thus, you must adjust the heights of the acceptance regions to account for the fact that the categories are shown in decreasing order.

Using simulation to create the acceptance intervals

You can use simulation to create the acceptance intervals, which are the ranges of the first-ranked, second-ranked, and third-ranked (etc.) categories when the categories are selected uniformly at random from a set of size N. In this example, N=40 and there are K=6 categories. You can use the RANDMULTINOMIAL function in SAS IML software to simulate a large number (B=5000) of samples. For each sample, sort the counts in descending order to mimic the method used by the Pareto chart. You can do these two steps efficiently by simulating all B samples into a matrix that has B rows and K columns. You then sort each row independently.

After this procedure, the first column of the sorted matrix contains the distribution of counts for the first-ranked category in each sample. In general, the i_th column contains the distribution of counts for the i_th-ranked category. You can therefore compute the 2.5th and 97.5 percentiles to form a 95% confidence interval for the counts. This procedure is implemented in the following IML program:

proc iml;
call randseed(1234);
N = 40;                /* Total number of observations */
K = 6;                 /* Number of categories */
alpha = 0.05;          /* Significance level */
B = 5000;              /* Number of Monte Carlo iterations */
 
/* Monte Carlo simulation from the Null distribution, which assumes uniform probability */
prob = j(K, 1, 1/K);
Y = RandMultinomial(B, N, prob);
/* Sort each row. Columns contain frequencies for each rank (1st, 2nd, 3rd,...).
   You could also use the SortMat function; see 
   https://blogs.sas.com/content/iml/2026/06/24/sort-rows-cols-matrix.html
*/
do i = 1 to nrow(Y);
   v = colvec(Y[i,]);
   call sort(v, 1, 1); /* sort descending */
   Y[i, ] = rowvec(v);
end;
w = (alpha/2) // (1-alpha/2);
call qntl(CL, Y, w);              /* 95% CL of counts for ranked categories */
/* transpose CL to make it easier to visualize */
CL = CL`;
print CL[c={'LowerCL' 'UpperCL'} r=('R1':'R6')];
QUIT;

The output shows the confidence intervals for counts when a set of K=6 categories are chosen uniformly at random (with replacement) N=40 times. The first-ranked category (the most-frequent one) usually has 8-13 counts. This is larger than the 40/6 = 6.67 counts that you might have expected. Similarly, the second-ranked category usually has between 7-10 counts. The last two categories usually have 6 or fewer counts. Obviously, you can divide these counts by N=40 to get the corresponding percentages. You can then write the percentages to a SAS data set and use a HIGHLOW statement in PROC SGPLOT to reproduce Wilkinson's Pareto chart.

A SAS macro to create acceptance regions for a Pareto chart

If you have a license for SAS IML software, the following SAS macro computes the same alternative Pareto chart as the CHARTTYPE=INTERVALS option in PROC PARETO.

/* Plot the observed counts as a dot plot and overlay an acceptance band
   composed of the lower and upper 95% limits for the null distribution.
   This macro requires license for SAS IML.
   EXAMPLES:
   %AcceptancePareto(sashelp.cars, Type);   * 95% acceptance regions;
   %AcceptancePareto(sashelp.cars, Cylinders, alpha=0.1);   * 90% acceptance regions;
*/
%macro AcceptancePareto(DSName, VarName, alpha=0.05);
   %local conf;
   %let conf = %sysevalf(100*(1-&alpha));
 
   proc freq data=&DSName order=Freq noprint;
      where not missing(&VarName);
      table &VarName / out=_FreqOut outcum;
   run;
 
   /* compute upper and lower CL */
   proc iml;
   call randseed(1234);
   use _FreqOut;   read all var "count";   close;
 
   N = sum(count);                /* Total number of observations */
   K = nrow(count);               /* Number of categories */
   alpha = α                /* Significance level */
   B = 5000;                      /* Number of Monte Carlo iterations */
 
   /* Monte Carlo simulation from the Null distribution, which assumes uniform probability */
   prob = j(K, 1, 1/K);  /* H0: All categories are equally probable */
   Y = RandMultinomial(B, N, prob);
   /* Sort each row. Columns contain frequencies for each rank (1st, 2nd, 3rd,...) */
   do i = 1 to nrow(Y);
      v = colvec(Y[i,]);
      call sort(v, 1, 1); /* sort descending */
      Y[i, ] = rowvec(v);
   end;
   w = (alpha/2) // (1-alpha/2);
   call qntl(CL, Y, w);              /* 95% CL of counts for ranked categories */
   /* transpose CL and convert to percentage */
   CL = 100 * CL` / N;
   create _CL from CL[c={'_Lower' '_Upper'}];
      append from CL;
   close;
   QUIT;
 
   data _AcceptancePareto;
   merge _FreqOut _CL; 
   label _Lower = "Lower &conf.% Acceptance Percent"
         _Upper = "Upper &conf.% Acceptance Percent"
         Percent = "Observed Percent";
   run;
 
   proc sgplot data=_AcceptancePareto;
      highlow x=&VarName low=_Lower high=_Upper / 
              type=bar barwidth=0.8 legendlabel="&conf.% Acceptance Interval";
      scatter x=&VarName y=Percent;
      yaxis  min=0 offsetmin=0 grid label="Percent";
      xaxis type=discrete discreteorder=data;
      keylegend / location=Inside position=NE opaque across=1;
      format Percent best4.;
   run;
%mend;
 
title "Revised Pareto Chart with 95% Acceptance Bands";
title2 "H0: Uniform Probability for Causes";
footnote J=L "Wilkinson (2006, TAS)";
%AcceptancePareto(Failure, Cause);

I've added grid lines and a legend, but the information in this plot is the same as the one produced by PROC PARETO. The graph shows that the empirical probability of the Contamination defects is significantly different from what would be expected if the defect causes were equally probable. Thus, the quality engineer should strive to reduce the contamination defects. The observed count for Oxide Defect is not significantly higher than expected. It is on the boundary of the acceptance region, so the engineer might want to gather more data before dedicating resources to address the Oxide Defect.

Summary

This article shows how to create an alternative to the standard Pareto charts (Wilkinson, 2006). The alternative chart provides acceptance intervals that indicate whether each category appears more often than would be expected if all categories are equally probable. This enables you to assess whether the tallest bars are "sufficiently tall" and are, in fact, occurring more often than chance. If not, there is little reason to spend time and money trying to address that cause.

The post An alternative Pareto chart in SAS appeared first on The DO Loop.

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.