The ideas behind the orthogonalization process are simple: Project a vector onto a linear subspace and compute the residual vector. The residual vector is orthogonal to the subspace. Thus, these two operations decompose any vector into two parts: A vector that lies in the subspace and a vector that is orthogonal to the subspace. This article shows two subroutines in the SAS IML language that are useful for projecting a vector into a subspace and finding an orthogonal vector to a subspace. The GSORTH subroutine implements the Gram-Schmidt orthogonalization of a matrix, whereas the ORTVEC subroutine performs one step of the process and is useful for finding a orthogonal vector to a subspace.

The Gram-Schmidt orthogonalization of a matrix

Suppose V1, V2, ..., Vk are k linearly independent vectors. Let S = span(V1, V2, ..., Vk) be the linear subspace that they span. The Gram-Schmidt process uses these vectors to construct a set of orthogonal vectors for S. In SAS IML, you can put the vectors into the columns of a matrix and call the GSORTH subroutine. The subroutine takes one input argument (a matrix, V) and returns three output matrices. The output are as follows:

1. Q: An orthonormal matrix whose columns contain an orthonormal basis for the subspace, S .
2. T: An upper triangular matrix that transforms the orthonormal basis into the original basis: V = Q*T.
3. lindep: A binary scalar value. If lindep=1, the columns of V are linearly dependent. If lindep=0, the columns of V are linearly independent, which implies dim(S) = k.

Here's an example. The columns of V span a 3-D linear subspace in the R⁴ Euclidean space.

proc iml;
/* the columns of V span a 3-D subspace in R^4 */
V = { 1 1  0,
     -1 0  2,
     -1 1 -1,
      1 1  2 };
 
call gsorth(Q, T, lindep, V);  /* output: Q, T, and lindep from the input matrix, V */
print Q[F=BestD6.], T[F=BestD6.];
/* verify that V = Q*T by computing the difference V - Q*T */
maxDiff = max(abs(V - Q*T));
print maxDiff;

The original three basis vectors are V1={1, -1, -1, 1}, V2={1, 0, 1, 1}, and V3={0, 2, -1, 2}. They span a subspace, S, in R⁴. After calling the GSORTH routine, the columns of the matrix Q contain an orthogonal set of basis vectors that also span S. (To standardize the output, the columns of Q have unit norm, which means they form an orthonormal basis.) The T matrix maps one set of basis vectors onto the other. Specifically, V = Q*T is the equation that represents each column of V as a linear combination of the columns of Q. For this example, the columns of V are linearly independent, so lindep=0 and dim(S) = 3.

Notice that this Q matrix is not square, so it cannot be called an orthogonal matrix. Instead, a rectangular matrix with orthonormal columns is called a semi-orthogonal matrix. If, in addition, the columns have unit length, it is called a semi-orthonormal matrix.

Computing orthogonal vectors to a linear subspace: The manual way

The GSORTH subroutine enables you to compute an orthonormal basis with a single call. In linear algebra texts, the Q matrix is usually constructed sequentially. At each step, you create a new orthonormal basis vector from one of the original spanning vectors. To understand how the GSORTH and ORTVEC subroutines work, it helps to look at the manual calculations. In my linear algebra class, the Gram-Schmidt orthogonalization process was used so often that it is indelibly etched into my brain from performing endless calculations by hand! This section shows the tedious manual computations. The subsequent section uses the ORTVEC subroutine to simplify these computations.

The Gram-Schmidt process builds an orthogonal basis incrementally. Let's look at the SAS IML code. The original vectors are V₁, V₂, and V₃, which are columns of a matrix, V. We use them to construct a new set of orthogonal vectors (U₁, U₂, and U₃), which are columns of a semi-orthogonal matrix.

The process works step-by-step, projecting vectors onto the existing subspace and calculating the residual:

• Step 1: The first orthogonal vector, U₁, is simply the first original vector.
  U₁ = V₁
  If you want to, you can standardize U₁ construct an orthonormal basis, but let's keep it simple for now. We'll standardize the U vectors at the end.
• Step 2: The second orthogonal vector, U₂, is the second original vector minus its projection onto the first orthogonal vector. Geometrically, you decompose V₂ into two components, one lying in the span of U₁ and the other orthogonal to it. The orthogonal component becomes U₂:
  U₂ = V₂ - proj(V₂, U₁)

• Step 3: The third orthogonal vector is the third original vector minus its projections onto the first two orthogonal vectors. Geometrically, you decompose V₃ into two components, one lying in span(U₁, U₂) and the other orthogonal to it. The orthogonal component becomes U₃:
  U₃ = V₃ - proj(V₃, U₁) - proj(V₃, U₂)

In these equations, the projection operator proj(y, x) calculates the projection of a vector y onto the line spanned by vector x. Let's implement these three steps of the G-S orthogonalization process in the SAS IML language. To make sure we get the same answer as the GSORTH subroutine, we can standardize the columns of U to form Q:

/* manual operation: project y onto span(x) */
start proj_vec(y, x);
   z = x/ norm(x);        /* unit vector in direction of x */
   return( (y`*z) * z );
finish;
 
U = j(nrow(V), ncol(V), .);
U[,1] = V[,1];
U[,2] = V[,2] - proj_vec(V[,2], U[,1]);
U[,3] = V[,3] - proj_vec(V[,3], U[,1]) - proj_vec(V[,3], U[,2]);
print U;
 
/* optionally normalize the columns of U at each step or at the end. */
norm_vec = j(1, ncol(U), .);
Q = j(nrow(U), ncol(U), .);
do j = 1 to ncol(U);
   norm_vec[j] = norm(U[,j]);
   Q[,j] = U[,j] / norm_vec[j];
end;
print norm_vec[F=Best5.], Q[F=BestD6.];

Notice that the Q matrix from these steps is identical to the Q matrix that is produced by the GSORTH subroutine. The U matrix is an "un-standardized" version of Q.

Computing orthogonal vectors to a linear subspace: The easier way

As described in the previous section, the Gram-Schmidt process is a series of operations that projects vectors onto a linear subspace, computes a vector orthogonal to the subspace, and then iterates. The ORTVEC subroutine encapsulates one step in this iterative process. It standardizes the orthogonal vectors at each step, thus forming the Q matrix directly instead of first forming U and then standardizing the columns.

The syntax for the ORTVEC routine is a little complicated. Let n be the dimension of the vectors (in our example, n=4). The input arguments for the ORTVEC subroutine are an n x 1 vector, v, and a semi-orthonormal n x k matrix, Q, where k ≤ n. This syntax is used for finding the (k+1)th orthonormal basis vector, based on the previous k orthonormal vectors. The vector v is the (k+1)th vector in the original basis. The ORTVEC routine returns four values: two vectors, w and r, a scalar rho, and a binary flag that tells you whether v is in the span of the columns of Q. The return values provide the decomposition of v into a component inside span(Q) and a component orthogonal to span(Q). In symbols, v = Q*r + rho*w, where norm(w)=1. The vector w is orthogonal to the span of the columns of Q. The vector w is the most important output from the ORTVEC subroutine.

To get the first orthonormal basis vector, you omit the Q matrix. This is shown in the following example, which uses the ORTVEC subroutine to compute the same Q vector as in the previous section:

/* the ORTVEC function in SAS IML performs several linear algebra operations
   related to G-S orthogonalization */
Q= j(nrow(V), ncol(V), .);
call ortvec(w1,r,rho1,lindep, V[,1]);  /* Step 1: No Q matrix yet; rho = norm(v1); v1 = rho1*w1 */
Q[,1] = w1;
/* Check: v1 = rho1*w1; print (v1-V[,1]); */
 
call ortvec(w2,r2,rho2,lindep, V[,2], Q[,1]); /* Step 2: v2 = Q1*r2 + rho2*w2 */
Q[,2] = w2;
/* Check: v2 = Q[,1]*r2 + rho2*w2; print (v2-V[,2]); */
 
call ortvec(w3,r3,rho3,lindep, V[,3], Q[,1:2]); /* Step 3: v3 = [Q1 Q2]*r3 + rho3*w3 */
Q[,3] = w3;
/* Check: v3 = Q[,1:2]*r3 + rho3*w3; print (v3 - V[,3]); */
print Q;

You now have three equivalent ways to construct the Q matrix. If you are interested in a programming challenge, you could put these ORTVEC operations into a DO loop (or a function) that performs the entire Gram-Schmidt orthogonalization algorithm for an arbitrary matrix, V. For an extra challenge, print out the sequence of 'r' vectors and the 'rho' scalars. Can you determine how they are related to the matrix T that is provided by the GSORTH routine?

Summary

The Gram-Schmidt orthogonalization algorithm enables you to construct an orthogonal basis for the span of other linearly independent vectors. The orthogonalization process decomposes every vector into a component that lies inside a subspace and a vector that is orthogonal to the subspace. In the SAS IML language, you can use the GSORTH subroutine to perform the full algorithm, or you can use the ORTVEC subroutine to perform one step. There are many applications that require finding a vector that is orthogonal to a subspace, so the ORTVEC subroutine is useful in these applications.

The post Constructing orthogonal vectors in SAS appeared first on The DO Loop.

Recently, I read a paper by Patrick Hoyer (2004, "Non-negative matrix factorization with sparseness constraints"), that imposes constraints on the NMF method. The constraints are designed to make the factors more sparse, which is why Hoyer introduced his sparseness measure for vectors. Briefly, Hoyer's algorithm results in nonnegative factors that might be more interpretable.

This article shows an example by revisiting the Scotch whisky data, which I analyzed by using NMF in a previous article. On paper, Hoyer's idea seemed both interesting and powerful. However, after running a sparseness-constrained NMF on an example, I was not impressed. Read on to find out why. To try out the Hoyer algorithm yourself, download the SAS program that implements a sparseness-constrained NMF factorization.

An unconstrained feature set for Scotch whiskies

The previous article used the NMF to extract a four-dimensional basis of vectors for Scotch whisky data. It is known that whiskies often have regional characteristics, such as the "earthy" or "peaty" whiskies from the island of Islay. Ideally, the NMF features would reflect these regional flavor characteristics.

The whisky data can be represented as an 86x12 nonnegative data matrix, X. Each row represents flavor characteristics of a whisky. Each whisky is assigned a flavor profile that consists of the values 0-4 for the 12 flavors. The 12 flavors are Tobacco, Medicinal, Smoky, Body, Spicy, Winey, Nutty, Honey, Malty, Fruity, Sweetness, and Floral.

The NMF creates a rank-k approximation to X. This article uses k=4. Thus, X is approximated as a product W*H, where W is 86x4 and H is 4x12. The rows of H are the "flavor profile" feature vectors. They form a basis for a 4-D linear subspace. They are analogous to the principal components (eigenvectors) in a PCA analysis. Each row of the W matrix is a set of weights (or coefficients) for the 4-factor approximation of the corresponding row of X.

The following image shows the transposed H matrix. The columns in the image visualize a 4-D basis for the unconstrained NMF.

The columns of H` are the features found by the NMF. The NMF method is not unique, but in this image:

1. The first NMF factor is a combination of six flavors: Smoky, Spicy, Malty, Fruity, Sweetness, and Floral variables.
2. The second NMF factor is a combination of 10 flavors.
3. The third factor is primarily composed of the Tobacco, Medicinal, Smoky, and Body flavors.
4. The fourth factor is primarily composed of Nutty flavor, with smaller amounts of Malty and Sweetness.

In terms of Hoyer's sparseness measure, the first and second features are not sparse. Hoyer's measure ranges from 0 (not at all sparse) to 1 (very sparse). The following table shows that the first NMF feature has moderate sparseness (0.459). The second feature has very little sparseness (0.148). The third and fourth features are strongly sparse (0.601) and very strongly sparse (0.822), respectively.

Before leaving the visualization of the unconstrained NMF, I want to draw your attention to the first row of H`. The third factor contains a small amount of the Tobacco variable. The cell is light orange. If you print out the H matrix, you can see that the first and second factor also contain a trace amount of the tobacco variable. However, compared to the other variables, this rank-4 NMF approximation does not represent the tobacco variable well. It seems that the tobacco variable probably points in a direction that is primarily orthogonal to the span of these four NMF factors.

A sparseness-constrained NMF algorithm

There are 12 variables and k=4 factors, so it would be nice if features were combinations of no more than four or five variables. If a 12-D vector has 4 nonzero components and 8 zero components, Hoyer's sparseness measure is about 0.6. Let's use that value as a target value. What will the value 0.6 do? The unconstrained factors naturally ranged from 0.15 to 0.82. The sparseness constraint forces the dense factor (0.15) to become much sparser. However, it will actually reduce the sparseness of the fourth factor!

I implemented Hoyer's sparseness-constrained NMF in SAS IML. The following image shows the result of running this constrained NMF algorithm on the Whisky data and constraining the rows of the H matrix to have sparseness equal to 0.6:

In the new NMF factorization, the sparseness of EVERY feature (rows of the H matrix) is exactly 0.6. In the new factorization, all features depend on 4-5 variables. Notice that the second feature, which used to be a blend of 10 flavors, is very different. Some of the flavors that were previously captured in the second feature are now part of other factors. For example, the Winey contributions have moved to the fourth factor, and the Honey contribution has moved to the first factor.

If you look closely at the first row in the image, you'll see that no factor contains the Tobacco variable! In this rank-4 sparseness-constrained factorization, the sparseness constraint forces the columns of H` to be formed by using only a few variables. The Tobacco variable is not included in any sparse factor. Thus, this model predicts that no Scotch whisky has a Tobacco flavor.

If you want to include the Tobacco variable in the NMF factors, you need to relax the sparseness constraint (or use a higher-rank approximation). For example, setting the sparseness target to 0.5 will result in factors that include Tobacco, but the resulting factors each contain six or more variables and are not as interpretable.

This is a specific instance of a general truth: The unconstrained NMF factorization will always approximate the data matrix better than a sparseness-constrained factorization. Sometimes, the lack of fit might be qualitatively noticeable, as in this example where the Tobacco variable is not included in any NMF factor.

The tradeoffs between sparseness and goodness of fit

The sparseness-constrained NMF is always a poorer approximation to the data matrix. Thus, you have to decide whether the trade-offs are worth it:

• Is it better to use the unconstrained NMF? The fit is better, but some factors might be formed as a linear combination of many variables, which hampers interpretability.
• Is it better to use the sparseness-constrained NMF? The fit is worse, but you can force the factors to be sparse and therefore highly interpretable. On the other hand, the sparseness constrained NMF is a Procrustean solution to the interpretability of the factors. By forcing every factor to have the same sparseness, it ignores the natural data-induced correlations between variables.

My general preference is to let the data decide how many variables are contained in each NMF factor. The unconstrained NMF approximation seems more faithful to the data. The factors in the unconstrained NMF are mostly interpretable. Their sparseness ranges from 0.148 (not sparse) to 0.822 (strongly sparse). It isn't terrible to have a factor whose interpretation is "an overall blend of flavors." After all, every observation is a Scotch whisky, so it isn't surprising to have a common factor that describes flavors that are common to many whiskies.

In the sparseness-constrained algorithm, every factor has exactly the same sparseness. Every factor depends strongly on four or five flavors. This does not seem realistic. If regional distilleries produce whiskies that have one or two special flavor characteristics, it seems better to have factors that are linear combinations of only one or two flavors. But the sparseness of such a factor would be too high for the specified constraint.

Summary

In multivariate statistics, the best algorithms are often those that let the data decide the features. By adding a sparseness constraint to the NMF algorithm, we let the tail wag the dog. The sparseness essentially dictates how many variables are important in each factor. If the sparseness value is high, then every factor will contain a small number of important variables, which might not be appropriate for the data.

There is also the mathematical fact that a constrained solution will always have a poorer fit than an unconstrained solution. By itself, that doesn't bother me too much, but you should get something valuable for the loss of fit. In a sparseness-constrained NMF, it isn't clear that the benefits offset the poorer approximation to the data.

Appendix: Sparseness and goodness of fit

It is an interesting exercise to plot the fit of rank-4 NMF approximation as you let the sparseness constraint vary from small to large values. The following image summarizes that experiment for the Whisky data. I have overlaid a few plots of the H` matrix so that you can see the effect that sparseness has on the NMF factors. The horizontal line is the relative error of the unconstrained NMF. The relative error compares the NMF error E=||X - W*H||_F with the smallest possible error, which is produced by using a rank-k SVD.

The post Nonnegative matrix factorization with sparseness constraints appeared first on The DO Loop.

Sparseness versus sparsity

First, I want to clarify the difference between Hoyer's definition of sparseness (for a vector) and the traditional notion of sparsity (for a matrix).

In numerical analysis, a matrix is said to be sparse if most of its elements are exactly zero and only a few elements are nonzero. For example, a large diagonal (or tridiagonal matrix) is sparse because the only nonzero elements are on (or near) the diagonal, and other elements are exactly zero. The sparsity of a matrix is reported as the proportion of elements that are zero. For an n x m, this is card(A[i,j] = 0) / (n*m), which is a number in [0,1]. Very sparse matrices have sparsity close to 1, whereas very dense matrices have sparsity close to 0. An n x n diagonal matrix has sparsity 1 – 1/n.

Hoyer's measure of sparseness is for vectors, not matrices. It does not count zeros and non-zeros. Rather, it looks at the magnitudes of the elements and measures how many are relatively large and how many are relatively small. If all elements are approximately the same magnitude, the vector has a sparseness measure that is close to 0. If there are a small number of large elements and the other elements are small, the vector has a measure that is close to 1.

Hoyer's sparseness measure

Hoyer's sparseness measure is defined as follows. If v ≠ 0 is an n-dimensional vector, n ≥ 2, the sparseness of v is defined as
      $S(v) = \frac{\sqrt{n} - \|v\|_1 / \|v\|_2}{\sqrt{n} - 1}$
where $\|v\|_1 = \sum_{i=1}^n |v_i|$ is the L₁ vector norm, and $\|v\|_2 = \sqrt{ \sum_{i=1}^n v_i^2 }$ is the L₂ or Euclidean norm.

Clearly, the ratio of the L₁ norm to the L₂ norm is an important quantity. There are a few easy mathematical results about the ratio $\|v\|_1 / \|v\|_2$:

1. The ratio is invariant under a change of signs for any component. That is, if S is a diagonal matrix whose diagonal elements are ±1, then $\|S v\| = \|v\|$.
2. The ratio is scale invariant. For any c ≠ 0, scalar multiplication by c does not change the ratio: $\|c v\|_1 / \|c v\|_2 = \|v\|_1 / \|v\|_2$.
3. The ratio is permutation invariant because the norms are permutation invariant. That is, if you permute the coordinates of v, you do not change the norm or the ratio of norms. In terms of matrices, if P is any permutation matrix, $\|P v\| = \|v\|$.

The first property tell you that it suffices to understand the ratio of norms for "positive vectors" in the first octant (that is, all components are positive). The second property tells you that you can standardize the vectors to have unit L₂ norm. The third property tells you that you can assume that |v₁| ≥ |v₂| ≥ … ≥ |v₂|, if you want to. These properties are apparent in the subsequent visualizations.

The range of the Hoyer measure

For any nonzero n-dimensional vector, the Hoyer measure of sparseness is in the interval [0, 1].

The Hoyer measure achieves its maximum value of 1 when (n-1) components of v are zero and the remaining component is nonzero. The properties in the previous section enable us to assume, without loss of generality, that the first component is nonzero. Then, ||v||₁ = ||v||₂, so the ratio of norms is 1. It follows that S(v) = 1 when v has only one non-zero component.

The Hoyer measure achieves its minimum value of 0 when all components of v are equal. The properties in the previous section enable us to standardize the vector so that its Euclidean norm is 1. If all components are equal, then ||v||₂ = sqrt(Σ a²) = 1 for some value, a. Thus, a = ±1/sqrt(n). Consequently, ||v||₁ = n / sqrt(n) = sqrt(n), so the numerator of Hoyer's measure vanishes. It follows that S(v) = 0 when all components of v are equal.

Examples of Hoyer's sparseness

Let's illustrate this concept with an example. Consider the following 6-D vectors. To within four decimal places, the Euclidean (or L2) norm of each vector is 1.

• q = {0.4082, 0.4082, 0.4082, 0.4082, 0.4082, 0.4082}. The L1 norm is 2.4492, and the Hoyer sparseness is 0.
• s = {0.5774, 0.5774, 0.5774, 0, 0, 0}. The L1 norm is 1.7322, and the Hoyer sparseness is 0.4950.
• u = {0.9948, 0.08, 0.04, 0.04, 0.02, 0.02}. The L1 norm is 1.1948, and the Hoyer sparseness is 0.8656. Notice that the sparseness is high even though no element is exactly 0.

The panel of graphs to the right visualizes the meaning of Hoyer's sparseness measure. Vectors that have a few elements that are larger than the other have the most sparseness. Vectors whose components are approximately uniformly distributed have the least sparseness.

Visualize Hoyer's sparseness for 2-D vectors

Since Hoyer's measure is invariant under scaling, it suffices to visualize the measure for unit vectors. If the vectors are two-dimensional, the unit vectors are on the unit circle centered at the origin. From the previous section, we know that the measure will obtain its maximum value of 1 for vectors that are aligned with a coordinate axis. It will obtain its minimum value of 0 when the components are equal in magnitude. In the first quadrant, this occurs when the vector makes an angle of 45 degrees or π/4 radians. By the permutation invariance, there are also minima for vectors at 3π/4, 5π/4, and 7π/4 radians.

The following SAS DATA step visualizes the Hoyer sparseness measure for 2-D vectors on the unit circle:

/* A measure of sparsity of a vector from:
   Patrick O. Hoyer. Non-negative matrix factorization with sparseness constraints. 
                     Journal of Machine Learning Research, 5:1457–1469, 2004.
   The "most sparse" vector is one for which only a single component is non-zero.
   That is, it aligns with a coordinate axis. These vectors have a sparseness measure of 1.
   The "least sparse" vector is one for which elements are equal. These vectors have a sparseness of 0.
*/
/* first, visualize for v on the unit circle in 2-D */
data circle_Hoyer;
   sn = sqrt(2);
   pi = constant('pi');
   do theta = 0 to 2*pi by pi/100;
      x = cos(theta);
      y = sin(theta);
      L1 = lpnorm(1, x, y);  /* L1 norm of v=(x,y) */
      Sparseness = (sn - L1) / (sn-1);
      output;
   end;
run;
 
/* use Spectral color ramp from ColorBrewer */
%let spectral = CX3288BD CX99D594 CXE6F598 CXFFFFBF CXFEE08B CXFC8D59 CXD53E4F;
/* Visualize Hoyer's sparseness measure on the unit circle */
title "Sparseness on the Unit Circle";
proc sgplot data=circle_Hoyer aspect=1;
   scatter x=x y=y / colorresponse=Sparseness colormodel=(&spectral) markerattrs=(symbol=CircleFilled);
   gradlegend / title="L1 Norm";
   xaxis grid; yaxis grid;
run;

As expected, the visualization shows that the Hoyer sparseness measure is 0 (blue) for the vectors (±1/√2, ±1/√2), which have components that are equal in magnitude. The measure is 1 (red) for vectors on a coordinate axis. Notice that I used the LPNORM function in Base SAS to compute the L1 norm of the vectors.

Visualize Hoyer's sparseness for 3-D vectors

In a similar way, you can visualize the sparseness of 3-D unit vectors. The measure will obtain its maximum value of 1 for vectors that are aligned with a coordinate axis. It will obtain its minimum value of 0 when the vectors components are equal in magnitude: (±1/√3, ±1/√3, ±1/√3).

The following SAS DATA step visualizes the Hoyer sparseness measure for 3-D vectors on the upper hemisphere of the unit sphere:

/* Generate the grid for the upper hemisphere */
data sphere_L1;
   sn = sqrt(3);
   delta = 2**(-5);     /* Set the resolution of the grid */
   /* Loop over (x,y) in unit circle, subset of [-1,1]x[-1,1] */
   do x = -1 to 1 by delta;
      do y = -1 to 1 by delta;
         if (x**2 + y**2) <= 1 then do; /* restrict to unit circle */
            z = sqrt(1 - x**2 - y**2);  /* upper hemisphere */
            L1 = lpnorm(1, x, y, z);    /* L1 norm of v=(x,y,z) */
            Sparseness = (sn - L1) / (sn-1);
            output;
         end;
      end;
   end;
   drop delta;
run;
 
/* use Spectral color ramp from ColorBrewer */
%let spectral = CX3288BD CX99D594 CXE6F598 CXFFFFBF CXFEE08B CXFC8D59 CXD53E4F;
/* Visualize the L1 norm as a heatmap */
title "Heatmap of Sparseness on the Upper Hemisphere of the Unit Sphere";
proc sgplot data=sphere_L1 aspect=1;
   heatmapparm x=x y=y colorresponse=Sparseness / colormodel=(&spectral);        
   gradlegend / title="L1 Norm";
run;

Because of the symmetries in the Hoyer sparseness measure, it suffices to visualize the measure on vectors in first octant where all elements are nonnegative. In the following image, I used a DROPLINE statement to locate the vector that has minimal Hoyer sparseness, which is (1/sqrt(3), 1/sqrt(3), 1/sqrt(3)).

A SAS function to compute Hoyer's sparseness

In SAS, the IML language is the natural place to implement a vector measure because the language supports vector operations in a natural way. The following PROC IML statements define a function that computes Hoyer's sparseness for a vector:

proc iml;
/* Hoyer's sparseness measure.
   For a column vector, v, this function returns 
   the Hoyer (2004) measure of sparseness
   s = (sqrt(n) - ||v||_1 / ||v||_2) / (sqrt(n)-1)
*/
start Sparseness(v);
   u = colvec(v);    /* always treat input as a column vector */
   L2 = norm(u, "L2");
   if L2=0 then return(.);
   L1 = norm(u, "L1");
   sn = sqrt( countn(u) );
   sparseness = (sn - L1/L2) / (sn - 1);
   return( clip01(sparseness) );
finish;
/* Helper function: clip values into [0,1]
   If x is any matrix, return a new matrix whose i_th element has the value 
   0    if x[i] < 0
   x[i] if 0 <= x[i] <= 1
   1    if x[i] > 1
   See https://blogs.sas.com/content/iml/2026/02/04/clip-values.html */
start clip01(x);
   return( 0 <> (x >< 1) );
finish;
 
/* Examples: Compute the Hoyer measure for these columns */
/*   q       r  s       t    u       v */  
A = {0.4082  1  0.5774  3.3  0.9948  1,
     0.4082  1  0.5774  1    0.08    0,
     0.4082  1  0.5774  1    0.04    0,
     0.4082  1  0       0.3  0.04    0,
     0.4082  1  0       0.2  0.02    0,
     0.4082  1  0       0.2  0.02    0};
 
S = j(1, ncol(A), .);
do i = 1 to ncol(A);
   S[i] = Sparseness( A[,i] );
end;
print S[F=6.4 c=('q':'w') L="Hoyer's Sparseness"];

The table shows the sparseness for each column of the A matrix.

• The first and second columns (q and r) both have sparseness 0 because both columns are constant in value. The first column has unit L2 norm whereas the second column does not. But it doesn't matter, because the Hoyer measure is normalized by using the ratio ||v||₁ / ||v||₂.
• The third column (s) has three zero elements and three equal and nonzero elements. The sparseness measure is 0.495.
• The fourth column (t) has one large element, two medium-sized elements, and three small but nonzero elements. The sparseness measure is 0.5445.
• The fifth column (u) has one large element and five small but nonzero elements. The sparseness measure is 0.8656.
• The last column (v) has one large element and five elements that are zero. The sparseness measure is 1, which is the largest possible value of sparseness.

An earlier section visualized three of these columns in a panel of needle plots.

Summary

This article introduces the Hoyer sparseness measure, S(v), for nonzero vectors v. The measure indicates whether a vector has many components that are the same magnitude (low sparseness) as opposed to vectors that have many small components and only a few components that are large in magnitude. The sparseness measure is invariant to scaling, permutation of components, and switching the signs of components. Accordingly, you can visualize the measure by using vectors whose components are all positive and whose Euclidean length is 1.

A future article shows an application of the Hoyer measure to the nonnegative matrix factorization (NMF) method.

The post Hoyer's sparseness measure appeared first on The DO Loop.

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.