It is straightforward to modify any correlation matrix to change the number of positive and negative values in the off-diagonal locations. This article shows the mathematics and a simple SAS IML function. The technique enables you to create a new correlation matrix from one that you already have. The entries in the new correlation matrix have the same magnitude, but different signs.

A simple transformation that changes correlation

Let X and Y be data vectors. (Or, they could be random variables if you prefer a more rigorous mathematical treatment.) If the correlation ρ = Corr(X,Y) is not zero, then Corr(-X,Y) has the opposite sign from ρ. In other words, multiplying a variable by -1 changes the sign of the correlation with other variables in the data set.

Of course, you could also change the sign of Y. If R is a 2x2 correlation matrix with off-diagonal element ρ, then you can construct a new correlation matrix by using the similarity transformation Q = S*R*S`, where S is a diagonal 2x2 "sign matrix." The S matrix contains ±1 on the diagonal. If S has one -1 on the diagonal, then Q has -ρ on the off-diagonal. Otherwise, Q has +ρ for the off-diagonal element.

You can extend this result to more variables. If R is a dxd correlation matrix, and S is any diagonal sign matrix with values ±1, then Q = S*R*S` is also a correlation matrix. The (i,j)th element of Q is Q[i,j] = S[i]*R[i,j]*S[j]= ±R[i,j]. Mathematically, Sylvester's Theorem ensures that Q is a valid correlation matrix. If R is symmetric and positive semi-definite, then so is Q.

Of course, if S is a symmetric matrix, then the similarity transformation simplifies to S*R*S, which is the form that I use in the next section.

A SAS IML program that changes correlation

The following program defines a function that computes a similarity matrix. Assume that R is a dxd matrix and v is a dx1 vector. The function returns the matrix S*R*S` where S = diag(v). I immediately call this function for the special case where R is a correlation matrix, and v[i]= ±1.

I use a programming technique that is worth noticing. You should never multiply by a large diagonal matrix. Instead of diag(v)*R*diag(v), you can use the more efficient computation v#R#v`, where '#' indicates elementwise multiplication. This trick speeds up the computation and uses much less memory.

proc iml;
/* Given a dxd matrix, R, and a dx1 vector v, return the similar matrix
   Q = diag(v)*R*diag(v).  
   An important application is to let v be a vector of +/-1. Then, 
   if R is a valid correlation matrix, Q is also a correlation matrix
   where some correlations have changed signs.
*/
start SimilarMat(R, v);
   S = colvec(v);
   Q = S # R # S`;
   return(Q);
finish;
 
/* --- Example Usage --- */
/* Create AR(1) correlation matrix:
   https://blogs.sas.com/content/iml/2012/11/05/constructing-common-covariance-structures.html */
start AR1Corr(dim, rho);
   u = cuprod(j(1,dim-1,rho)); /* cumulative product */
   return( toeplitz(1 || u) );
finish;
 
/* Create a 5x5 AR(1; 0.5) matrix */
R = AR1Corr(5, 0.5);
v = {1, -1, -1, 1, -1};  /* change correlations for the 2nd, 3rd, and 5th variables */
Q = SimilarMat(R, v);
 
origNames = 'X1':'X5';
newNames  = {'X1' '-X2' '-X3' 'X4' '-X5'};
print R[c=origNames r=origNames F=BestD7.], 
      Q[c=newNames r=newNames F=BestD7.];

The output shows that the original matrix, R, has all positive correlations. However, in the new matrix, the correlation in the (i,j)th cell has the sign v[i]*v[j]. This creates a correlation matrix that has negative values.

Summary

This article shows a simple trick for changing the signs of elements in a correlation matrix. If R is any correlation matrix and S = diag(v) is a diagonal sign matrix with values ±1, then the product S*R*S` is a correlation matrix that has different signs than R. You can compute the product efficiently without ever forming a diagonal matrix. This trick is useful in simulation studies and in creating matrices to test statistical procedures.

The post An easy way to create a correlation matrix with negative correlations appeared first on The DO Loop.

Correlation matrices that have a common spectrum

A matrix is positive semi-definite if all eigenvalues are greater than or equal to zero. All correlation matrices are symmetric and positive semi-definite.

The spectrum of the matrix is the set of eigenvalues: Λ = {λ₁, λ₂, ..., λ_d}. For a positive definite matrix, λ_i > 0 for all i. For a semi-definite matrix, λ_i ≥ 0. For any square matrix, the sum of the eigenvalues equals the trace of the matrix. Consequently, for an d x d correlation matrix, Σ_i λ_i = d.

As shown in the previous article, there are many different correlation matrices that have the same spectrum.

The Marsaglia and Olkin method for random correlation matrices

The Marsaglia-Olkin (1984) paper is behind a paywall, so I have not read it. Fortunately, Gentle (2003, p. 200) presents pseudo-code for their algorithm. Gentle's pseudo-code contains two typos (Steps 5 and 6), which I have corrected in the following description.

The input to the Marsaglia-Olkin routine is Lambda, which is the spectrum of a d x d SPD correlation matrix. Lambda is a d x 1 vector {λ₁, λ₁, ..., λ_d}, where λ_i > 0 and Σ_i λ_i = d. The algorithm outputs R, which is a random correlation matrix such that the eigenvalues of R are the specified spectrum. The algorithm is as follows:

0. Set E = I_d and k=1.
1. Generate a d x 1 vector, w, of i.i.d. standard normal random variates.
   Form x = E*w and compute the scalar a = (1-Lambda)` * x##2
2. Generate a d x 1 vector, z, of i.i.d. standard normal random variates.
   Form y = E*z and compute the scalars:
   • b = (1-Lambda)` * (x#y)
   • c = (1-Lambda)` * y##2
   • e² = b##2 - a*c
3. If e² < 0, goto Step 2.
4. Choose a random sign, s ∈ {-1, 1}. Set r = (b+s*e)/a * x - y.
5. Choose another random sign, s ∈ {-1, 1}. Set p_k = s*r/norm(r). (Fixes typo in Gentle.) This vector will be a row in the P matrix. It is a unit vector in the hypercone defined by the equation p_k`(I - *Lambda;)p_k = 0.
6. Perform a rank-one update: E → E - p_k`*p_k (fixes typo), and increment k=k+1. At each step, E is the projection matrix onto the subspace orthogonal to the rows of P that have been constructed so far.
7. If k < d, goto Step 1.
8. Generate a d x 1 vector, w, of i.i.d. standard normal random variates.
   Form x = E*w and set p_d = x / norm(x).
9. Construct the matrix P by using the vectors p_k as its rows.
   Return P*diag(Lambda)*P` as the random correlation matrix.

I always say that the main benefit of a matrix language such as MATLAB, R, and SAS IML is the ease of converting a high-level pseudo-code description of an algorithm into a ready-to-run program. The Appendix to this article contains my implementation of the algorithm in a SAS IML moduled called CorrWithEigen_MO. The next section shows how to use the function to generate a random correlation matrix.

Generate random correlation matrices in SAS

First, run the program in the Appendix, which stores the CorrWithEigen_MO function. You can then load the function and use it in other programs. The following call to PROC IML generates two random 6 x 6 correlation matrices that have the eigenvalues Lambda = {1.8, 1.5, 1, 1, 0.5, 0.2}.

proc iml;
load module=(CorrWithEigen_MO);
call randseed(123);
reset fuzz;   /* print tiny numbers as 0 */
 
/* Specify the spectrum for a 6x6 corr matrix */
Lambda = {1.8, 1.5, 1, 1, 0.5, 0.2};
/* get random 6x6 correlation matrix with this spectrum */
Corr1 = CorrWithEigen_MO( Lambda );
/* call again, make sure we get a different answer */
Corr2 = CorrWithEigen_MO( Lambda );
print Corr1[F=best6.], Corr2[F=best6.];
 
/* verify the spectrum of the random matrices */
v1 = eigval(Corr1);
v2 = eigval(Corr2);
print Lambda v1 v2;

The output shows two different random 6x6 correlation matrices. Each matrix has the same set of eigenvalues, as demonstrated by the calls to the EIGVAL function.

Summary

This article provides an IML implementation of the Marsaglia and Olkin (1984) algorithm for the direct generation of random correlation matrices that have a specified spectrum. You can use the algorithm to generate random correlation or covariance matrices in SAS. You can use this function in simulation studies where you need a positive definite matrix.

Appendix: The Marsaglia and Olkin Algorithm in IML

This section defines an IML function (and some helper functions) that implement the Marsaglia and Olkin (1984) algorithm for the direct generation of random correlation matrices that have a specified spectrum.

/* From Gentle (2003, p. 200) Random Number Generation and Monte Carlo Methods
   Algorithm 5.9 
   Marsaglia-Olkin (1984) Method for Random Correlation Matrices with Given Eigenvalues
   Ref: Marsaglia, G. and Olkin, I. (1984) "Generating Correlation Matrices", 
        SIAM J. on Sci. and Stat. Comp., 5(2), pp 470-475. doi 10.1137/0905034.
 
   INPUT:
   Lambda = the spectrum of a dxd SPD correlation matrix. Lambda is a vector
            (lambda_1, ..., lambda_d), where lambda_i > 0 and \sum lambda_i = d.
   OUTPUT:
   R = dxd correlation matrix, R, such that the eigenvalues of R are the specified spectrum:
        Lambda(R) = (lambda_1, ..., lambda_d)
   ALGORITHM:
   0. Set E = I_d and k=1.
   1. Generate d-vector, w, of i.i.d. std normal deviate.
      Form x = E*w and compute the scalar a = (1-Lambda)` * x##2
   2. Generate a d-vector, z, of i.i.d. std normal variates. 
      Form y = E*z and compute the scalars 
      b = (1-Lambda)` * (x#y)
      c = (1-Lambda)` * y##2
      e^2 = b##2 - a*c
   3. If e^2 < 0, goto Step 2.
   4. Choose a random sign, s \in {-1, 1}. Set r = (b+s*e)/a * x - y.
   5. Choose another random sign, s \in {-1, 1}. Set p_k = s*r/norm(r). (fixes typo in Gentle)
   6. Perform a rank-one update: E -> E - p_k`*p_k (fixes typo), and increment k=k+1.
   7. If k < d, goto Step 1.
   8. Generate a d-vector, w, of i.i.d. std normal deviate.
      Form x = E*w and set p_d = x / norm(x).
   9. Construct the matrix P by using the vectors p_k as its rows.
      Return P*diag(Lambda)*P` as the random correlation matrix.
*/
 
proc iml;
/* return k random signs with values +1 or -1.
   Let B be a binary random variable chosen uniformly at random from {0,1}.
   Then C = 2*(B - 1/2) is random uniform in {-1,1}.
*/
start RandSign(k);
   return( 2*(randfun(k, "Bernoulli", 0.5) - 0.5) );
finish RandSign;
 
/* for column vector v, standardize so sum(v)=nrow(v) */
start StdizeEigenval(v);
    L = colvec(v);
    return ( L * nrow(L)/sum(L) );  /* make sum(v)=dimension */
finish StdizeEigenval;
 
/* The  Marsaglia-Olkin (1984) method for creating a random correlation 
   matrix with a specified set of eigenvalues */
start CorrWithEigen_MO(targetLambda, maxIters=100);
   lambda = StdizeEigenval(targetLambda);
   d = nrow(Lambda);
   if any(Lambda < 0) then do; 
      print "ERROR: The spectrum must contain only nonnegative values"; 
      return( J(d,d,.) ); 
   end;
   /* If the spectrum is {1,1,...,1}, the only solution is I(d).
      The do-while logic in the M-O method will enter an infinite loop
      in this case, so detect this edge case and return the identity matrix */
   if norm(1 - Lambda, "LInf") < 1e-8 then return (I(d));
 
   /* use column vectors for Lambda, w, and z */
   w = j(d, 1, .);
   z = j(d, 1, .);
   P = j(d, d, .);
   /* 0. Set E = I_d and k=1. */
   E = I(d);
 
   /* compute the k_th row of P */
   do k = 1 to d-1;
      /* 1. Generate dx1 vector, w, of i.i.d. std normal variates.
            Form x = E*w and compute the scalar a = (1-Lambda)` * x##2. */
      call randgen(w, "Normal");
      x = E*w;
      a = sum((1-Lambda) # x##2);
      /* 2. Generate a dx1 vector, z, of i.i.d. std normal variates. 
            Form y = E*z and compute the scalars 
            b = (1-Lambda)` * (x#y)
            c = (1-Lambda)` * y##2
            e2 = b##2 - a*c
         3. If e2 < 0, goto Step 2.     */
      e2 = -1;
      do cnt = 1 to maxIters until(e2 > 0);
         call randgen(z, "Normal");
         y = E*z;
         b = sum((1-Lambda) # x # y);
         c = sum((1-Lambda) # y##2);
         e2 = b##2 - a*c;
      end;
      if e2 <= 0 then do;
         print "ERROR:  Marsaglia-Olkin method failed after 100 iterations.";
         return (J(d, d, .)); 
      end;
      /* 4. Choose two random signs, s \in {-1, 1} */
      s = RandSign(2);
      r = (b+s[1]*sqrt(e2))/a * x - y;   /* Set r = (b+s*e)/a * x - y (Note Pr(a=0)=0) */
      p_k = s[2] * r/norm(r);            /* 5. Set p_k = s/norm(r) * w */
 
      /* 6. Perform a rank-one update: E -> E - r*r` */
      E = E - p_k*p_k`;
      P[,k] = p_k;
   end;   /* 7. If k < d, goto Step 1. */
 
   /* 8. Generate a dx1 vector, w, of i.i.d. std normal variates.
         Form x = E*w and set p_d = x / norm(x). */
   call randgen(w, "Normal");
   x = E*w;
   P[,d] = x / norm(x);
   /* 9. Construct the matrix P by using the vectors p_k as its rows.
         Return P*diag(Lambda)*P` as the random correlation matrix.
         (Note: Instead, I will construct P by using p_k as columns, then return P`*Lambda*P) */
   G = P` # sqrt(Lambda`);
   Corr = G*G`;
   return( Corr );
finish CorrWithEigen_MO;
store module=(RandSign StdizeEigenval CorrWithEigen_MO);
QUIT;

The post A direct method to generate correlation matrices with specified eigenvalues appeared first on The DO Loop.

My implementation is based on a numerical technique called the Genz transformation (A. Genz, 1992; Genz and Bretz, 2009, pp. 49-50). It leverages the geometric fact that the Cholesky factor of a correlation matrix provides a transformation that maps uncorrelated MVN variables to correlated MVN variables. I previously wrote about using the Cholesky transformation to correlate variables, and Genz's method is a clever and useful application of this linear algebraic technique. The technique goes back to multivariate statistical methods developed by Mahalanobis in the 1930s.

This article discusses the Genz transformation method and shows how to compute and visualize the integrand for a bivariate normal distribution. This low-dimensional example demonstrates the main ideas that are used in higher dimensions.

The computation of MVN probabilities

The goal is to estimate the probability that a random multivariate normal vector falls inside a specific rectangular region:
     $P( L_1 < X_1 < U_1, L_2 < X_2 < U_2, \dots, L_d < X_d < U_d )$
where (X₁, X₂, ..., X_d) ~ MVN(0, Σ), and Σ is the covariance matrix. The lower and upper limits are defined by the vectors L = (L₁, L₂, ..., L_d) and U = (U₁, U₂, ..., U_d). We'll adopt the convention that each L_i is less than U_i and we'll allow elements of L to represent negative infinity and elements of U to represent positive infinity. In SAS, you can use the special missing value .M to represent minus infinity and .I to represent positive infinity.

The probability is formally computed by evaluating a multidimensional integral:
$\int_{L_1}^{U_1} \int_{L_2}^{U_2} \cdots \int_{L_d}^{U_d} \phi(x; \Sigma)\, dx$ where φ(x; Σ) is the MVN density function.

Does that integral look scary to you? Because it looks scary to me! In numerical analysis, there are not many good ways to solve high-dimensional integrals. For low to moderate dimensions, you can use ordinary Monte Carlo simulations to estimate this probability. There are two common Monte Carlo techniques:

• The Region method (a.k.a., "naive Monte Carlo"): Generate B random variates X_i ~ MVN(0, Σ) and calculate the proportion of the X_i inside the rectangular region defined by the limits. For completeness, the Appendix shows a traditional Monte Carlo implementation of the region method.
• The Average Function method: Evaluate a certain function for B random variates in the domain of the function. Compute the average of the values. I am intentionally leaving "the function" and "the domain" unspecified for now. As we will soon see, Genz transforms the domain of the problem, and the function that gets evaluated is somewhat complicated.

For the Average Function method, you can generate quasi-random variates and perform a quasi-Monte Carlo analysis. A quasi-Monte Carlo estimate has a smaller standard error than a naive Monte Carlo estimate that uses the same number of points.

Unfortunately, QMC methods are limited to finite regions because quasi-random sequences are generated inside the unit hypercube. Genz solved this issue by transforming the general MVN region into an equivalent integration problem strictly on the unit hypercube. Because the innermost integral can be evaluated analytically, a d-dimensional probability reduces to an integral over the (d-1)-dimensional hypercube:
     $\int_{(0,1)^{d-1}} g(w) \, dw$
where the integrand, g, depends on the covariance matrix and the limits of integration.

You can use QMC techniques to efficiently estimate the integral with high accuracy. You generate N quasi-random points in the (d-1)-dimensional hypercube, evaluate the function g at each point, and take the average, as follows:
     $p \approx \frac{1}{N} \sum_{i=1}^N g(w_i)$

Of course, there is no such thing as a free lunch, so the Genz integration function, g, is more complicated than the standard MVN density function. But Genz's paper was seminal because it shows how to evaluate the integral by using linear transformations and the inverse CDF transformation.

An overview of the Genz transformation

Genz's method consists of three main steps:

1. The Cholesky transformation: Transform the problem from correlated variables into independent standard normal variables. This transformation has been known since the 1930s.
2. Conditional probability: When integrating the i_th variable, you must condition the bounds on the exact values drawn for the previous (i-1) dimensions. This step uses conditional probability to update the lower and upper limits of integration for each variable.
3. The inverse CDF Transformation: You can use the inverse CDF transformation to map these updated bounds into the unit hypercube.

A 2-D Implementation in SAS IML

Although Genz's algorithm works in arbitrary dimensions, let's use it to integrate the MVN probability in 2-D. SAS already has functions for evaluating the bivariate normal probability, so we can use those functions to check whether we implemented Genz's method correctly.

Assume L and U are row vectors. The following SAS IML program defines two helper functions to handle infinite bounds for the CDF("Normal") function. The main function is called GenzIntegrand_2D. As input parameters, it takes two-element vectors L and U and a 2x2 correlation matrix, R. It also takes a vector of values, w, in the interval (0,1). There are several ways you can use that input argument. I will use it to visualize the function by passing in a grid of values on the interval (0,1). However, if you want to average the function to estimate the 2-D MVN probability, you would send in a set of quasi-random numbers in (0,1) and average the output.

The Genz transformation "integrates out" the innermost variable in the problem. As a result, the Genz integrand is a function of one fewer variables than the original problem. For a 2-D probability, the Genz integrand is a 1-D function. For a 3-D probability, the integrand is a function of two variables, and so forth.

The following SAS IML program begins with three helper functions. The main function is GenzIntegrand_2D. The comments in the function indicate the three main steps for transforming the bivariate normal probability into a 1-D function that can be integrated.

proc iml;
/* Helper functions: Extend the CDF function so that 
             / 0 if x = -Infinity
   cdf(x) = {  cdf(x) if x is finite 
             \ 1 if x = +Infinity
*/
start cdf_at_lower(L);
   if L=. then return( 0 );    /* = cdf("Normal", -Infinity) */
   return( cdf("Normal", L) );
finish;
start cdf_at_upper(U);
   if U=. then return( 1 );    /* = cdf("Normal", +Infinity) */
   return( cdf("Normal", U) );
finish;
/* The Title2_from function sets the TITLE2 statement dynamically at runtime. See 
   https://blogs.sas.com/content/iml/2015/01/14/global-statements-loops.html
  The string looks like
  "rho=0.5; P(. < X1 < 2 & -2 < X2 < 1)" 
*/
%macro char(x);  choose(&x=., ".", char(&x))  %mend;
start Title2_from(L, U, rho);
   rhoStr = cats("rho=",char(rho),";");
   range1 = cats(%char(L[1]),"< X1 <",%char(U[1]));
   range2 = cats(%char(L[2]),"< X2 <",%char(U[2]));
   parmStr = catx(" ", rhoStr, "P(", range1, "&", range2, ")");
   call execute("title2 '" + parmStr + "';" );
finish;
 
/* The Genz method transforms a vector of independent standard normal variables 
   Y ~ MVN(0, I) into the observed correlated variables X~MVN(0, R)
   by using the Cholesky factor, C. In 2-D, relates X to Y via the linear equations
   X1 = C11*Y1
   X2 = C21*Y1 + C22*Y2
*/
start GenzIntegrand_2D(w, L, U, R);
   /* The ROOT function returns an upper triangular matrix. Transpose to get lower triangular */
   C = t(root(R)); 
 
   /* First dimension (outer integral) */
   v1 = 0;            /* (un)conditional mean of X1 */
   c11 = C[1,1];
 
   /* Standardize limits and calculate probability mass M1 */
   alpha = cdf_at_lower( (L[1]-v1)/c11 );
   beta  = cdf_at_upper( (U[1]-v1)/c11 );
   M1 = beta - alpha;   /* 1-D marginal probability */
 
   /* The inverse CDF maps each w in (0,1) to normal variable y1.
      Note: y1 is a vector in (-Infinity, Infinity) */
   y1 = quantile("Normal", alpha + w * M1);
 
   /* Second dimension (inner integral) */
   c21 = C[2,1];
   c22 = C[2,2];
 
   /* Conditional mean depends on the y1 from the previous step.
      Subtract the conditional mean and standardize limits for inner integral. */
   v2 = c21 * y1;
   alpha = cdf_at_lower( (L[2]-v2)/c22 );
   beta  = cdf_at_upper( (U[2]-v2)/c22 );
   M2 = beta - alpha;
 
   /* Final integrand: M1 is a scalar, M2 is a vector, so g is a vector. */
   g = M1 * M2;
   return (g);
finish;
 
/* Visualize the 1-D integrand for the following 2-D problem */
rho = 0.5;
R = (1   || rho) //
    (rho || 1  );
L = {.M  -1};
U = { 1   2};
 
/* Evaluate and plot the Genz integrand, which is defined on the open interval (0,1) */
w = {1E-4, 5E-4, 0.001, 0.005} // 
    T(do(0.01, 0.99, 0.01))    //
    {0.995, 0.999, 0.9995, 0.999}; 
g = GenzIntegrand_2D(w, L, U, R);
 
title  "Genz Integrand for 2-D MVN";
call Title2_from(L, U, rho);
xLabel = "w (Uniform QMC Coordinate)";
yLabel = "g(w) = Probability Mass";
call series(w, g) grid={X Y} label=xLabel
                  other=cat("yaxis grid min=0 label='",yLabel,"';");

The resulting curve represents the integrand for the transformed integration problem. To find the true MVN probability in the region, you integrate the function over the interval (0,1). You can estimate the integral by evaluating the function on a 1-D QMC sequence, then taking the average. Visually, the probability is the average height of this smooth positive function.

Since this is a bivariate problem, you can use SAS to find that the true probability is about 0.677. Notice that this value is the average height of the function on (0, 1).

The Genz integrand for uncorrelated variables

If the correlation matrix for a MVN distribution is the identity, the problem decomposes into a product of univariate marginal probabilities. The Genz integrand for this situation is a straight line. The height of the line is the probability value. To see this fact in action, let's use a very small value of rho so that the function will be almost a horizontal line.

/* If you use a tiny value for rho, the integrand approaches a constant 
   function. When rho=0 exactly, the integrand is the product of the 
   two marginal probabilities. The function becomes a constant.
*/
rho = 0.01;  /* prob = 0.6882697 for this value of rho */
R = (1   || rho) //
    (rho || 1  );
g = GenzIntegrand_2D(w, L, U, R);
call Title2_from(L, U, rho);
call series(w, g) grid={X Y} label=xLabel
                  other=cat("yaxis grid min=0 max=0.75 label='",yLabel,"';");

As expected, the Genz integrand is very close to a horizontal line. The height of the curve is very close to 0.688, which is the probability value for this problem. The curve becomes a horizontal line for rho=0.

The Genz integrand for small probabilities

One of the problems with a Monte Carlo simulation that uses the "Region method" is that the method requires many random variates to get an accurate estimate. For example, if you set rho=-0.5, the bivariate normal probability P( X1 < -1 & -3 < X2 < -2 ) = 0.000145. If you use the Monte Carlo "Region method" to estimate the probability, you would expect only 145 points per million to be inside the region of interest. Most variates are "wasted" because they are outside of the region. In contrast, the Genz transformation of this problem does not waste any evaluations. The Genz integrand is very short, but every value in (0,1) contributes to the average value of g on (0, 1). This is shown in the following program:

/* What does the integrand look like if we use a negative rho
   and specify a region in which the probability is small?
*/
rho = -0.5;
R = (1   || rho) //
    (rho || 1  );
L = {.M  -3};
U = {-1  -2};
g = GenzIntegrand_2D(w, L, U, R);
call Title2_from(L, U, rho);
call series(w, g) grid={X Y} label=xLabel
                  other=cat("yaxis grid min=0 label='",yLabel,"';");

As expected, the integrand is very small. The maximum value of the function on (0, 1) is approximately 0.0003. The function is approximately linear and g(0) ≈ 0, so a quick back-of-the-envelope calculation of the probability is 0.00015, which is pretty doggone close to the true value! You could estimate the probability by using only a few thousand quasi-random points in (0,1), which is a huge savings over the millions of points that are required for a naive Monte Carlo estimate.

Summary

This article presents a general computational framework for estimating a probability for a rectangular region of a multivariate normal distribution. Three bivariate examples are given: a problem where the variables are moderately correlated, a problem where the variables are almost independent, and an example where the probability is very small. For each case, I visualize the Genz integrand, which is a 1-D function. By integrating the Genz integrand, you obtain the desired probability. You can use quasi-Monte Carlo integration, which amounts to evaluating the integrand many times and taking the average.

The algorithm generalizes to higher dimensions. For a d-dimensional MVN distribution, the Genz integrand is a function of d-1 variables that each have (0, 1) as a domain. Thus, the problem reduces to estimating the integral of a function on the unit hypercube, which can be computed efficiently by using quasi-Monte Carlo methods.

Appendix: The traditional Monte Carlo estimate

As mentioned in the article, the Genz transformation is more efficient than the traditional Monte Carlo estimates. For example, the following PROC IML program simulates N=1000 points from the bivariate normal distribution with ρ = 0.5 that is shown in the first graph. For this run, using 1000 data points results in an estimated probability of 0.664. The true probability of 0.677, so the estimate is not bad. However, the standard error of the Monte Carlo method is proportional to 1/sqrt(N), whereas the Genz transformation and quasi-random integration has a much smaller standard error proportional to 1/N.

proc iml;
rho = 0.5;
R = (1   || rho) //
    (rho || 1  );
L = {.M  -1};
U = { 1   2};
 
/* Visualize the Monte Carlo estimate P(X1< 1 && -1 < X2 < 2) */
N = 1000;
call randseed(123, 1);
X = randnormal(N, {0 0}, R);
hit = (X[,1] < U[1]) && (L[2] < X[,2] && X[,2] < U[2]);
p_est = mean(hit);
print p_est;
 
result = X||hit;
create MCViz from result[c={'x1' 'x2' 'inRegion'}];
append from result;
close;
quit;
 
title "Area Method Monte Carlo Computation";
title2 "Estimated Probability = 0.664; N = 1000";
proc sgplot data=MCViz;
   scatter x=x1 y=x2 / group=inRegion;
   refline 1 / axis=x;
   refline -1 2 / axis=y;
run;

The post A visual introduction to the Genz method for computing multivariate normal probabilities appeared first on The DO Loop.

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.