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

The Childs-Sun formula for dimension 4

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

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

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

Conditional probability and the geometry of the formula

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

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

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

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

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

Implementing the P4 formula in SAS IML

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

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

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

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

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

Summary

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

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

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

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

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

Appendix 2: The sum notation

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

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

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

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

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

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

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

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

Orthant probabilities

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

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

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

Orthant probabilities in dimensions 1 and 2

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

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

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

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

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

Orthant probabilities in dimension 3

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

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

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

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

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

The computation for other 2-D quadrants

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

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

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

The computation for other k-D orthants

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

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

Summary

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

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

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

Appendix: How the correlation matrix changes under coordinate reflections

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

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

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

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

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

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

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

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

The 1-D Z-score transformation

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

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

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

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

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

Standardizing the multivariate normal distribution

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

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

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

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

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

Converting from covariance to correlation in SAS IML

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

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

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

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

An Example: The 5-D Min Matrix

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

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

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

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

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

Summary

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

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

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

The EFFECT statement provides a convenient way to specify knots and splines from within a SAS regression procedure. However, certain applications (such as simulation studies) require more control of knot placement. For these applications, SAS IML provides the BSPLINE function, which enables you to evaluate B-spline basis functions at a set of X values, given the knots. This article demonstrates the BSPLINE function and shows how to generate knots according to the common schemes that are supported by the EFFECT statement.

The example data

Unlike other statistical modeling methods, knots and splines can be defined without using any data. If you specify an interval, [a, b], and a set of knots, you can construct a B-spline basis on the interval. However, in practice, the interval is [min(X), max(X)] for some data variable, X. Furthermore, a common method to place the knots inside the interval is to use percentiles of the data. To accommodate these two facts, let's use the Horsepower variable in the Sashelp.cars data set as "X".

So that you can easily run the examples on your own data, I have put the name of the data set and the name of the variable into macro variables (DSName and XVar, respectively). By modifying those two macro variables, the code in this article should work for any data set and any numerical variable.

/* example data */
data Have;
set sashelp.cars;
keep Horsepower;
run;
 
/* define these macro variables to point to your data set and numerical variable */
%let DSName = Have;
%let XVar = Horsepower;

The BSPLINE function

PROC IML supports the BSPLINE function, which enables you to evaluate a B-spline basis function on an arbitrary set of X values.

First, let me demonstrate that you don't need any data to generate a B-spline basis. For the Horsepower variable, X, min(X) = 73, max(X) = 500, and the 25th, 50th, and 75th percentiles of the data are 165, 210, 255. By using only these numbers, you can construct a B-spline basis on the interval [73, 500]. The following IML statements place internal knots at the values 165, 210, 255. There are two "clamped" boundary knots at the endpoints of the interval.

/* Use the BSPLINE function in PROC IML to reproduce an example from 
   https://blogs.sas.com/content/iml/2026/04/22/viz-knots-splines.html
   which uses X=Horsepower variable from Sashelp.Cars
*/
proc iml;
deg = 2;
/*       |-min(X)-|-25th 50th 75th Pctl -|-max(X)-|  */
knots = {73, 73,     165, 210, 255,       500, 500};
t = T(do(73, 500, (500-73)/99));   /* evaluate at 100 points in range of X */
basis = bspline(t, deg, knots);    /* 100x6 matrix. Columns are spline functions */
 
/* columns of BASIS are the splines evaluated at t */
title "Second Spline Basis";
call series( t, basis[,2] ) grid={x y} other="refline 165 210 255/axis=x;";

The BASIS variable is a matrix that has 100 rows and six columns. Each column is the result of evaluating a B-spline basis function at a set of evenly spaced values in the interval [73, 500]. The second basis function has support (that is, nonzero values) on the interval [73, 210]. The graph of that basis function is shown above.

The BASIS matrix contains the spline values in "wide form." If you want to visualize the basis by plotting all six functions on the same graph, you can use the WIDETOLONG function to reshape the data, then use the GROUP= option on the SERIES call to plot the functions, as follows:

/* if you want to plot all columns, convert wide to long and use the GROUP= option 
   See https://blogs.sas.com/content/iml/2015/02/27/multiple-series-iml.html
*/
run WideToLong(tt, Value, Spline_Number,              /* output (long) */
               basis, t, "Spl1":"Spl6");      /* input (Basis is wide) */
 
title "All B-Splines (deg=2)";
call series(tt, Value) group=Spline_Number grid={X Y} 
                       other="refline 165 210 255/axis=x;" label={"&XVar" "Spline Value"};

This graph shows the same B-spline basis as for the PERCENTILES(3) example in my previous article. The main difference is that the previous graph shows the splines evaluated at the observed values of X, whereas the graph in this article evaluates the splines on an evenly spaced grid of points. The second basis function is now displayed in a red color.

A SAS IML function for common knot-placement schemes

Although B-splines can be defined without any reference to data, in practice they are often used as regression splines. As such, the interval on which they are defined corresponds to the range of the X data. For some knot-placement schemes, the internal knots are placed at percentiles of the X variable. Furthermore, to create a design matrix for fitting a regression model, you must evaluate the splines at the observed values of X. So, in practice, it is useful to extract the knot positions from data.

As discussed in a previous article, the EFFECT statement in SAS regression procedures handles these details automatically. It is useful to define a SAS IML function that can generate the same knot-placement schemes that the EFFECT statement uses. The schemes that depend on the data are shown in the following table:

It is straightforward to define a SAS IML function that returns the same set of knots as the EFFECT statement. (You can skip the "LIST" method, which assumes that the user already knows the positions of the knots.) The following function shows one implementation, which uses separate helper functions for generating internal and external knots. You can test the function by loading the Horsepower variable and sending that data vector into the function, along with a parameter that specifies the number and location of knots. So that the main ideas are evident, this implementation does not perform any argument validation.

/* reproduce the knot-placement methods in the EFFECT statement:
   Method   Parm        Equivalent KNOTMETHOD= option            External Knots
   ======   ========    ======================================   ==============
   "PCTL"   scalar n    PERCENTILES(n)                           Boundary
   "PCTL"   vector v    PERCENTILELIST(v), where 0 < v < 100     Boundary
   "EQUAL"  scalar n    EQUAL(n)                                 Evenly spaced
   "RANGE"  vector v    RANGEFRACTIONS(v), where 0 < v < 1       Boundary
*/
start Knots(x, deg, method, parm);
    int_knots = Knots_Internal(x, method, parm);
    if upcase(method)="EQUAL" then do;
       delta = int_knots[2] - int_knots[1];
       knots = Ext_Knots("EQUAL_L", deg, min(x), delta) //
               int_knots //
               Ext_Knots("EQUAL_R", deg, max(x), delta);
   end;
   else
      knots = Ext_Knots("BOUNDARY", deg, min(x)) //
              int_knots //
              Ext_Knots("BOUNDARY", deg, max(x));
   return( knots );
finish;
 
/* Helper: Function to create internal knots for B-splines.
   Valid methods are "PCTL", "EQUAL", or "RANGE".
*/
start Knots_Internal(x, method, _parm);
   parm = colvec(_parm);
   if upcase(method)="PCTL" then do;
      if nrow(parm)=1 then 
         call qntl(knots, x, T(1:parm)/(parm+1));   /* PERCENTILES(n) */
      else 
         call qntl(knots, x, parm/100);             /* PERCENTILELIST(v) */
   end;
   else if upcase(method)="EQUAL" then
      knots = min(x) + range(x)*T(1:parm)/(parm+1); /* EQUAL(n) */
   else if upcase(method)="RANGE" then 
      knots = min(x) + range(x)*parm;               /* RANGEFRACTIONS(v) */
   else
      knots = .;                                    /* unsupported method */
   return( knots );
finish;
 
/* Helper: Function to create external or boundary knots for B-splines.
   Valid methods are "BOUNDARY", "EQUAL_L", "EQUAL_R" */
start Ext_Knots(method, deg, value, delta=);
   if upcase(method)="BOUNDARY" then
      ext_knots = repeat(value, deg);
   if upcase(method)="EQUAL_L" then
      ext_knots = value - T((deg-1):0)*delta;
   if upcase(method)="EQUAL_R" then
      ext_knots = value + T(0:(deg-1))*delta;
   return( ext_knots );
finish;
 
/* read in example data */
use &DSName;
   read all var "&XVar" into x;
close;
deg = 2;
 
/* get knots for common placement schemes */
k_Pctl3     = Knots(x, deg, "PCTL", 3);
k_PctlList  = Knots(x, deg, "PCTL", {20, 50, 80});
k_eq3       = Knots(x, deg, "EQUAL", 3);
k_RangeList = Knots(x, deg, "RANGE", {0.25, 0.50, 0.75});
 
/* display the results in a single table */
results = k_Pctl3 || k_PctlList || k_eq3 || k_RangeList;
rowHeader = {'L Ext 1','L Ext 2','Int 1','Int 2','Int 3','R Ext 1','R Ext 2'};
colHeader = {'Pctl(3)','Pctl List','Equal(3)','Range List'};
print results[c=colHeader r=rowHeader L="Knot-Placement Methods"];

I reviewed the previous article and verified that these are the same knot positions that are generated by the EFFECT statement.

A design matrix for B-splines in SAS IML

If you plan to use B-splines in a regression model, you need to create the design matrix for the X variable. To generate a design matrix, you first specify a knot-placement scheme, which gives you a vector of knot locations. You can then call the BSPLINE function to produce the spline basis and evaluate the basis function at each observed value of X. The following statements show an example:

/* to generate a design matrix, first get the knots, then evaluate the splines on the X data */
knots = Knots(x, deg, "PCTL", {20, 50, 80});
basis = bspline(x, deg, knots);    /* 428x6 matrix. Columns are spline functions */

The BASIS matrix is the design matrix for the splines. It does not include an intercept column. For a multivariable regression model, you would horizontally concatenate this design matrix with the data vectors or design matrices for additional regressors.

Visualize the B-spline basis

If you want to visualize the spline basis, you can add X to the design matrix, sort the columns, and then repeat the wide-to-long visualization from earlier in this article:

/* if you want to visualize, sort the data by the X variable, then convert wide to long */
M = x || basis;   /* concatenate design matrix with X, then sort by X */
call sort(M, 1);
t = M[,1];
basis = M[,2:ncol(M)];
run WideToLong(tt, Value, Spline_Number,              /* output (long) */
               basis, t, "Spl1":"Spl6");      /* input (Basis is wide) */
title "B-Splines (deg=2)";
title2 "Knots Placed by using PERCENTILELIST({20, 50, 80})";
call series(tt, Value) group=Spline_Number grid={X Y} 
                       label={"&XVar" "Spline Value"};

Summary

This article shows how to write a function in SAS IML that reproduces common knot-placement methods for B-splines. The B-splines use both Internal and external (or boundary) knots. Internal knots are placed within an interval, which is often the observed range of the data variable. Boundary knots are stacked at the endpoints of the interval. True external knots are placed outside the interval. After you generate the knots, you can then call the built-in BSPLINE function in IML to generate a spline design matrix.

Since most SAS regression procedures support the EFFECT statement, why would you want to generate knots and splines in IML? I can think of a few reasons:

• Varying the knot-placement scheme is useful for understanding whether the parameter estimates in a spline regression model are sensitive to the knot-placement method.
• If you want to generate your own knot-placement scheme, you can modify the IML functions in this article.
• By moving the generation of the spline design matrix into IML, you can easily run simulation studies and bootstrap analyses.
• B-splines are not the only type of spline. If you want to build a regression model that uses a new type of spline, you can reuse the standard knot-placement schemes for the new spline type.

The post Specify knots for a spline basis appeared first on The DO Loop.

Splines are useful tools for fitting regression models to data. A spline replaces a single variable (call it X) with several other variables, which are a spline basis for X. When using a spline basis, the shape and location of the basis functions depend on the placement of knots. Knots are breakpoints that partition the range of a variable into subintervals. The regression model can fit the behavior of the data separately on each subinterval.

This article explains, compares, and visualizes common knot-placement schemes for regression splines.

Previous articles about regression splines have visualized natural cubic splines and truncated power functions (TPF). This article focuses on B-splines. B-splines not only require placing "internal" knots which are within the range of the data for X, but also use "external" or "boundary" knots. Consequently, this article also discusses how the placement of external knots affects the spline basis. A typical visualization is shown to the right for a B-spline basis. The vertical reference lines indicate the locations of knots.

For previous articles about splines, see

• "Regression with restricted cubic splines in SAS": How to use the EFFECT statement in SAS to define and use splines in regression models. This article discusses the KNOTMETHOD= option, which enables you to specify the placement of the knots.
• "Visualize a regression with splines": How to visualize the spline basis and interpret the parameters in a regression model that includes splines.

Why use splines for regression?

Why are splines useful in regression modelling? Briefly, splines enable you to use linear regression models to fit data that have local or even nonlinear effects. In an elementary linear model, you model the response variable (Y) as a linear combination of one or more independent variables. In the simplest case of one independent variable, the model is Y = b0 + b1*X1 + error. This model is only reasonable if Y and X are linearly related across the entire range of X, which is not always true. Splines enable you to fit a piecewise-linear model where the model changes on various subintervals of the X range. The knots define the subintervals. You can also use knots and splines to fit a piecewise-polynomial model. This is shown in subsequent sections. It is possible because the basis functions for B-splines (and other polynomial-based splines) are formed by piecewise polynomials of degree d that are joined together so that they are continuously differentiable up to order d-1.

Sample data

I'll demonstrate the knot placement and the spline visualization by using the Horsepower variable in the SasHelp.Cars data set, which has 428 observations. In this article, the phrase "the X variable" refers to the Horsepower variable. So that you can easily run the examples on your own data, I have put the name of the data set and the name of the variable into macro variables (DSName and XVar, respectively). By modifying those two macro variables, the code in this article should work for any data set and any numerical variable.

The following SAS statements define the data and the macro variables:

/* example data */
data Have;
set sashelp.cars;
keep Horsepower;
run;
 
/* define these macro variables to point to your data set and numerical variable */
%let DSName = Have;
%let XVar = Horsepower;
 
/* you should not have to modify any code below this line */
data SplDS / view=SplDS;
   set &DSName;
   _Fake = 0;    /* add fake variable for the MODEL statement in PROC GLIMMIX */
run;

The code creates a DATA step view that adds a dummy variable (_Fake) to the data set. You can use PROC MEANS to obtain the five-number summary of the data:

proc means data=SplDS min Q1 median Q3 max ndec=0;
   var &XVar;
run;

The range of the data is [73, 500].

Overview of knot placement methods

In SAS regression procedures, the location of knots is determined by the KNOTMETHOD= option in the EFFECT statement. In this article, the examples always use the BASIS=BSPLINE and DEGREE=2 options. In practice, many people use DEGREE=3, which results in more basis functions.

The knot placement options are as follows:

• PERCENTILES(n) and PERCENTILELIST(percentile-list): The knot locations are placed at percentiles of the X variable. For a B-spline basis, the external knots are placed at min(X) and max(X).
• EQUAL(n): The knot locations are placed at n evenly spaced points inside the interval [min(X), max(X)]. For a B-spline basis, the external knots are placed at evenly spaced points less than or equal to min(X) and greater or equal to max(X). The distance between all knots (exterior and interior) is the same.
• RANGEFRACTIONS(fraction-list): If the fraction list is (f1, f2, f3,...), the knot locations are placed at min(X) + (max(X)-min(X))*f_i, where 0 < f_i ≤ 1. For a B-spline basis, the external knots are placed at min(X) and max(X).
• LIST(number-list) and LISTWITHBOUNDARY(number-list): The knot locations are placed at specified locations in data coordinates. For a B-spline basis, the LIST option places external knots at min(X) and max(X). The LISTWITHBOUNDARY option enables you to specify the location of the internal and external knots.

The PERCENTILES and PERCENTILELIST methods

The knot locations are placed at percentiles of the X variable. For a B-spline basis, the external knots are placed at min(X) and max(X). The following call to PROC GLIMMIX includes an EFFECT statement. The KNOTMETHOD= option specifies the PERCENTILES(3) option, which places knots at the 25th, 50th, and 75th percentile of the X variable. The OUTDESIGN= option in the PROC GLIMMIX statement creates a SAS data set that contains the spline bases evaluated at the values of the X variable. For a DEGREE=2 basis, the columns of the data set are named Spl1--Spl6. The NOINT option in the MODEL statement prevents the OUTDESIGN= data set from containing a constant column of 1s.

/* PERCENTILES(n) specifies the locations for internal knots as n even percentiles of the data. 
   For a B-spline basis, the extremes of the data are used as boundary knots.
*/
proc glimmix data=SplDS NOFIT
             outdesign(names X=Spl)=Pctl_BSpl;
   effect Spl=spline( &XVar / basis=BSpline details degree=2
                      knotmethod=percentiles(3) ); /* <== KNOTMETHOD */
   model _Fake = Spl / NOINT;
   ods select SplineKnots BSplineDetails;
   ods output SplineKnots=KnotPositions;
run;

The ODS SELECT statement selects two output tables, which are shown above. The SplineKnots table provides information about the placement of knots. For the PERCENTILES(3) method, there are d knots placed at the minimum value of the X variable, 3 interior knots placed at the 25th, 50th, and 75th percentiles of the X variable, and d knots placed at the maximum value of the X variable. The exterior or boundary knots are indicated by an asterisk. Notice that it is mathematically acceptable to stack two or more knots at the same position.

The BSplineDetails was created by using the DETAILS option in the EFFECT statement. It tells you that the basis consists of 6 spline functions. The first basis function has support (that is, nonzero values) on the interval defined by the first three knots, which is [min(X), 25_th_Pctl]. The second basis function has support on the interval defined by the knots 1-4, which is [min(X), 50_th_Pctl]. This information continues until the sixth basis function, which has support on the interval defined by the knots 5-7, which is [75_th_Pctl, max(X)].

You can visualize the spline basis on the range of the X variable. If X has widely spaced values, the splines might look "chunky" when evaluated at the values of X. However, the basis functions are piecewise polynomials in the intervals between the knots, so they are continuous and smooth. For a spline basis of degree d, the basis functions are continuously differentiable d-1 times at the knot locations. The following statements create the visualization by using the %SplineViz macro, which is defined in the Appendix.

title "Three Interior Knots: percentile(3)";
title2 "Exterior Knots on Boundary";
%SplineViz( Pctl_Knots );

The six spline basis functions are shown. The vertical lines indicate the placement of knots. The knots define subintervals on which the splines are defined. The splines look smooth on the left because X has many unique values there. The splines look "chunky" on the right because there are only a few unique values of X that are greater than 350.

The shape of the interior splines might remind you of probability distributions. Some are roughly "bell-shaped" whereas others are skewed. Because the exterior knots are clamped at the extremes of the data, the first and last spline functions look like truncated polynomials that have the values 0 and 1 at the endpoint of the subinterval on which they are defined.

The EQUAL method

The KNOTMETHOD=EQUAL(n) option is similar to the PERCENTILES(n) option, except the knot locations for EQUAL(n) are placed at n evenly spaced points inside the interval [min(X), max(X)]. For a B-spline basis, the external knots are also evenly spaced points. The leftmost knots are less than or equal to min(X); the rightmost knots are greater than or equal to max(X). The distance between all knots (exterior and interior) is the same.

For example, the following call to PROC GLIMMIX uses the EQUAL(3) option. Otherwise, the syntax to PROC GLIMMIX is unchanged.

/* EQUAL(n) specifies that n equally spaced knots be positioned between the 
   extremes of the data. For a B-spline basis, any needed boundary knots 
  continue to be equally spaced. You can use the DATABOUNDARY option to override this behavior.
*/
proc glimmix data=SplDS NOFIT
             outdesign(names X=Spl)=Equal_BSpl;
   effect Spl=spline( &XVar / details basis=BSpline degree=2
                              knotmethod=EQUAL(3) ); /* <== KNOTMETHOD */
   model _Fake = Spl / NOINT; 
   ods select SplineKnots;
   ods output SplineKnots=KnotPositions;
run;

The output shows that the knots are equally spaced. The distance between adjacent knots is Δx = 106.75. The exterior knots on the left are at min(X)-Δx and min(X). The interior knots are at min(X) + i*Δx for i=1,2,3. The exterior knots on the right are placed at max(X) and max(X)+Δx. The following statements visualize the spline basis for this set of knots:

title "Three Interior Knots: equal(3)";
title2 "Exterior Knots Are Equally Spaced";
%SplineViz( Equal_BSpl );

These basis functions are essentially identical, ignoring the "chunky" artifacts on the right. The functions differ by a shift.

The RANGEFRACTIONS method

The RANGEFRACTIONS method is similar to EQUAL method but enables you to set knots at unequal positions within the range of the X variable. If the fraction list is (f1, f2, f3,...), the knot locations are placed at min(X) + (max(X)-min(X))*f_i, where 0 < f_i ≤ 1. So, you can recreate the EQUAL(3) spacing by using RANGEFRACTIONS(0.25 0.5 0.75). If you want unequal spacing, you can use a syntax such as RANGEFRACTIONS(0.2 0.5 0.8) to position the knots at 20%, 50%, and 80% of the range.

For a B-spline basis, the RANGEFRACTIONS method places the external knots at min(X) and max(X). Let's place three interior knots at evenly spaced positions within the range but clamp the exterior knots at the extremes of the data. This will enable us to compare spline bases that differ only by the exterior knots:

/* RANGEFRACTIONS(fraction-list) specifies that internal knots be placed at each fraction of the range of X.
   For a B-spline basis, the data extremes are used as boundary knots.
*/
proc glimmix data=SplDS nofit
             outdesign(names X=Spl)=Range_BSpl;
   effect Spl =spline(&XVar / basis=BSpline details degree=2
                                  knotmethod=rangefractions(.25 .5 .75) );
   model _Fake = Spl / NOINT; 
   ods select SplineKnots;
   ods output SplineKnots=KnotPositions;
run;

The output confirms that the interior knots are placed the same as for the EQUAL(3) option, but the exterior knots are clamped at the extreme values of the data. The following statements visualize the spline basis:

title "Three Interior Knots: rangefractions(.25 .5 .75)";
title2 "Exterior Knots Are Clamped at Min/Max";
%SplineViz( Range_BSpl );

This spline basis differs from the EQUAL(3) basis, which was highly symmetric. Placing the exterior knots at the boundary of the data means that the range of the first and last spline is [0,1]. The interior splines are not the same size and shape. (Or, perhaps I should quip that they are KNOT the same shape!)

The LIST and LISTWITHBOUNDARY methods

The LIST and LISTWITHBOUNDARY methods enable you to place the knots by using the data scale. I will not show an example. You merely list the data values at which you want the knots.

Summary

This article shows how to position knots for regression splines by using the KNOTMETHOD= option in the EFFECT statement, which is supported by many SAS regression procedures. You can use the OUTDESIGN= option in the PROC GLIMMIX statement to output the splines (evaluated at the data locations) to a data set. The output from PROC GLIMMIX includes the location of the knots. This article showed several common ways to place knots and, for each method, visualized the resulting set of spline basis functions.

Appendix: The SplineViz macro

The following macro visualizes the spline bases that are generated by each call to PROC GLIMMIX. The name of the design matrix is passed in; inside the macro it is referred to as &DS. The columns for the spline basis are named Spl1, Spl2, Spl3, .... The positions of the knots are read from the KnotPositions data set, which was created from the ODS OUTPUT statement, which writes the SplineKnots table to a data set. The macro creates a series plot that shows the spline basis functions evaluated on the values of the data variable, X.

/* macro to visualize the spline bases that are generated by each call to PROC GLIMMIX. 
   The original data should be sorted by the X variable.
   The name of the design matrix is passed in as &DS.
   The columns for the spline basis are named Spl1, Spl2, Spl3, ...
*/
%macro SplineViz( DS );
/* 1. Transpose design matrix from wide to long. The new variables are 
      SplNum : identifies the i_th basis function, i=1, 2, 3, ...
      Value  : the values of the X variable at which the i_th spline is evaluated
      Y      : the value of the i_th spline function evaluated at X
*/
data BSplineLong;
   set &DS;
   array basis Spl:;
   value = &XVar;
   do SplNum = 1 to dim(basis);
      Y = basis[SplNum];
      output;
   end;
   keep SplNum Value Y;
   label SplNum="Spline Number" Value="&XVar" Y="Spline Value";
run;
 
proc sort data=BSplineLong;
   by SplNum Value;
run;
/* add a label variable. You will need to change the format of these values
   if the scale of your data is larger or smaller than the example data.
*/
data SplineViz;
   set BSplineLong KnotPositions(keep=&XVar rename=(&XVar=KnotPos) );
   length KnotPosLabel $3;
   KnotPosLabel = putn(KnotPos, "3.0");
run;
 
proc sgplot data=SplineViz;
   series x=Value y=Y / group=SplNum nomissinggroup;
   fringe Value;
   refline KnotPos / axis=x labelloc=inside label=KnotPosLabel labelattrs=(size=8);
   yaxis max=1;     * for Viz, put all plots on same vertical scale;
run;
%mend;

The post Visualize the placement of knots for regression splines appeared first on The DO Loop.

Is bias a big deal?

For data analysis tasks, I don't often worry about biased statistics. Bias is defined as the difference between the expected value of an estimator and the population parameter. For most statistics, the standard error of the statistic is much larger than its bias, thus bias isn't very important in practice. Furthermore, many statistics are asymptotically unbiased, so the bias is small for large samples.

However, knowing whether a statistic is biased is important for Monte Carlo simulation studies. In a simulation study, you know the exact value of a population parameter. When I develop a simulation, I often verify its correctness by running many simulations and checking that the Monte Carlo estimate converges to the population parameter. However, a biased statistic does not converge to the population parameter! The Monte Carlo estimate converges to the expected value of the statistic. For a biased estimator, the expected value is different from the population parameter, by definition. As the number of Monte Carlo iterations (B) increases, the standard error of the simulation shrinks toward zero, but the bias remains constant.

The median of an exponential distribution

I recently wrote an article about the convergence of Monte Carlo estimates. The example in that article studied the sample mean for random samples drawn from an exponential distribution that has a scale parameter (λ) with the value 10. You can look up in Wikipedia or a textbook that the median of the exponential distribution is log(2)*λ. So, when λ=10, the median of the Expo(10) distribution is 6.9315.

In the previous article, the sample size was N=100. This article uses N=99 instead. Since the sample median is the "middle value" in the sample, using an odd sample size makes the theoretical formulas a little easier to understand.

A Monte Carlo simulation for the median

To demonstrate that the sample median is a biased statistic, let's run the following Monte Carlo simulation:

• Choose B=10,000 as the number of simulations.
• Draw B random samples of size N=99 from an Expo(λ=10) distribution.
• For each random sample, compute the sample median.
• Average the sample statistics to obtain the Monte Carlo estimate of the population median.
• Compare the Monte Carlo estimate to the known value of the distribution's median. The difference is the bias of the sample median in samples of size 99.

/* Monte Carlo simulation for median of EXPO(10) distribution. Estimate bias in the 
   sample median for samples of size N=99.
   The true median is log(2)*lambda ~ 6.93147. */
proc iml;
lambda = 10;                   /* generate random sample of size N from Expo(lambda) */
expo_median = log(2)*lambda;
 
call randseed(12345);
N = 99;                        /* sample size is odd ==> median is middle order statistic */
B = 10000;                     /* number of Monte Carlo samples */
 
sample_median = j(B,1,.);
x = j(N,1);
do i = 1 to B;
   call randgen(x, "Expo", lambda);
   sample_median[i] = median(x);
end;
 
MC_est = mean(sample_median);
MC_SE_est = std(sample_median) / sqrt(B); /* std error for Monte Carlo estimate */
width = 3*MC_SE_est;                      /* 3*SE ==> 99.7% CI */
Lower_CL = MC_est - width;
Upper_CL = MC_est + width;
print MC_est expo_median (MC_est-expo_median)[L="Est Bias"];
print MC_est MC_SE_est width Lower_CL Upper_CL;

The simulation shows two facts. First, the Monte Carlo estimate is larger than the median of the distribution. The estimated bias is approximately 0.0457. Is that bias large? Well, "large" or "small" must be compared to some standard unit, and for a statistic the "unit" is the standard error of the statistic. The asymptotic standard error of the median statistic is approximately λ/sqrt(N) ~ 1.005. So, the bias is less than 5% of the standard error. In other words, if you are analyzing one single data set of size N=99, the standard error of the median is 20 times larger than the bias. This is why we typically can ignore bias.

The second fact is the accuracy of the Monte Carlo estimate. After 10,000 simulations, a 99.7% confidence interval for the estimate has a half-width of 0.03. We are fairly confident that the expected value of the sample median is within 0.03 of the computed Monte Carlo estimate. The next section shows that the expected value for the sample median is 6.982, which is indeed inside the confidence interval for this set of random trials.

Order statistics for the exponential distribution

For the exponential distribution, you can use order statistics to explicitly compute the expected value of the sample median. To simplify matters, let's assume that N=2*k - 1 is an odd number. To simplify notation, let's sort the data and denote the values by
X₁ ≤ X₂ ≤ ...≤ X_k ≤ ...≤ X_N.

The median is the "k_th order statistic," which is simply the k_th element in the sorted list of values. As shown in a Wikipedia article on order statistics, the expected value of the k-th smallest observation is known to be:
E[X_k] = λ ∑_i=1^k 1 / (N – i + 1)

For N=99, k=50. The following SAS IML function evaluates the expected value of the sample median for an Expo(10) sample:

proc iml;
/* For a sample of size N drawn from an Expo distribution,
   the expected value of the k_th smallest obs (the k_th order statistic)
   is given by an explicit formula. 
*/
start EExpOrderStat(k, N, lambda=1);
   i = T(1:k);
   return( lambda * sum( 1 / (N - i + 1) ) );
finish;
 
N = 99;
lambda = 10;
/* expected value of the median for an ODD sample ~ Expo(lambda) */
k = (N+1)/2;
Emedian_N99 = EExpOrderStat(k, N, lambda);
print Emedian_N99;

The expected value is 6.98172. Notice that this value is inside the 99.7% confidence interval from our Monte Carlo simulation.

Convergence of a biased statistic

What does the convergence look like if you graph it? You will see an estimate that converges, but not to the parameter value. Instead, it converges to the expected value of the statistic. The following graph shows the Monte Carlo estimate of the sample median after k simulations, k=1..B. You can see that the estimate converges, but there is a big gap (the bias) between the estimate and the population parameter.

Conclusion

For symmetric distributions, the sample median is an unbiased estimator of the population median. Furthermore, the estimator is asymptotically unbiased, which means the bias is negligible for very large samples. However, for right-skewed distributions like the exponential distribution, the sample median is biased upwards. This bias is important for Monte Carlo simulations because the Monte Carlo estimate converges to the expected value of a statistic. When the statistic is biased, the Monte Carlo estimate converges to a value that is different from the population parameter. This is most apparent when you are simulating small samples.

Appendix: The standard error of the median for an exponential distribution

In the text, I claimed that the asymptotic standard error of the sample median of Expo(λ) data is approximately λ/sqrt(N). But wait, isn't that the well-known formula for the standard error of the mean? We want the standard error of the median!

The formula is correct. The formula relies on the fact that the sampling distribution of the median is asymptotically normal for continuous distributions. This asymptotic normality is true for all percentiles. It is the basis for the CIPCTLNORMAL option in PROC UNIVARIATE, which provides confidence intervals, assuming the sampling distribution of the percentiles are normal. If you like math, you can read the proof in these lecture notes from McGill University.

Asymptotically, the distribution of the sample median is normally distributed with standard error
SE_med = 1 / (2*sqrt(N)*f(M))
where f is the PDF of the continuous distribution and M is the true population median. (Obvious, we need f(M) ≠ 0.) For the exponential distribution, f(x) = (1/λ) exp(-x/λ) and M = λ*log(2). Thus, a several terms cancel and
f(M) = 1/(2λ)

Consequently, the asymptotic standard error is SE_med = λ/sqrt(N), which is the same formula as the standard error of the mean.

The post The sample median is a biased estimator for skewed distributions appeared first on The DO Loop.

The Goldilocks Principle for Monte Carlo methods

The Goldilocks Principle is a name given to efforts to achieve a balance between extremes. For example, suppose my goal is to ensure that my house is clean. There are two extremes. If I spend every hour of every weekend cleaning, dusting, and vacuuming, my house will look pristine. However, I expend a lot of time and effort. In contrast, if I spend no time cleaning, the house will be dirty. Somewhere between those two extremes is the Goldilocks strategy: I exert enough effort to make the house look nice, but the effort doesn't consume all my resources. The optimal amount of effort depends on how clean I want the house to be.

The same trade-offs apply in many areas of computational statistics. In Monte Carlo simulations, you can get better estimates if you run more simulations, but more simulations require more time and computing resources. A traditional Monte Carlo simulation looks like this:

1. Choose the number of simulations, B, that you plan to run. This is essentially a GUESS!
2. Generate a random sample of size N and compute statistics for that sample.
3. Repeat this B times. Store all sample statistics in an array.
4. Compute the Monte Carlo mean and report a measure of accuracy, such as the standard error or a confidence interval.

The problem is we do not know a priori how large to choose B. If B is too small, the simulation study is underpowered, and the final Monte Carlo estimate is not very accurate. Consequently, you might need to repeat the simulation with a larger value of B to obtain a more accurate estimate. On the other hand, if B is too large, the Monte Carlo estimate is very accurate, but you wasted time and CPU costs needlessly.

By using Welford's update formula inside a simulation loop, you don't have to guess a value for B, number of Monte Carlo simulations. Instead, you treat the statistic from each random sample as a piece of streaming data. You keep track of the mean and variance of the Monte Carlo estimate. When a confidence interval for the estimate is sufficiently small, you stop the simulation loop. This is sometimes called an adaptive Monte Carlo simulation.

Example of the traditional Monte Carlo simulation in SAS

Before looking at the adaptive Monte Carlo method, let's run a traditional simulation study. Suppose you have data (N=100) that you believe to be exponentially distributed. You compute the mean of the sample to estimate the population mean. A classic problem in statistics is to infer how close the estimate is to the true population mean. To do that, you need to understand the sampling distribution of the mean in a random exponential sample of size N=100. There are several ways to do that, but I want to use a Monte Carlo simulation.

The following Monte Carlo simulation uses PROC IML to generate B=10,000 random samples from the Expo(λ=10) distribution, where λ is a scale parameter. The program computes the sample mean for each random sample. It then computes the Monte Carlo estimate of the statistic by taking the mean of the 10,000 samples statistics. Lastly, it estimates the standard error of the Monte Carlo estimate and uses that to form a 99.7% confidence interval for the estimate.

/* traditional Monte Carlo simulation: Compute the Monte Carlo estimate of the 
   mean of a random sample of size N. Do this a large number of times (B) to ensure
   that the estimate is accurate. 
 
   For this example, we sample from the Exp(lambda) distribution where lambda=10 
   is a scale parameter. The true mean is lambda. 
*/
proc iml;
call randseed(1234);
N = 100;                       /* sample size */
B = 10000;                     /* number of Monte Carlo samples */
lambda = 10;                   /* scale parameter for Expo distrib */
 
sample_mean = j(B,1,.);
x = j(N,1);
do i = 1 to B;
   call randgen(x, "Expo", lambda);
   sample_mean[i] = mean(x);
end;
 
MC_est = mean(sample_mean);           /* MC estimate of statistic */
SE_est = std(sample_mean) / sqrt(B);  /* MC standard error */
width = 3*SE_est;                     /* half-width of CI: 3*SE ==> 99.7% CI */
Lower_CL = MC_est - width;
Upper_CL = MC_est + width;
print MC_est SE_est width Lower_CL Upper_CL;
 
/* we saved all the Monte Carlo estimates, so plot the sampling distribution */
title "MC Distribution";
title2 "N=100;  B=10,000";
call histogram(sample_mean) other="refline 9.9860983 9.9560683	10.016128/axis=x;";

What do we learn from this output about the accuracy of the Monte Carlo estimate? First, the histogram shows that the sample mean in random samples of size N=100 is quite variable, but approximately normally distributed. The Monte Carlo estimate of the population mean is the average of all these sample values. The Monte Carlo estimate and a 99.7% confidence interval are overlaid on the histogram. After B=10,000 simulations, a 99.7% confidence interval for the Monte Carlo estimate has a half-width of about 0.03.

If we want the Monte Carlo estimate to be more accurate, we can run additional simulations. How many simulations would we need for the 99.7% CI to have a half-width of 0.01 units? Theory tells us that the Monte Carlo standard error is proportional to 1/sqrt(B). The standard interpretation of this result is that if you want to cut the standard error in half, you must quadruple the number of simulations. In this case, we want to reduce the standard error from 0.03 to 0.01 (a ratio of 1/3), so we expect this to happen if we perform 9 times as many simulations, or roughly B=90,000.

These calculations are based on having run an initial simulation study with low power. This is often a good idea. But if you monitor the standard error during the simulation loop, you might not need the initial study. Instead, you can keep iterating until a confidence interval for the standard error is below some pre-determined tolerance.

Adaptive Monte Carlo simulation in SAS

This section implements an adaptive Monte Carlo simulation in SAS IML. The main difference between an adaptive program and a standard Monte Carlo simulation is that the adaptive program uses a DO-UNTIL statement to repeat the trials until a confidence interval for the estimate is sufficiently small.

In the previous example, the accuracy of the MC estimate after 10,000 iterations is about 0.03, as estimated by the half-width of a 99.7% confidence interval. Instead of running the entire simulation and then checking the accuracy at the end, you can use Welford's formula to provide a running standard error. You can use DO-UNTIL logic to loop until the Monte Carlo estimate is sufficiently accurate

The following IML program implements this adaptive method. Inside the iteration loop, Welford's formula updates the mean and variance. The program always generates at least 1,000 random samples (minB) but never computes more than 100,000 samples (maxB).

/* Adaptive Monte Carlo simulation using Welford's formula */
proc iml;
call randseed(1234);
N = 100;                       /* sample size */
maxB = 100000;                 /* maximum number of Monte Carlo samples */
minB =   1000;                 /* run at least this many trials */
lambda = 10;                   /* scale parameter for Expo distrib */
 
x = j(N,1);
m = 0; varsum = 0;             /* Welford's running totals */
tol = 0.01;                    /* desired accuracy */
err = 2*tol;                   /* make sure err > tol to enter loop */
do k = 1 to maxB until(err < tol);
   /* you can put ANY simulation here */
   call randgen(x, "Expo", lambda);
   y = mean(x);
   /* use Welford's online mean and variance formulas to monitor convergence */
   varsum = varsum + (k-1) * (y - m)**2 / k;  /* update the sum of squared differences */
   m = m + (y - m) / k;                       /* update the mean */
   if k > minB then do;
      sd = sqrt(varsum/(k-1));
      err = 3*sd / sqrt(k);                 /* 3*SE ==> 99.7% confidence */
    end;
end;
 
print k m err[L='error'];
quit;

By using Welford's streaming mean and variance formulas inside the DO-UNTIL loop, the simulation runs for exactly as long as required. When the tolerance is 0.01, the simulation (for this random-number seed) generates 90,374 samples. This is an amazingly close agreement with the theoretical estimate, which predicted 90,000 samples.

Summary

This article shows a modern application of Welford's formulas for the streaming mean and variance of data. The application is to iterate an algorithm until the accuracy of an estimate is less than some specified tolerance. Here, the algorithm is a Monte Carlo simulation, and the estimate is the Monte Carlo estimate of a sample mean.

I recently used this approach while computing the Monte Carlo approximation of an integral. I wanted to ensure that my integral estimate was accurate to within a certain tolerance. By monitoring the estimate and its standard error, I was able to guarantee the accuracy of the integral calculation.

Appendix: The width of the confidence interval

I didn't want this article to be too long, but I wanted to show a graph of the running Monte Carlo estimate and its 99.7% confidence interval. The Monte Carlo estimate converges quickly, as shown in the graph below. However, the estimate bounces around. You can see that when B=5000, the MC estimate is very close to the population mean, which is 10. However, the confidence interval for the estimate (shown by the gray-blue band) has a half-width of about 0.025. Indeed, as we add more simulations, the MC estimate jitters up and down, so we cannot say with confidence that we have an accurate estimate.

You can make a similar graph all the way out to the 90,000th simulation, but the half-width changes so slowly that it becomes hard to discern that it is decreasing at all. Instead, the following graph shows ONLY the half-width of the CI versus the number of simulations. In this graph, you can see the slow rate at which the CI decreases. The half-width doesn't drop below 0.01 until after 90,000 simulations.

The post On monitoring the convergence of Monte Carlo simulations appeared first on The DO Loop.

In statistical programming, we often assume that the data are stored in a matrix or a data set and can be read at will. For example, if you want to compute a sample mean or standard deviation, you simply pass a vector to a built-in function, and the software spits out the statistics. The function that computes the mean passes through the data once. The textbook formula for the standard deviation requires two passes through the data.

Modern computers are so fast, and RAM is so cheap and plentiful, that today's data analysts probably don't think twice about using multi-pass algorithms to compute statistics. But it didn't used to be that way!

This year SAS is turning 50 years old. In the 1970s, RAM was expensive. It was also small: A mainframe computer might have 64KB or less! Furthermore, if the data was on a physical media such as a tape, rewinding the tape to read the data again was a time-consuming operation. Consequently, many research papers in the 1960s and '70s focused on how to compute statistics with one pass through the data. One such paper was by Welford (1962, Technometrics), who showed a simple method for calculating corrected sums of squares and products. This one-pass method could be used to compute the mean, variance, and other central moments without using any arrays or extra memory.

An old proverb states, "everything old is new again," so perhaps it is not surprising that Welford's one-pass formula for the mean and standard deviation has applications in the modern era. One application is streaming data. If you want to compute statistics for a utility meter or a weather station that monitors conditions continuously, it would be inefficient to store all the data. It is better to compute statistics on-the-fly, if you can, and maybe store the hourly mean value.

Another modern application is for Monte Carlo simulations. Often, we guess how many iterations we need, run the simulation, and hope that the resulting estimate is accurate enough to use. Welford's method provides a way to compute the standard error while the simulation runs, which means that we can program logic to exit the simulation when the Monte Carlo estimate satisfies a predetermined level of accuracy.

This article discusses how to use Welford's one-pass formula to compute moments for streaming data. A subsequent article shows how to use this technique to monitor and stop a Monte Carlo simulation.

The traditional "offline" mean and variance

In the old days, algorithms were sometimes described as "online" or "offline." Since these terms are before the creation of the internet, they sound a little strange to the modern programmer. In computer science, an "offline" algorithm requires that the complete data are available before the algorithm is run. In contrast, an "online" algorithm processes its input serially, one datum at a time, and doesn't store the raw data. I prefer to call this type of algorithm a one-pass algorithm that requires little or no auxiliary storage.

The familiar "textbook" method for computing means and variances is an offline algorithm. To compute a mean, you sum the full set of N numbers then divide by N. For the variance, you pass through the data again, this time computing the sum of squared deviances (SSE) from the mean. You divide the SSE by N-1 to obtain the variance. If desired, take the square root to obtain the standard deviation. For example, the following SAS IML program shows the traditional offline algorithm:

proc iml;
y = {48, 51, 46, 37, 45, 50, 46, 51, 51, 55};
 
/* A traditional or "offline" calculation of the mean and standard deviation
   uses two passes through the data: 
   One for the mean, and a second pass for the sum of squared deviations, 
   which gives the variance or standard deviation. */
trad_m = mean(y);
trad_sse = ssq(y - trad_m);
trad_sd = std(y);
print trad_m trad_sse trad_sd;   /* sample mean = 48; StdDev ~ 5 */

This two-pass approach is perfectly fine for everyday data analysis. It requires that the data vector y is known and fits in memory. But it is possible to compute these statistics without using extra memory. The next section shows how to compute the mean and standard deviation by using a one-pass algorithm.

The Welford update formula for "online" data

Welford (1962) published a one-pass algorithm for computing the mean and variance. This method was popularized when Knuth featured it in his seminal book, The Art of Computer Programming, Volume II: Seminumerical Algorithms (1981). The Welford method maintains two running totals: one for the mean, and another for the sum of squared differences. For each new datum, the algorithm updates these running totals. The following IML program implements Welford's algorithm by explicitly looping over the elements of the y vector:

/* The Welford update formula:
   First, initialize the Welford variables for the one-pass mean and std dev */
m = 0;
varsum = 0;
 
/* Update the running totals for each new datum y[k]. */ 
do k = 1 to nrow(y);
   varsum = varsum + (k-1) * (y[k] - m)**2 / k;  /* update the sum of squared differences */
   if k > 1 then sd = sqrt(varsum/(k-1));
   m = m + (y[k] - m) / k;                       /* update the mean */
end;
print m varsum sd;

In the program, the value of m on the right side of the mean-update formula is the sample mean after seeing k-1 observations, so denote it by m_k-1. The formula overwrites m with the mean of k observations, so denote the left side as m_k. It's not hard to see why Welford's update formula is true. You can rewrite the mean of k observations by writing the sum of k observations as the k_th observation plus the sum of the earlier k-1 observations, as follows:
$\begin{aligned} m_k &= \frac{1}{k} \sum_{i=1}^{k} x_i \\ &= \frac{1}{k} \left( x_k + \sum_{i=1}^{k-1} x_i \right) \\ &= \frac{1}{k} \left( x_k + \frac{k-1}{k-1} \sum_{i=1}^{k-1} x_i \right) \\ &= \frac{1}{k} \left( x_k + (k-1) m_{k-1} \right) \\ &= m_{k-1} + ( x_k - m_{k-1} ) / k \end{aligned}$

The derivation of the formula for the variance is only slightly more complicated. You can find the derivation online.

Other useful one-pass methods

Since SAS is celebrating its 50th birthday, I should mention the importance of one-pass methods in early SAS procedures. In the 1970s, SAS was able to leverage one-pass methods to compute basic statistics, regression analyses, and multivariate analyses on data sets that had a huge number of observations. The regression and multivariate analyses were possible because SAS computes the uncorrected sum of squares and crossproducts (USSCP or X`X) matrix in one pass through the data.. From the USSCP matrix, you can use the SWEEP operator to obtain the centered SSCP, which then yields the covariance and correlation matrices. SAS used these and other one-pass methods to serve customers who needed to analyze Big Data quickly.

Summary

This article discusses the value of Welford's method for the mean and variance. Welford's method is an example of a one-pass formula. Today, memory is cheap and plentiful, so one-pass formulas are not as necessary as they were in the 1970s when SAS was founded. However, they are still useful. Modern sensors generate hundreds of kilobytes of data every minute, so it makes sense to compute statistics "on the fly" rather than storing all the raw data. In addition, the next article shows how you can use Welford's formulas in Monte Carlo simulations to monitor convergence.

The post One-pass formulas for mean and variance appeared first on The DO Loop.

Previous articles in the series provide information about the whisky data, the PCA analysis, and issues related to two-factor low-rank factorizations such as the NMF:

• PCA Analysis of whisky: What is the Scotch whisky data? How can you use a PCA to identify whiskies that have similar flavor characteristics?
• A two-factor rank-k decomposition has inherent ambiguities: You can arbitrarily scale and permute the order of the factors.
• A comparison of the PCA and NMF decompositions on simulated data, which reveals why the NMF might be preferred for nonnegative data.

In SAS, you can obtain the NMF by using PROC NMF in SAS Viya, or by implementing the algorithm directly in SAS IML. Sharad Saxena has written three blog posts about how to use PROC NMF in SAS Viya:

1. How to Understand the NMF
2. How to Extract Topics from a Set of Documents
3. How to Use the NMF to Build a Recommender System

Dimension reduction by using the NMF

A previous article uses classical principal component analysis (PCA) to analyze the flavor profiles of 86 Scotch whiskies. Each whisky is assigned a value 0-4 according to 12 flavors: Tobacco, Medicinal, Smoky, Body, Spicy, Winey, Nutty, Honey, Malty, Fruity, Sweetness, and Floral. The PCA and NMF analyses were discussed in an article by Young, Fogel, and Hawkins (2006, "Clustering scotch whiskies using non-negative matrix factorization", SPES Newsletter).

By using the PCA, you can find a low-dimensional linear subspace in the 12-dimensional variable space such that the subspace captures most of the variance in the data. This is called dimension reduction. For the whisky data, you can use PCA to choose a four-dimensional subspace that explains 67% of the variability in the data. The PCA factorization has two important properties: it maximizes the explained variance of the data within a lower-dimensional linear subspace, and it provides orthogonal basis vectors for the subspace.

Whereas the PCA can analyze any centered data matrix, the NMF factorization is designed for dimension-reduction of a nonnegative data matrix. The NMF seeks to decompose the data matrix into profile vectors (which are basis vectors for a lower-dimensional subspace) and a matrix of weights so that the product approximates the data matrix. In addition, the NMF requires that the profile vectors and weights be nonnegative. This last condition is why the NMF factors are often more interpretable than the PC vectors. However, the NMF decomposition does not have the maximum-variance or orthogonal-basis properties that make the PCA useful.

How to interpret the NMF factors

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 NMF creates a rank-4 approximation to X. The approximation is represented as a product W*H, where W is 86x4 and H is 4x12. The rows of H form a basis for a 4-D linear subspace. These are the "flavor profile" vectors that are analogous to the principal components (eigenvectors) in a PCA analysis.

The rows of the W matrix are weights (or coefficients) for the 4-factor approximation of the rows of X. If W_i is the i_th row of W, the product W_i * H approximates the flavor profile for the i_th whisky. Equivalently, if H₁-H₄ are the rows of H, then the linear combination
W_i1H₁ + W_i2H₂ + W_i3H₃ + W_i4H₄
is an approximation of the flavor profile for the i_th whisky. In this sense, W provides the weights for the approximation of observations in terms of the NMF profile vectors.

You can also interpret H as a set of weights for an approximation to the variables of X. The columns of W can be considered "latent factors" for X. The linear combination W*H[,j] is an approximation of the j_th variable, X_j.

Be aware that the role played by W and H depends on how the data are arranged in X. For example, it is equally valid to represent the data by using a 12x86 data matrix, Y, where the rows are the flavors and the columns are the whiskies. In that case, Y` = X, and Y = H`*W` is a valid NMF factorization.

An NMF analysis of the whisky data

I'm a statistician, not a Scotch connoisseur. Nevertheless, let's perform an NMF analysis of the whisky data. I used SAS IML to implement a simple version of the NMF algorithm. You can download the SAS programs in this article from GitHub and run the analysis yourself. After reading the whisky data into a matrix, the following IML statements perform an NMF analysis and visualize the H matrix. Since the rows of H are the basis vectors for a 4-D subspace, I transpose the matrix and display the columns. As shown in a previous article, you can scale the NMF factors so that they are approximately on the same scale as the PCA factors.

%let varNames = Tobacco Medicinal Smoky Body Spicy Winey Nutty Honey Malty Fruity Sweetness Floral;
proc iml;
/* Apply NMF to the Scotch whisky data */
varNames = propcase({&varNames});
use Whisky;
   read all var varNames into X;
   read all var {"Distillery" "Selected"};
close;
 
/* The PCA analysis of the whisky data kept four PCs.
   Load the nmf_mult function and compute a four-factor NMF */
load module = _ALL_;
names = 'NMF 1':'NMF 4';
run nmf_mult(W1, H1, X, ncol(names));
 
/* scale the H matrix to use the same scale used in the PCA article */
maxValPCA = 0.8506;             /* max value in the PCA eigenbasis */
scale = maxValPCA / H1[,<>];    /* scale each row by maximum, then multiply by maxVal */
W = W1 / scale`;
H = H1 # scale;
print  (W[<>,])[L="Max W Cols"],   (H[,<>])[L="Max H Rows"];
 
%let Reds = CXF7F7F7 CXFDDBC7 CXF4A582 CXD6604D CXB2182B ;
ods graphics / width=360px height=480px;
call heatmapcont(H`) title="NMF Factors (H`) for Whiskies"
                     colorramp={&Reds} range=(0||max(H))
                     xvalues=names yvalues=varNames;

For comparison, I also show a heat map of the factors (basis elements) for the four-component PCA analysis from my previous article. The colors in the left image correspond to the interval [0, 0.8506]. The colors in the right image correspond to the interval [-0.8506, 0.8506], where white represents 0, so that the red colors represent the same range in each image. The images enable you to compare the NMF and PCA factors. Here are some comparisons:

• Both images show latent factors for a four-factor approximation of the whisky flavors. That is, the flavor characteristics of all whiskies are approximated by linear combinations of these four factors.
• The factors for the NMF are all positive. Because the weights are also positive, the flavor of each whisky is approximated by using an additive sum of these NMF components. In contrast, the PCA factors have both positive and negative values. The flavor of each whisky is approximated by using a sum of these contrasts.
• The first NMF factor is a strong combination of the Fruity, Sweetness, and Floral variables, with smaller amounts of Smoky, Spicy, and Malty. At the risk of promoting stereotypes, these dominant flavors are sometimes characterized as "feminine" flavors.
• The second NMF factor is a combination of 10 flavors, with an emphasis on Body. Unlike the first factor, it includes the Winey, Nutty, and Honey flavors. This factor might be described as modeling "overall" flavors. Or perhaps as a "base" or "average" profile for Scotch whiskies. In a subsequent section, you will see that the second factor has a substantial weight for almost all whiskies.
• The third NMF factor emphasizes "masculine" flavors: Tobacco, Medicinal, Smoky. This factor also includes the Body and Spicy variables.
• The fourth NMF factor emphasizes the Nutty flavors, but also includes the Malty and Sweetness variables.
• It is difficult to compare the NMF factors with the PCA factors. The first PCA factor is a contrast between the masculine and feminine flavors. The second PCA factor seems to emphasize overall flavor. The fourth factor is primarily a contrast between Nutty and Sweetness. The fact that the PCA factors are harder to interpret is one reason why NMF is used in applications where interpretability is important.

Approximating whisky flavors in terms of factors

The goal of the PCA and NMF analyses is to choose a small number of factors that can be used to reduce the dimensionality of the data. Instead of associating each whisky with 12 variables (flavors), the PCA and NMF analyses both choose four new "compound" factors. The flavor profile of each whisky is approximated as a linear combination of these factors. A previous article discusses how to visualize the weights for the four-component PCA analysis. In a PCA analysis, the weights can be positive or negative.

For the NMF analysis, the weights are always nonnegative. The following heat map visualizes a selected subset of whiskies in terms of the NMF components.

ods graphics / width=480px height=640px;
/* choose a subset of whiskies */
selIdx = loc(Selected);
W_sel = W[selIdx, ];
Whisky_name = Distillery[selIdx];
call heatmapcont(W_sel) title="NMF Characteristics of Selected Whiskies (W)"
         colorramp={&Reds} range=(0||max(W))
         xvalues=names yvalues=Whisky_name;

For each whisky (each row), the heat map shows the weights for the approximation of the whisky's flavor in terms of the four NMF factors. You can use this heat map in place of scatter plots to visualize the four-dimensional scores for each whisky. A few features are apparent:

• The second NMF factor is used for almost all whiskies. In that sense, the second factor is a "base" set of flavors. Most individual whiskies start with this flavor profile and then add additional features by using other factors. The Aberlour, Balvenie, Glenlivet, Linkwood, and Macallan whiskies have similar flavors, which is based on the second factor. The similarity in these whiskies can be explained by noting that these distilleries are all located in the same region (Speyside) Scotland.
• The first NMF factor is dominant in the Glenfiddich whisky and five others (Auchentoshan, GlenGarioch,...). Recall that the first NMF factor emphasizes "feminine" flavors such as Fruity, Sweetness, and Floral.
• The third NMF factor is dominant in the Ardbeg, Lagavulin, and Laphroaig whiskies. Recall that the third NMF factor emphasizes "masculine" flavors. Indeed, all three whiskies are often described as "earthy" or "peaty," and the distilleries are all on the island of Islay.
• The fourth NMF factor is primarily related to the Nutty variable. This factor is not the sole component of whisky. Rather, this factor is added to other components. The Edradour whisky is primarily composed of the "base" flavor with Nutty enhancements.

The orthogonality (not!) of NMF factors

Unlike the PCA factors, the NMF factors are not orthogonal to each other. The NMF factors are correlated, as shown by the following computation:

/* show that the profile vectors are not orthogonal */
R = corr(H`);
print R[c=names r=names F=Bestd5.];

How much variance is captured by the NMF factors?

When you use a k-factor PCA for dimension-reduction, the factors form a basis for a k-dimensional linear subspace that explains the most variance in the data. Since the NMF factors are not orthogonal, the subspace spanned by the NMF factors must explain less variance. But how much less? That answer depends on the data, but you can write a short IML function that provides the answer.

The math is standard linear algebra, but a little long to explain. Let the columns of A represent a set of linearly independent basis vectors for an arbitrary subspace. To compute the variance explained by the subspace, you can project the data onto that subspace, then compute the covariance matrix of the projected data.

Let S be the covariance matrix of the original data. Briefly, you can use the following method to compute the proportion of the variance in S that is explained by the column space of A.

• Let Q be an orthonormal basis for the column space of A.
• The projection matrix onto this subspace is P = Q*Q`.
• The covariance of the projected data is P*S*P. Therefore, the explained variance is trace(P*S*P).
• The trace of a product of square matrices is invariant under a cyclic shift of the matrices in the product. If you substitute P = Q*Q`, shift the order, and simplify, you discover that the explained variance is trace(Q`*S*Q).

This method is implemented in the following SAS IML function, which accepts the covariance matrix, S, and a matrix of basis elements, A. The function returns the proportion of the variance in S that is explained by the subspace defined by the columns of A.

/* Let S be an nxn covariance matrix. 
   Let A be an nxk matrix with k linearly indep columns.
   This function returns the proportion of the variance in S that is explained by column space of A.
*/
start VarExplained(A, S);
   call svd(U, D, V, A);   /* Find orthonormal basis, Q, for the span of A */
   k = ncol(A);
   Q = U[, 1:k]; 
 
   var_explained = trace(Q`*S*Q);   /* variance explained by colspace(A) */
   total_var = trace(S);            /* total variance in S */
   prop_var = var_explained / total_var; /* proportion of var explained */
   return prop_var;
finish;
 
/* For the PCA, find the proportion of variance explained by the first 4 PCs */
call eigen(D, U, S);                 /* PCA decomp */
U = U[, 1:4];
prop_var_PCA = VarExplained(U, S);   /* somewhat redundant, since Q=U */
/* For the NMF, find the proportion of variance explained by the rowspace of H */
prop_var_NMF = VarExplained(H`, S);
print prop_var_PCA[F=PERCENT7.1] prop_var_NMF[F=PERCENT7.1];

The output shows that the maximum variance that can be captured in the 4-D linear subspace is 67.1%. By comparison, the NMF factors explain 60.3% of the variance.

Summary

The purpose of this article is to use the NMF factorization on the Scotch whisky data set, which classifies the flavors in 86 Scotch whiskies according to 12 flavor characteristics. These data were analyzed by Young, Fogel, and Hawkins (2006). The NMF decomposes the data matrix into two matrices: A matrix of factors (H) and a matrix of weights (W). The NMF enables you to reduce the dimensionality of the data by replacing the 12 original variables with four factors that are nonnegative combinations of the variables. The rows of the weights matrix provides the coefficients for a rank-4 approximation of the data. The flavor profile of each whisky is approximated by using a linear combination of factors. If you visualize the weights matrix, you can identify which whiskies have similar flavor profiles.

You can download the SAS programs in this article from GitHub.

The post An NMF analysis of Scotch whiskies appeared first on The DO Loop.

The purpose of this article is to use an example to compare the classical singular value decomposition (SVD) and the nonnegative matrix factorization (NMF). An excellent paper that presented the NMF to the statistical community is Fogel, Hawkins, Beecher, Luta, and Young (2013), "A tale of two matrix factorizations." One section of the paper that caught my eye is Section 4.2, "Realistic synthetic mixture," which compares and contrasts the NMF and the classical SVD. This article presents a similar example. It also improves the visualizations in the Fogel et al. paper. The images in this blog post are better for visualizing how the SVD and the NMF methods differ.

Simulating patient-gene data for a disease

The Fogel et al. example is based on clinical study that involves gene expression levels for 30 genes across 40 patient samples. (Actually, the paper uses 40 genes, but the example is more enlightening if the number of variables and observations are different.) Assume a disease (such as diabetes) has two distinct genetic etiologies. The L matrix ("L" for "left side") represents 40 patients. The first column of L represents a measure of one etiology. The second column of L represents a measure of a second etiology. In this example, 10 patients have only the first etiology, 10 have only the second etiology, and 10 have both.

The R matrix ("R" for "right side") represents expression levels for 30 genes. This example assumes the first 10 genes are related to the first etiology, and the next 10 genes are related to the second etiology. The last 10 genes are not related to the disease.

You can use SAS IML to simulate these factors and to construct the data matrix as their product. Because the data are simulated, we know that the data matrix is exactly rank-2, and we know a true set of rank-2 matrices whose product equals the data matrix. You can simulate the L and R matrices by randomly assigning a value between 1 and 2 to specific cells. Notice that these factors are entirely nonnegative, so the product matrix, X = L*R, is also nonnegative. The following SAS IML statements simulate the L and R matrices and create heat maps to visualize them.

proc iml;
reset fuzz;
/*** Simulation of L and R for clinical example of two etiologies ***/
call randseed(123);
nrow = 40;
ncol = 30;
L = j(nrow, 2, 0);
idx1 = (11:20) || (31:40);
L[idx1, 1] = randfun( ncol(idx1), "Uniform", 1, 2); /* 1st col */
idx2 =  21:40;
L[idx2, 2] = randfun( ncol(idx2), "Uniform", 1, 2); /* 2nd col */
 
R = j(2, ncol, 0);
idx1 = 1:10;
R[1, idx1] = randfun( 1|| ncol(idx1), "Uniform", 1, 2); /* 1st row */
idx2 = 11:20;
R[2, idx2] = randfun( 1||ncol(idx2), "Uniform", 1, 2);  /* 2nd row */
 
/* simplify L and R by rounding to 2 decimal places */
L = round(L, 0.01);
R = round(R, 0.01);
/*** end simulation of L and R ***/
 
/* create blue-white-red color ramp from ColorBrewer */
%let BlueReds = CX2166AC CX4393C3 CX92C5DE CXD1E5F0 CXF7F7F7 CXFDDBC7 CXF4A582 CXD6604D CXB2182B ;
%let Reds = CXF7F7F7 CXFDDBC7 CXF4A582 CXD6604D CXB2182B ;
 
/* Visualize L (40 rows, transposed for viewing) and R (30 columns) */
ods graphics / width=640px height=84px; 
call heatmapcont(L`) colorramp={&BlueReds} yvalues={'Real 1' 'Real 2'} range={-2.2 2.2} title="Original L`";
 
ods graphics / width=480px height=84px;
call heatmapcont(R) colorramp={&BlueReds} yvalues={'Real 1' 'Real 2'} range={-2.2 2.2} title="Original R";

The columns of the L matrix are the factors for gene activation. The rows of the L matrix represent patients. (Note: the image is for the transposed matrix, L`). You can visualize the etiologies for the four groups. For the R matrix, the rows are the factors and the columns are the genes. The heat maps are entirely red or white because all values of L and R are nonnegative. The product, X = L*R, represents the gene expression levels for 30 genes across 40 patient samples. The following statements visualize the data matrix:

X = L*R;
ods graphics / width=360px height=480px; 
/* Note: Dark red now represents 4, not 2.2 */
call heatmapcont(X) colorramp={&Reds} range={-0.01 4.04} title="Original (Synthetic) Data Matrix";

In the real world, the order of the patients and the order of the genes would be permuted arbitrarily. Also, there would be random noise that appears in the 40x30 panel of measurements. Regardless, when presented with these data, a researcher would like to recover the underlying factors that generated them. If researchers know or suspect that the disease has two etiologies, they might try to approximate the matrix by using a rank-2 decomposition of factors. In this example, we know the true factors whose product forms the X matrix, but in the real world, those factors are unknown and must be estimated from X.

The SVD factorization

The classical way to decompose X into rank-2 factors is to use the singular value decomposition (SVD). A newer decomposition is the nonnegative matrix factorization (NMF). Following Fogel et al. (2013), the remainder of this article compares these two factorizations.

A previous article discusses the fact that a low-rank factorization is not unique, but can be arbitrarily rescaled and permuted. Fogel et al. visualized the results by using heat maps. However, their heat maps are individually scaled to use the minimum and maximum value in each matrix. Because the SVD factors have positive and negative elements whereas the NMF factors do not, it is hard to compare the two factorization methods unless the colors in the heat maps represent the same scale. In this article, I carefully scale the factors so that red colors always represent positive values, white represents zero, and blue always represents negative values.

Let's look at an SVD analysis of the simulated data matrix, X. The following statements call the SVD subroutine in SAS IML to compute X = U*D*V`. To get the rank-2 approximation of X, the first two columns of U and V are extracted. For this example, the rank-2 product exactly equals X because we constructed X as the product of two rank-2 matrices and did not add any noise.

k = 2;
call SVD(A, D, B, X);
i = 1:k;   
L_SVD =   A[, i] # sqrt(D[i])`;   /* scale each factor evenly */
R_SVD = T(B[, i] # sqrt(D[i])`);
X2 = L_SVD*R_SVD;
 
/* show that the SVD recovers X exactly */
maxDiff = max(abs(X-X2));
print maxDiff;

Interpreting the SVD factors

Although the product of the SVD exactly equals X, the SVD factors do NOT look like the nonnegative factors L and R that were used to create X. The reason is that the factors in an SVD are orthogonal, which means that the factors must contain both positive and negative elements when k ≥ 2. (Technically, if X is block-diagonal, then the SVD can recover the original features, but that is not the case in this example.) Because X contains strictly nonnegative values and isn't block-diagonal, it is mathematically impossible for the SVD to return orthogonal factors that are entirely nonnegative.

This fact is shown by using the following heat maps, which display the SVD factors on the same color scale as the true factors.

/* Visualize true L` (40 columns) and the left factor from the SVD */
ods graphics / width=640px height=84px; 
call heatmapcont(L`) colorramp={&BlueReds} yvalues={'Real 1' 'Real 2'}   range={-2.2 2.2} title="Real Factor L`";
call heatmapcont(L_SVD`) colorramp={&BlueReds} yvalues={'SVD 1' 'SVD 2'} range={-2.2 2.2} title="SVD Factor";

The heat maps show (the transpose of) the original factor, L, and the SVD factor, U. Whereas the true factors for X have strictly nonnegative values, the second factor in the SVD contains negative values (the blue cells). The negative values are required because of orthogonality: the inner product of the columns must equal 0! The negative values represent contrasts, which can be difficult to interpret when the true underlying components of the data are strictly additive.

You can use similar statement to visualize the right-hand SVD factor, V, and compare it to the original factor, R. The results are shown below:

The NMF factorization

The SVD cannot return purely nonnegative factors, which makes interpretation of the factors difficult when the data are nonnegative. Researchers developed the NMF to address this issue. The NMF decomposes a nonnegative matrix into a product W*H, where the elements in W and H are constrained to be nonnegative. The decomposition is not unique.

Let's use NMF to produce rank-2 factors of X. In SAS, you can obtain the NMF by using PROC NMF in SAS Viya, or by implementing the algorithm directly in SAS IML. This article computes the NMF by running an IML module that implements the Lee and Seung multiplicative update algorithm. The factorization you get might depend on the way that you specify an "initial guess" for the factors. Consequently, the result that I present is just one possible NMF factorization.

The NMF returns two rank-2 matrices, W and H, such that X ≈ W * H. As with the SVD, the product is a rank-2 approximation to X. However, due to the algorithm and initialization that I used, the product of the NMF factors does not exactly equal X. The NMF algorithm attempts to optimize a non-convex objective function. It can converge to local minima.

/* Run the Lee and Seung "multiplicative update" algorithm for NMF.
   The IML implementation can be downloaded from GitHub. */
run nmf_mult(W1, H1, X, 2);
 
/* Let's see whether the product W1*H1 is a reasonable approximation of X. */
X_est = W1*H1;
ods graphics / width=360px height=480px; 
/* Note: Dark red now represents ~4, not 2.2 */
call heatmapcont(X_est) colorramp={&Reds} range={-0.01 4.04} title="NMF Estimate of Matrix";

According to the heat map, the product W1*H1 does a very good job approximating X. An obvious difference is seen for Patients 32 and 39. Whereas these patients actually have positive gene expression levels for genes 1-10, the rank-2 NMF approximation does not fully recover that fact.

Comparing the NMF factors

As explained in the previous article, there are two inherent ambiguities in the factors of the NMF:

1. The columns of W and the rows of H can be arbitrarily scaled, provided they are scaled reciprocally.
2. The columns of W and the rows of H can be arbitrarily permuted, provided they are permuted together.

To compare the factors from the NMF with the original factors, it is best that the colors represent the same scale. We could standardize the factors L and R, but since we have already visualized them, I will instead adjust the scaling of W and H so that their maximum elements are approximately equal to 2, which was the maximum value in the original L and R factors. This is not a standard way to standardize (pun intended!), but it is useful for this example.

/* The following statements rescale W and H to try to put all values in [0,2]. */
maxVal = 2;
scale = maxVal / W1[<>,];    /* scale each col by maximum, then multiply by maxVal */
W = W1 # scale;              /* scale columns of W */
H = H1 / scale`;             /* reciprocally scale rows of H */
 
/* Now, visualize the NMF factors, W and H, and compare with L and R from the original simulation. */
ods graphics / width=640px height=70px; 
call heatmapcont(L`) colorramp={&BlueReds} yvalues={'Real 1' 'Real 2'} range={-2.2 2.2} title="Real Factor L`";
call heatmapcont(W`) colorramp={&BlueReds} yvalues={'NMF 1' 'NMF 2'}   range={-2.2 2.2} title="NMF Factor W`";

For some data, you might need to permute the columns of the left NMF factor and the rows of the right factor before you can compare them to real factors. (An easy way to compare: sort the columns of all left factors from largest to smallest by using the Euclidean norm.) In this example, we got lucky and can immediately compare the NMF and real factors. The similarity between the left factors is impressive. The NMF factor captures the structure of the real factor except for Patients 32 and 39. In this case, the NMF algorithm converged to a solution that is not the true factor. This happens frequently. One of the drawbacks of NMF is that the solution you get depend on an initial guess for the factors, and there are many competing ideas about how to provide an initial guess.

The agreement between the right factors is also impressive, as shown below:

This output highlights two critical features of the NMF algorithm:

1. Both factors contain only nonnegative elements. It immediately follows that the product W*H must also be nonnegative.
2. The NMF decomposition does an excellent job of recovering the true underlying factors. If you compare the true factors with the NMF factors, you can see that the patterns are remarkably similar. The most obvious difference is that the first row of the transposed W factor does not correctly estimate the values for Patients 32 and 39. The NMF estimate indicates that these patients only exhibit the second etiology, whereas in reality, they exhibit both.

Summary

This article uses SAS IML to present an example that is similar to Section 4.2 in Fogel, et al. (2013). The example simulates a data matrix that represents the expression profiles of 30 target genes for 40 clinical study participants. A classical SVD analysis produces factors whose product perfectly equals X. The SVD uses deterministic linear algebra, and is guaranteed to find an optimal rank-k approximation. However, the factors contain both positive and negative values because the SVD enforces orthogonality of the factors.

In contrast, the NMF produces factors that contain strictly nonnegative values. When the true underlying features of the data involve positive quantities whose effects are additive, the NMF factors are more natural and interpretable. However, the NMF attempts to optimize a non-convex objective function. The algorithm's solution depends on how the algorithm is initialized. Some results might correspond to local minimum instead of a global one.

You can download the complete SAS program that I used to simulate the data, run the NMF algorithm, and visualize the results.

The post A comparison of SVD and NMF factorizations appeared first on The DO Loop.