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.

This article looks at a two-factor approximation, where an n x m data matrix, Y, is approximated by the product L*R, where L is an n x k matrix and R is a k x m matrix. (As a mnemonic, I use L as the "left" matrix and R as the "right" matrix.) Usually k is a small integer. This is motivated by a particular two-factor decomposition, called the nonnegative matrix factorization (NMF),which I write about in another article. Although the SVD and the NMF share similarities, they are also very different. For most data matrices, you can standardize the factors in the SVD to obtain a unique decomposition. However, two-factor NMF decomposition is more complicated. This article takes a closer look at the nonuniqueness of the factors in a two-factor decomposition.

The nonuniqueness of the NMF explains why different algorithms might give NMFs that are correct but look completely different. The tips in this article make it easier to compare factorizations that are produced by different software or different algorithms.

The Fundamental Theorem of Linear Algebra

Nonuniqueness of factors is very common in mathematics. For example, the integer 12 can be expressed in multiple ways as the product of two nontrivial factors: 2*6, and 3*4. A way to deal with nonuniqueness is to standardize the factors so that they have special properties. For example, the Fundamental Theorem of Arithmetic says that every composite integer can be represented uniquely as a product of prime numbers, up to the order of the factors. This allows us to represent 12 = 2²*3 in a standard way. Note that the representation has two properties: the factors are "standardized" to be prime, and the prime factors are listed in increasing order (2 before 3) to eliminate the order ambiguity.

In a similar way, the SVD decomposition of a matrix is sometimes called the Fundamental Theorem of Linear Algebra. The singular value decomposition states that every n x p matrix can be written as the product of three matrices: A = U Σ V^T where U and V are orthogonal matrices, and Σ is a diagonal matrix that contains the singular values. The SVD is not unique, but if the singular values of Σ are distinct, then we can obtain a factorization that is "mostly unique," in the following sense:

• Order the singular values of Σ in decreasing order. You can always do this by using a permutation that re-orders the columns of U and V.
• The columns of U and V are unit length, but they can vary in sign. Make some arbitrary choice such as "the first row of U is positive" or "the sum of the columns of U are positive." Rules like these remove the sign ambiguity.

All statistical libraries order the singular values. However, most do not standardize the signs of the columns of U and V, which means that you might see slightly different SVDs if you compare the results from different software. However, if you apply a rule to standardize the columns of U, then the results from different software will agree, provided that the singular values are distinct.

The SVD has a strong connection to low-rank approximations. Namely, if you choose some integer k ≥ 1, you can use the SVD to find the best rank-k approximation to the original matrix, Y. The new decomposition is often denoted as U_kΣ_kV`_k, where U_k is defined as the first k columns of U, Σ_k is the upper-left k x k submatrix of Σ, and V_k is the first k columns of V.

Nonuniqueness of a two-factor decomposition: The scaling ambiguity

As explained above, the SVD provides a three-factor rank-k decomposition that can be standardized so that it becomes unique. However, if you want to express the SVD by using only two factors, you must sacrifice uniqueness. The SVD represents the linear transformation as the composition of a rotation, a scaling, and another rotation. You can try to distribute the scaling transformation (Σ) into one or both of the orthogonal matrices, but there are infinitely many ways to do that. This is called the scaling ambiguity. If α is any value in [0, 1], you can form a two-factor decomposition in infinitely many ways:
Y = (U*Σ^α) (Σ^1-αV`)

You can easily construct examples to demonstrate this nonuniqueness of a two-factor decomposition, as shown in the following SAS IML program:

proc iml; 
reset fuzz;
/* Define a simple matrix with an underlying rank of 2 */
Y = { 8 11 5,
     11 12 5,
     12  9 3,
     17 14 5 };
 
/* Use SVD to find the factors */
call svd(U, Sigma, V, Y);
 
/* Extract the first k=2 columns to get a rank-2 product */
i = 1:2;
U2 = U[,i];
S  = Sigma[i];
V2 = V[,i];
 
/* You can distribute the scaling of S into the left factor, 
   the right factor, or both factors by varying alpha */
alpha = {0, 0.5, 1};
do j = 1 to nrow(alpha);
   L = U2 * diag(S##alpha[j]);
   R = T( V2 * diag(S##(1-alpha[j])) );
   /* The product L * R is identical for all choices of alpha! */
   prod = L*R;
   print prod;
end;

I do not display the output from the program, but for every assignment inside the DO loop, the product L*R is exactly the same and is equal to Y (because Y is rank-2).

A problem with two-factor decompositions

In general, suppose you decompose a matrix Y as Y = L*R, where L is n x k and R is k x m. This decomposition is not unique because you can also write
Y = (L*A) (A^-1R)
for any invertible matrix, A. Thus, let's examine ways to standardize the problem to resolve as many ambiguities as we can.

The order ambiguity

A permutation matrix is an invertible matrix that changes the orders of columns or rows. The inverse of a permutation matrix equals its transpose. Thus, for any two-factor decomposition Y = L*R, you can permute the columns of L and simultaneously permute the rows of R by using a permutation matrix, P, as follows:
Y = (L*P) (P`*R).
This is known as the order ambiguity or the permutation ambiguity. For example, the following IML program interchanges the order of the columns of L and the rows of R, and the product is the same:

/* Permute the columns of L and rows of R */
P = {0 1,
     1 0};
L_new = L * P;
R_new = P` * R;
Y_check = L_new * R_new; /* perfectly reconstructs Y */
print Y_check;

The mixing ambiguity

Recall that the SVD decomposition is (practically) unique is because the left and right factors are orthonormal. If you change the scaling of the U or V matrices, or you change the order of their columns and rows, they are still orthogonal. Consequently, if we adopt certain conventions (for example, sort the columns by the magnitude of the singular values and require that V is orthonormal), then we can standardize the SVD, provided that the singular values are distinct.

If we do not insist that the left and right factors be orthogonal, then we might have to sacrifice uniqueness of the factorization. As mentioned above, for any invertible matrix, A, it is true that Y = (L*A) (A^-1R) is a valid two-factor decomposition. The following IML program shows an example:

/* A can be ANY invertible matrix. */
A = { 0.9  0.1,
      0.1  0.9};
L_new = L * A;
R_new = inv(A) * R;
Y_check = L_new * R_new; /* the product is preserved */
print Y_check;

Why is this important?

As mentioned earlier, a popular algorithm in machine learning and text mining is the nonnegative matrix factorization (NMF) algorithm. Given an n x m data matrix, Y, that has nonnegative elements, you can find factors L (an n x k matrix) and R (a k x m matrix) that also have nonnegative elements, and for which L*R is the best rank-k approximation to Y subject to the constraints that L and R are nonnegative. The constraints on L and R might make us hopeful that we cannot arbitrarily form new factors (L*A) and (inv(A)*R). That is often correct: Not all invertible matrices A will result in the new factors being nonnegative!

You can prove that if A is the product of a permutation matrix and a positive diagonal matrix, then the new factors are nonnegative, but, unfortunately, there are other choices of A that also preserve nonnegativity of factors.

For example, the following rank-2 example defines L and R as nonnegative factors. The A matrix is the same as in the previous section. The inverse of A is NOT nonnegative, but the products (L*A) and (inv(A)*R) retain their nonnegative properties. Consequently, the NMF factorizations are both correct. Furthermore, you cannot transform one into the other simply by scaling and sorting the columns of L or the rows of R.

L1 = {1 2 ,
      3 8 ,
      5 4 };
R1 = {22 21 20 19,
      12 22 11 10};
Y = L1 * R1;
print Y, L1, R1;
 
A = {0.9  0.1, 0.1  0.9};
invA = inv(A);
print A, invA;     /* Note: inv(A) is not nonnegative! */
 
/* nevertheless, the products (L*A) and (inv(A)*R) can still be nonnegative */
L2 = L1*A;
R2 = inv(A)*R1;
Y_check = L2 * R2;
print Y_check, L2, R2;

This example motivates an important question: is there is a way to check whether two decompositions (L1, R1) and (L2, R2) represent the same solution up to multiplication by an invertible matrix? The answer is yes! If all the matrices are full rank and L1*R1 = L2*R2, then you can solve a least-squares system to find the matrix A that converts one solution into the other.

Here is how to find the matrix A. Assume that there is a k x k matrix A for which
L2 = L1 * A
R2 = inv(A) * R1
You can use the generalized inverse (also called the Moore-Penrose inverse) to solve for A in each equation. In the first equation, multiply on the left by the Moore-Penrose pseudoinverse for L1. In the second equation, multiply on the right by the Moore-Penrose pseudoinverse for R1. If you get the same answer from each equation, then the solutions are equivalent.

/* What if we use different algorithms and one algorithm returns
   (L1,R1) whereas the other returns (L2,R2)? Can we find an invertible 
   matrix, A, that maps one solution to the other? */
A_from_L = ginv(L1) * L2;
A_from_R = R1 * ginv(R2);
/* Do the two LS solutions give the same matrix? */
max_diff = max(abs(A_from_L - A_from_R));
print max_diff;
print A_from_L, A_from_R;     /* Note: inv(A) is not nonnegative! */

In this example, each overdetermined system produces the same matrix, so the two factorizations are equivalent up to a linear transformation.

Summary

This article explores the uniqueness of various matrix factorizations related to rank-k approximation. If you use the classical SVD algorithm, you can obtain a rank-k approximation of a data matrix that is mostly unique. To get uniqueness, you must specify how to resolve the order ambiguity and the scaling ambiguity. If a matrix has distinct singular values, you can obtain a unique rank-k decomposition that you can use to compare SVDs from different software. The factors in the SVD are orthogonal matrices, which imposes enough restrictions on the factors that you can obtain a unique representation.

Unfortunately, uniqueness is not guaranteed for other decompositions. An example is the NMF, which I discuss in another article. In the NMF, the factors are not orthogonal, but they are composed of nonnegative elements. You can resolve the scaling ambiguity and the order ambiguity by standardizing the factors. However, two algorithms might give different factors whose product is the same. In that case, you can solve a least squares system to determine whether there is an invertible linear transformation that convert one set of factors into the other set.

The post On the nonuniqueness of low-rank matrix factors appeared first on The DO Loop.

This year's article shows a relationship between pi and simple game that involves a series of coin tosses. In a short, readable paper, Jim Propp (2026) describes a simple Monte Carlo method for estimating pi by tossing a coin. In this blog post, I use SAS to demonstrate the convergence by using a Monte Carlo simulation. I also explore the distribution of the statistic that estimates pi. But before you get all your friends together and start tossing coins, let me state clearly that this estimation method is not efficient in practice! This is an example of a mathematical result that is correct, but not efficient in practice.

The game: Toss a coin until there are more heads than tails

How can you flip a fair coin to estimate pi? Imagine a stadium filled with people. Each person has a fair coin and plays the following coin-tossing game:

• Toss a coin and keep track of how many heads and tails appear.
• Stop tossing as soon as the number of heads exceeds the number of tails.
• Divide the total number of heads by the total number of tosses. Transmit this fraction to a moderator.

At the last step, each person sends their result to a moderator. (Nowadays, you would enter the numbers into a phone app!). The moderator (or app) computes the average of all the fractions of all the people in the stadium. The mathematical result by Propp is that the average of the fractions is π/4. Thus, if the moderator multiplies the average by 4, she obtains an approximation of pi.

Sequences and probabilities

Let's examine the coin-tossing process for each person in the stadium. First, notice that the total number of tosses is always odd. This is because the game stops as soon as the number of heads is one more than the number of tails. Here are some possibilities for the people in the auditorium:

• About 1/2 (or 50%) of the people in the stadium toss the coin only once. They see heads on the very first coin toss. They are finished. They transmit the fraction 1/1 to the moderator.
• 1/8 (or 12.5%) of the people toss the coin three times. Their heads-and-tails sequence is THH. They transmit the fraction 2/3 to the moderator because they saw 2 heads after 3 tosses.
• 2/32 (or 6.25%) of the people toss the coin five times. Their sequence is either TTHHH or THTHH. They transmit the fraction 3/5 to the moderator.
• 5/128 (or 3.9%) of the people toss the coin seven times. They transmit the fraction 4/7 to the moderator. Their sequence is one of the following: TTTHHHH, TTHTHHH, TTHHTHH, THTTHHH, or THTHTHH.
• In general, the fractions are all of the form (k+1) / (2k+1) for k=0, 1, 2, .... For a fair coin, Propp shows that the exact probability of the k_th fraction is $\binom{2k}{k} / (2 \cdot 4^k (k+1))$. (In practice, a coin-flipper needs to transmit only the total number of tosses (t) because the number of heads is always (t+1)/2.)

The mathematical result: The expected value of a distribution

This section is very "mathy," so feel free to skip it and proceed to the next section where I show a SAS program that simulates the coin-tossing game to estimate pi.

Let X denote the random variable that represents the ratio "number-of-heads divided by number-of-tosses" for a random trial. (A "trial" is the complete sequence of tosses.) That is, values of X are the values x_k = (k+1) / (2k+1) with the above probabilities. Recall that the expected value of the discrete random variable is the sum E[X] = Σ x_k P(X=x_k). Propp noticed that the expected value for X is the Taylor series for (1/2)arcsine(1), which is π/4.

In other words,
$\pi/4 = 1/2 \cdot (1/1) + 1/8 \cdot (2/3) + 2/32 \cdot (3/5) + 5/128 \cdot (4/7) + ...$
and the probabilities 1/2, 1/8, 2/32, 5/128, ... can be estimated by using a large number of physical coin tosses or by using a computer simulation. This is a very cool result!

This is shown graphically in the plot to the right. (Click to enlarge.) Each vertical bar is located at a ratio that represents the number of heads divided by the number of coin tosses. The height of the bar is the probability that the ratio will occur in the coin-tossing game. The red vertical line is the quantity π/4. The mathematical result is that the π/4 is the weighted sum of all the ratios.

At this point a pure mathematician is satisfied. The series converges, the result is proved, end of discussion! But an applied mathematician is interested in the rate of convergence, and on that topic, there is bad news. As a general rule, Taylor series converge fastest near the center and slowest near the boundary of their disk of convergence. The Taylor series (about 0) for arcsine(x) converges for |x| ≤ 1, which means that the series is very slow to converge when x=1. Furthermore, because it is a series of positive terms, the partial sums approach π from below, and any finite summation is less than π.

How slow is the convergence? The following approximations give you some idea:

• If you use four terms of the series (as written above), the estimate of π is 4*(0.62) = 2.483. This is not a very good approximation of pi!
• If you use 10 terms of the series, the estimate is π ≈ 4*(0.691) = 2.764, which is still a bad approximation of pi.
• If you use 100 terms of the series, the estimate is slightly more than 3: π ≈ 4*(0.757) = 3.028.
• Even if you use 100,000 terms of the series, the estimate is only good to one decimal place: π ≈ 4*(0.7845) = 3.138.

In conclusion, this series is mathematically interesting but is not an efficient way to compute pi.

A Monte Carlo simulation to estimate pi

There are several Monte Carlo simulations that estimate pi, but now we have a new method: simulate many people tossing a coin until the number of heads is greater than the number of tails. You can write a simple DATA step in SAS to perform the Monte Carlo simulation, as follows:

/* Play the coin-tossing game to estimate pi (Propp, 2026) */
%macro CoinIsHeads;
   (rand("Bernoulli", 0.5)=1)
%mend;
 
%let numPeople = 1000;
data TossCoins;
call streaminit(123);
cusum = 0;
do trialID = 1 to &numPeople;          /* each person tosses a coin */
   nTails = 0;
   nHeads = 0;
   /* toss until the number of heads is greater than the number of tails */
   do while (nTails >= nHeads);        
      if %CoinIsHeads then nHeads + 1;
      else                 nTails + 1;
   end;
   n = nTails + nHeads;                /* total number of tosses */
   trial_ratio = nHeads / n;           /* the ratio nHeads/n for each trial */
   cusum + trial_ratio;                /* cumulative sum of ratios */
   pi_estimate = 4*cusum / trialID;    /* multiply by 4 so running mean = estimate of pi */
   output;
end;
run;
 
title "Running Average of Ratio nHeads/n";
%let pi = 3.14159265358979;
proc sgplot data=TossCoins;
   scatter x=trialID y=pi_estimate;
   refline &pi / axis=y label='pi' labelpos=min lineattrs=GraphData2;
   yaxis grid label="Estimate of pi";
   xaxis grid label="Number of Trials";
run;

The graph shows the convergence of the estimate of π as you successively tabulate the results of 1,000 participants. For this random set of trials, the estimate starts large, but settles down after about 200 results. For the 1,000 participants in this trial, the final estimate is 3.14514, which is only two digits of accuracy. See Appendix 1 for a discussion of how the estimate might change if you repeat the experiment.

The number of trials is unbounded

There is a hidden inefficiency in this simulation. Notice that the innermost loop in the program is a DO-WHILE loop that simulates tossing a coin sequentially until the number of heads exceeds the number of tails. You might think that this inner loop runs a few dozen times, at most. But that's not true! A participant might get "unlucky" and start out by tossing a large number of tails, it might take hundreds or thousands of tosses before the number of heads exceeds the number of tails! You can use PROC UNIVARIATE to reveal the five largest number of trials in the simulation:

proc univariate data=TossCoins;
   var n;
   ods select ExtremeObs;
run;

The column outlined in red shows the number of coin tosses for the most unlucky participants. These unfortunate coin-flippers tossed their coins many thousands of times until the number of heads exceeded the number of tails. The most unlucky participant tossed his coin 526,919 times! Clearly, this experiment not feasible in real life! See Appendix 2 to learn more about why you should expect arbitrarily large stopping times.

Summary

Propp (2026) proves an interesting mathematical result: You can use a series of coin flips to estimate π. The method involves a mathematical game where you toss a coin until the number of heads exceeds the number of tails. You can estimate π as the average of the ratio "number-of-heads divided by number-of-tosses." Asymptotically, as many people play this game, the estimate converges to π/4, which is the expected value of the distribution of the ratio. Just multiply the average by 4 to estimate π.

However, don't start flipping coins just yet! In practice, this process converges very slowly, and an unlucky coin-flipper might spend hours or days before the sequence of heads and tails shows more heads than tails.

Appendix 1: The variance of the estimate of pi

The main body of this article includes a graph of the running total of the estimate of pi when 1,000 people play the coin-tossing game. Of course, that is only one possible realization. If you ask these same 1,000 people to repeat the experiment, the result will be different because of random chance. The following graph shows 100 different instances of playing the same game.

The values at the right side of the graph (n=1000) represent the estimates from each game. The average of those estimates is called the Monte Carlo estimate. For these trials, the Monte Carlo estimate of π is 3.1438, which is accurate to two decimal places. An empirical 95% confidence interval for that estimate is [3.090, 3.187].

Appendix 2: The stopping time of a discrete symmetric random walk

The coin-tossing game is a classic example of a one-dimensional symmetric random walk. In a symmetric random walk, we imagine an ant that starts at x=0 on a number line. With probability 1/2, the ant takes a random unit step to the left or right. What is the expected number of steps until the ant reaches x=1?

Random walks are well-studied. It is mathematically certain that the ant will eventually reach x=1, but you can prove that the expected stopping time is infinite. Returning to coin tosses, the game is guaranteed to end, but the number of coin flips before it ends is an unbounded quantity. If enough people play, some unlucky soul will pass out from exhaustion long before he reaches the stopping criterion.

The post Estimate pi by tossing coins appeared first on The DO Loop.

The data and the PCA

This article uses the Scotch whisky data set. I previously showed how to use PROC PRINCOMP in SAS to perform a classical PCA on the Whisky data. The previous article includes the following call to PROC PRINCOMP. After performing the PCA, the procedure saves the table of the top eigenvectors to a SAS data set named 'PCs'. It also saves the principal component scores to a data set named 'PC_COV_Scores'. Notice that I use a covariance-based analysis for these data because all original variables are measured on the same scale and have the same units. If you are analyzing data that have different scales or units (for example, lengths and weights), you should omit the COV option so that the correlation matrix is analyzed.

/* PCA from PROC PRINCOMP. 
   The Whisky data are defined in https://blogs.sas.com/content/iml/2026/02/23/whisky-pca.html */
%let varNames = Tobacco Medicinal Smoky Body Spicy Winey Nutty Honey Malty Fruity Sweetness Floral;
proc princomp data=Whisky N=4 cov plots=all
              out=PC_COV_Scores;   /* 'scores' = projections onto first PCs */
   var &VarNames;
   ID ID;
   ods select Eigenvectors ScoreMatrixPlot '2 by 1' '4 by 3';
   ods output Eigenvectors = PCs;
run;

Visualize the PCs: The quick way

The output from PROC PRINCOMP includes a table that shows the most important eigenvectors. The table for the whisky data is shown to the right. The eigenvectors are linear combination of the original variables. The components of the eigenvectors are weights: values that have large magnitudes indicate that the principal component is strongly correlated with that original variable. Consequently, researchers often use the eigenvector table to assign meaning to the principal components. For example, in the previous article, I made the following interpretations of the PCs:

• The first column shows how the first (and most important) PC relates to the original variables. The first PC represents high values of Body, Smoky, Medicinal, and Tobacco flavors and low values of Honey, Fruity, Sweetness, and Floral. Thus, the first PC dimension contrasts harsh "masculine" flavors from sweeter "feminine" flavors.
• The second column represents whiskies that have many taste characteristics. The second PC is associated with Body, Winey, and Honey flavors. Because most of the values in the column are positive, this component also represents a general blend of flavors.

To interpret the PCs from the table requires scanning the column and mentally noting the numbers that are largest in magnitude. This is tedious and susceptible to error. A simpler approach is to color-code the table as a heat map. Color the heat map so that cells with high values are in one color and cells with low values are in another color. You can use PROC SGPLOT to create such a heat map for the table of eigenvectors, which was written to the 'PCs' data set. The following call to PROC TRANSPOSE converts the data from wide form to long form. You can then call the HEATMAPPARM statement in PROC SGPLOT to create a heat map, as follows:

/* Use a heatmap to visualize the PCs. 
   First, transform from wide format to long format. See
   https://blogs.sas.com/content/iml/2011/01/31/reshaping-data-from-wide-to-long-format.html
*/
proc transpose data=PCs 
   out=LongPC(rename=(Col1=Value)) name=PC;                    
   by Variable notsorted;       /* for each row */
   var Prin1-Prin4;             /* make a row for these variables */
run;
 
title "Principal Components";
title2 "Unsymmetric Range for Gradient Legend";
ods graphics / width=360px height=480px;
proc sgplot data=LongPC;
   heatmapparm x=PC y=Variable colorresponse=Value / colormodel=AltTwoColor outline;
   yaxis reverse;
   label PC="Principal Components" Variable="Flavor Characteristics";
run;

This visualization uses red for high values and blue for low values.

• The first column shows that the largest (positive) contributions for the first PC come from the Medicinal, Smoky, and Body variables. Those cells are in red. The first PC has "negative" amounts of Honey, Fruity, Sweetness, and Floral. Those cells are in blue.
• The second column shows that the second PC is formed by using positive contributions from the Winey, Body, and Honey variables, and "negative" amounts of the Medicinal and Floral variables.
• The third PC has a large contribution from the Spicy variable. The fourth PC that has a large contribution from the Nutty variable.

This graph is sufficient for many PCA analyses, but it has a problem. I do not like the fact that the range of the color scale is not symmetric about the value 0. In this heat map, there are blue cells that represent positive values. The darkest red corresponds to the value +0.85, whereas the darkest blue corresponds to -0.38.

I think the graph would be improved and more accurate if the neutral color (white) corresponded to 0, and the positive and negative ranges were equal in length. Then you could unequivocally say, "red cells are positive, and blue cells are negative." The code in the next section shows how to accomplish this modification. The resulting heat map follows:

This second heat map is easier to read because now red represents positive contributions, white represents zero, and blue represents negative contributions. For each column, you can use the pattern of dark reds and blues to interpret how the PC is formed as a linear combination of the original variables.

A symmetric range attribute map for diverging color ramps

A two-color ramp that passes through a third "neutral" color is called a diverging color palette. They are most informative when the neutral color corresponds to a meaningful value. Sometimes the neutral color corresponds to a baseline value or an average value. In this case, I want the neutral color corresponds to 0. I also want the range of colors to correspond to a symmetric interval about 0 so that dark shades or red and blue represent the same magnitudes but different signs.

For the heat map in the previous section, the maximum value in the table is about 0.85, and the minimum value is about -0.38. Thus, the middle value (0.235) is used for the neutral color. Instead, I want to define the color range to be the symmetric interval [-0.85, 0.85] so that the midpoint (0) is used for the neutral color.

In SAS, you can accomplish this by using a range attribute map. A range attribute map is a data set that assigns the range and colors for a color ramp. When you define a range attribute map for a variable that has both positive and negative values, you can use two special values to ensure that the range is symmetric about 0:

• Use the keyword _NEGMAXABS_ for the minimum value. This keyword is equivalent to the expression –max(abs(MIN), abs(MAX)).
• Use the keyword _MAXABS_ for the maximum value. This keyword is equivalent to the expression max(abs(MIN), abs(MAX)).

The following SAS DATA step creates the range attribute map. It uses a red-blue diverging palette of colors from the ColorBrewer system.

/* create a range attribute data set so that the value range is symmetric, such as 
   [-0.85, 0.85] with 0 in the middle. There is a trick to doing this:
   Use "_NEGMAXABS_" = –max(abs(MIN), abs(MAX))  for the MIN variable and 
   use "_MAXABS_" = max(abs(MIN), abs(MAX))  for the MAX variable.
   Use a 5-color diverging blue-white-red palette from the ColorBrewer system:
   CX0571B0 CX92C5DE CXF7F7F7 CXF4A582 CXCA0020 
*/
data SymRangeAttrs;
retain ID "SymRange";
length min max $12 color altcolor colormodel1-colormodel5 $15;
input min max color altcolor colormodel1-colormodel5;
datalines;
_NEGMAXABS_ _MAXABS_  .  .  CX0571B0 CX92C5DE CXF7F7F7 CXF4A582 CXCA0020
;
 
title "Principal Components";
title2 "Symmetric Range";
ods graphics / width=360px height=480px;
proc sgplot data=LongPC RATTRMAP=SymRangeAttrs;  /* <== DEFINE ATTR MAP DATA SET HERE */
   heatmap x=PC y=Variable / colorresponse=Value outline
                             RATTRID=SymRange;                 /* <== NAME MAP HERE */
   yaxis reverse;
   label PC="Principal Components" Variable="Flavor Characteristics";
run;

The resulting heat map is shown in the previous section.

Visualize the principal component scores

In addition to visualizing the principal components themselves, you can construct a heatmap that visualizes the PC scores for each observation. Normally, scatter plots are used to visualize the PC scores. For example, the scatter plot to the right shows a scatter plot that shows the scores for the first two PCs for the 86 whiskies in the data. (Click to enlarge.)

The scatter plot shows you which observations (whiskies) have high or low values when projected onto the principal components. Typically, you would look at multiple scatter plots to fully understand how each observations relates to the PC factors. However, if you do not have too many observations, you can summarize that information in a single heat map. With a single heat map, you can see the PC scores for each observation simply by scanning across a row. You can find the high and low values for each PC simply by scanning down each column. You could sort the data by the scores onto the first PC, but I have not done that here.

Because the heat map shows one row for each observation, this visualization is not appropriate when there are many observations. However, the heat map can be useful for visualizing a subset of the data. The following statements visualize a subset of 20 whiskies. PROC TRANSPOSE converts the scores from wide (four columns) to long format. The value for each dimension is shown by using the symmetric range attribute map, as explained in the previous section:

/* Now use a heatmap to visualize the projection of centered data onto dominant PCs. */
proc transpose data=PC_COV_Scores(keep=Distillery RowID Selected Prin1-Prin4)
               out=LongScore(rename=(Col1=Value)) name=PC;
   where selected=1;            /* only display a subset of the data */
   by Distillery notsorted;     /* for each row */
   var Prin1-Prin4;             /* make a row for these variables */
run;
 
ods graphics / width=480px height=640px ;
title "Principal Component Scores for 20 Whiskies";
proc sgplot data=LongScore RATTRMAP=SymRangeAttrs;  /* <== DEFINE ATTR MAP DATA SET HERE */
   heatmap x=PC y=Distillery / colorresponse=Value outline
                             RATTRID=SymRange;                 /* <== NAME MAP HERE */
   yaxis reverse;
   label PC="Principal Components";
run;

With this heat map, you can find the whiskies whose flavor profiles are positively or negatively correlated with each PC. If you like whiskies that are full-bodied, medicinal, and smoky (the primary attributes captured by the first PC), then scan down the first column and find all the dark red whiskies. Examples include Laphroaig, Lagavulin, and Ardbeg.

The "opposite" flavor profile for the first PC is floral, sweet, and fruity. To find whiskies with those flavors, scan down the first column and find all the dark blue whiskies. The deepest blue is found for Auchentoshan, GlenMoray, Glengoyne.

Repeat this process for other flavor characteristics that are captured by a principal component: Red for whiskies that have a positive correlation and blue for whiskies that are negatively correlated.

Notice that each row of this heat map enables you to visualize the four-dimensional scores for each whisky. This is an advantage of scatter plots, which show two coordinates at a time. You would need six scatter plots to show all pairwise combinations of PC scores for this analysis, but this one heat map contains the same information. Of course, the human eye cannot distinguish shades very well, so this display is best for providing qualitative information. The scatter plots are better at showing small differences in the PC scores.

Summary

Statistical procedures in SAS provide many graphs to help you visualize your analyses. The PRINCOMP procedure can produce a large number of informative graphs. However, this article shows that it is helpful to use a heat map to visualize the table of eigenvectors. The heat map is most useful when it uses a range attribute map that is symmetric about 0. You can use a similar heat map to visualize the scores of a subset of the data. This provides a single graph that reveals qualitative differences between whiskies.

The post Two new visualizations of a principal component analysis appeared first on The DO Loop.

Scotch whisky (spelled without an 'e') is a popular drink that makes up a multi-billion dollar industry. Scotch whisky accounts for almost 75% of Scotland's food and drink exports! Poets have lauded it. (See Appendix 2.) On television and in movies, amber-brown whiskies are consumed by upper-class characters who have (or crave) power, prestige, and wealth. But for those of us who are not so wealthy, it is scary to spend $50 or more on a bottle of whisky that we might not enjoy.

Fortunately, statistics can help us to identify whiskies that we might like or dislike based on others that we have tasted. Suppose a friend offers a "wee dram" of whisky from his private stash. Maybe we like it, maybe we don't. If we like it, we might want to find other brands that taste similar. If we hate it, we'd like to identify other brands that we should avoid.

There are several statistical methods that can identify whiskies that have a similar set of flavor characteristics, including classical techniques such as cluster analysis and factor analysis. In this article, I will use classical principal component analysis (PCA) to analyze the flavor profiles of 86 Scotch whiskies. I will show which flavor characteristics are important for differentiating whiskies. A future article compares the PCA to a non-negative matrix factorization (NMF) of the same data.

The data and analyses are motivated by an article by Young, Fogel, and Hawkins (2006, "Clustering scotch whiskies using non-negative matrix factorization", SPES Newsletter). The data are from the book Whisky Classified, Choosing Single Malts by Flavor by David Wishart (2002, revised 2018), who classified Scotch whisky according to 12 flavors: Tobacco, Medicinal, Smoky, Body, Spicy, Winey, Nutty, Honey, Malty, Fruity, Sweetness, and Floral.

Full disclosure: I'm a statistician, not a Scotch connoisseur! My experience tasting whisky is very limited. I have only tasted two Scotch whiskies in my life: I liked Glenfiddich, but I did not like Laphroaig.

The Scotch whisky data

Appendix 1 defines the data used by Young, Fogel, and Hawkins (2006). There are 86 brands of whisky (the rows) and 12 flavors (the columns) plus some ID variables. The following table shows the flavor profiles for a subset of 20 whiskies:

Each whisky is assigned values 0-4 for each flavor, which indicates the absence of a flavor (0) or an intense expression of that flavor (4). Let's look at the flavor profiles for the two whiskies that I have tasted:

• Glenfiddich is light-bodied and a little smoky. Its predominant characteristic is sweetness, and it has secondary characteristics of malty, fruity, and floral.
• Laphroaig is full-bodied and very smoky and medicinal. It has a little sweetness with hints of tobacco, winey, nutty, and malty.

Principal components: Extracting flavor profiles

Flavors are imparted during various steps in the distillation process such as malting, fermentation, and maturing the whisky in oak casks. They are interconnected and cannot be independently controlled by the blender. Scotch whisky can only contain water, grain, and yeast, so the Master Blender cannot add nuts or honey or fruit to affect the taste of the finished product.

Because the flavors are correlated, a principal component analysis (PCA) is a classical way to identify combinations of flavors that most explain the variance in the 12-dimensional flavor space. A common use of PCA is to reduce the dimensionality of the problem from 12 flavor characteristics to a smaller number of "principal components," which are linear combinations of the flavors. Following Young, Fogel, and Hawkins (2006), I choose four PCs, which explain almost 70% of the variability in the data. In SAS, you can use the following call to PROC PRINCOMP to perform a PCA and to generate many tables and graphs that explain the analysis. Because the flavors are measured on the same scale, you can specify the COV option to use the covariance matrix for the analysis:

/* PCA from PROC PRINCOMP */
%let varNames = Tobacco Medicinal Smoky Body Spicy Winey Nutty Honey Malty Fruity Sweetness Floral;
proc princomp data=Whisky N=4 cov plots=all
              out=PC_COV_Scores;   /* 'scores' = projections onto first PCs */
   var &VarNames;
   ID ID;
   ods select Eigenvectors ScoreMatrixPlot '2 by 1' '4 by 3';
   ods output Eigenvectors = PCs;
run;

The procedure outputs many tables and graphs, which are explained in the next sections.

Interpreting the principal components

The most important output from a PCA is the set of eigenvectors, which are used as the basis vectors. The following table shows four eigenvectors, which are called the first four principal components (PCs).

• The first column shows how the first (and most important) PC relates to the original variables. The first PC represents high values of Smoky, Medicinal, and Tobacco flavors and low values of Honey, Fruity, Sweetness, and Floral. Thus, the first PC dimension contrasts harsh "masculine" flavors from sweeter "feminine" flavors.
• The second column represents whiskies that have many taste characteristics. The second PC is associated with full-bodied, Winey, and Honey flavors. Because most of the values in the column are positive, this component also represents a general blend of flavors.
• The third column represents whiskies that are Spicy. It also contrasts whiskies that are high in Fruity and Floral flavors from those that are high in Sweetness and Winey.
• The fourth column represents whiskies that are Nutty. It contrasts them with Winey and Sweetness flavors.

A PCA is a "contrast-revealing" algorithm. In a PCA analysis, the eigenvectors are orthogonal to each other. This forces the principal components to contain both positive and negative weights. This allows PCA to create contrasts between flavors, but it is unintuitive to think about a whisky having a 'negative' amount of a flavor! A subsequent article shows how the nonnegative matrix factorization (NMF) drops the orthogonality constraintso that it can build flavor profiles additively.

Visualizing the principal components: Component pattern plots

These PCs are easier to visualize if we view the "component pattern" plots, which I have described in a previous article. They show the correlations between the principal components and the original variables. The following graph shows the correlations between the original variables and the first two PCs.

Notice that the horizontal axis separates the "masculine" flavors from the sweet and fruity and floral flavors. Similarly, the vertical axis is positively correlated with most of the flavors, but is highly correlated with Winey, Body, and Honey.

You can show a similar plot for the third and fourth PCs:

Interpreting the PC scores: Identifying similar and dissimilar whiskies

Because the PCs differentiate between various flavors, it is useful to project the original observations onto the basis of the PCs. This is called "scoring," and the resulting scatter plot are called score plots. It is useful to identify the points by the name of the distillery. However, since distillery names can be quite long (eg, "GlenDeveronMacduff"), I have used distillery names for 20 selected whiskies and have used row numbers to identify the remaining whiskies. The score plot for first two PCs is shown below:

In these plots, the origin (0,0) represents a hypothetical whisky that has the "average" flavor profile. Whiskies in the upper-right quadrant have greater-than-average flavors in both components, whereas whiskies in the lower-left quadrant have less-than-average flavor in both components.

This is the most important plot for a beginner to study. Why? Because the whiskies to the right are "strong" on the first PC, meaning that they are full-bodied with strong accents of Smoky, Medicinal, and Tobacco flavors, and low values of Honey, Fruity, Sweetness, and Floral. The Laphroaig whisky, which I have tasted and did not like, is far to right. Thus, I know to steer clear of similar whiskies such as Ardberg and Lagavulin.

In a similar way, the whiskies at the top of the plot (Macallan and Glendronach) represent whiskies that are full-bodied with predominately Winey and Honey flavors. Notice that Macallan and Glendronach have approximately 0 score on the horizontal axis, which means that they are not strong associated with the flavors represented by the first PC. The Glenfiddich whisky, which I like, is at the bottom of this graph. That means it is light-bodied and does not have Winey and Honey flavors. It is left-of-center on the horizontal scale, so it does not have the strong "masculine" flavors that I dislike about Laphroaig.

The score plot onto the third and fourth PCs provide additional information, primarily about the Spicy and Nutty flavors.

Whiskies on the right side of this plot include GlenGarioch and Talisker, which are spicy. To the left are Tomore and Laphroaig, which have no Spicy flavor. High on the vertical axis are Edradour and Longmorn, which are nutty. At the bottom of the plot are Linkwood and GlenDeveronMacduff, which have no Nutty flavor.

For the whiskies I have tasted, Laphroaig is to the left and in the middle, indicating that it has no Spicy flavor and a moderate amount of Nutty flavor. Glenfiddich is in the lower-left quadrant, so it has negligible Spicy and Nutty flavors.

Summary

The purpose of this article is to introduce 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), who used nonnegative matrix factorization (NMF). I describe that analysis in a subsequent article. However, the strengths of the NMF are best appreciated when compared with a conventional principal component analysis (PCA), which is a classical way to identify combinations of variables that explain the variance in the 12-dimensional flavor space. This article shows how a PCA can help someone who is ignorant about single-malt whiskies to understand which distilleries have similar flavor profiles.

Appendix 1: Scotch whisky data

Define a SAS data set that contains the original data from Wishart (2004) that was analyzed by Young, Fogel, and Hawkins (2006). Then add a new variable, SELECTED, for 20 whiskies that are either popular brands (eg, Glenlivet, Glenfiddich, Macallan,...) or have extreme scores in a PCA of the data.

/* Perform a principal component analysis (PCA) of the 
Scotch Whisky data set in Young, Fogel, and Hawkins 
(SPES, 2006), who analyzed the data by using a nonnegative 
matrix factorization (NMF). The article is available at 
https://www.niss.org/sites/default/files/ScotchWhisky.pdf
 
I downloaded the data from  
https://www.niss.org/sites/default/files/ScotchWhisky01.txt
and corrected a few typos:
- ID=12: Replaced 'Belvenie' with 'Balvenie' 
- ID=56: Replaced 'Laphroig' with 'Laphroaig' 
*/
data Whisky_Orig;
infile datalines dsd dlm=','; 
length Distillery $20;
input
RowID Distillery Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey Nutty Malty Fruity Floral;
datalines;
01,Aberfeldy,2,2,2,0,0,2,1,2,2,2,2,2
02,Aberlour,3,3,1,0,0,4,3,2,2,3,3,2
03,AnCnoc,1,3,2,0,0,2,0,0,2,2,3,2
04,Ardbeg,4,1,4,4,0,0,2,0,1,2,1,0
05,Ardmore,2,2,2,0,0,1,1,1,2,3,1,1
06,ArranIsleOf,2,3,1,1,0,1,1,1,0,1,1,2
07,Auchentoshan,0,2,0,0,0,1,1,0,2,2,3,3
08,Auchroisk,2,3,1,0,0,2,1,2,2,2,2,1
09,Aultmore,2,2,1,0,0,1,0,0,2,2,2,2
10,Balblair,2,3,2,1,0,0,2,0,2,1,2,1
11,Balmenach,4,3,2,0,0,2,1,3,3,0,1,2
12,Balvenie,3,2,1,0,0,3,2,1,0,2,2,2
13,BenNevis,4,2,2,0,0,2,2,0,2,2,2,2
14,Benriach,2,2,1,0,0,2,2,0,0,2,3,2
15,Benrinnes,3,2,2,0,0,3,1,1,2,3,2,2
16,Benromach,2,2,2,0,0,2,2,1,2,2,2,2
17,Bladnoch,1,2,1,0,0,0,1,1,0,2,2,3
18,BlairAthol,2,2,2,0,0,1,2,2,2,2,2,2
19,Bowmore,2,2,3,1,0,2,2,1,1,1,1,2
20,Bruichladdich,1,1,2,2,0,2,2,1,2,2,2,2
21,Bunnahabhain,1,2,1,1,0,1,1,1,1,2,2,3
22,Caol Ila,3,1,4,2,1,0,2,0,2,1,1,1
23,Cardhu,1,3,1,0,0,1,1,0,2,2,2,2
24,Clynelish,3,2,3,3,1,0,2,0,1,1,2,0
25,Craigallechie,2,2,2,0,1,2,2,1,2,2,1,4
26,Craigganmore,2,3,2,1,0,0,1,0,2,2,2,2
27,Dailuaine,4,2,2,0,0,1,2,2,2,2,2,1
28,Dalmore,3,2,2,1,0,1,2,2,1,2,3,1
29,Dalwhinnie,2,2,2,0,0,2,1,0,1,2,2,2
30,Deanston,2,2,1,0,0,2,1,1,1,3,2,1
31,Dufftown,2,3,1,1,0,0,0,0,1,2,2,2
32,Edradour,2,3,1,0,0,2,1,1,4,2,2,2
33,GlenDeveronMacduff,2,3,1,1,1,1,1,2,0,2,0,1
34,GlenElgin,2,3,1,0,0,2,1,1,1,1,2,3
35,GlenGarioch,2,1,3,0,0,0,3,1,0,2,2,2
36,GlenGrant,1,2,0,0,0,1,0,1,2,1,2,1
37,GlenKeith,2,3,1,0,0,1,2,1,2,1,2,1
38,GlenMoray,1,2,1,0,0,1,2,1,2,2,2,4
39,GlenOrd,3,2,1,0,0,1,2,1,1,2,2,2
40,GlenScotia,2,2,2,2,0,1,0,1,2,2,1,1
41,GlenSpey,1,3,1,0,0,0,1,1,1,2,0,2
42,Glenallachie,1,3,1,0,0,1,1,0,1,2,2,2
43,Glendronach,4,2,2,0,0,2,1,4,2,2,2,0
44,Glendullan,3,2,1,0,0,2,1,2,1,2,3,2
45,Glenfarclas,2,4,1,0,0,1,2,3,2,3,2,2
46,Glenfiddich,1,3,1,0,0,0,0,0,0,2,2,2
47,Glengoyne,1,2,0,0,0,1,1,1,2,2,3,2
48,Glenkinchie,1,2,1,0,0,1,2,0,0,2,2,2
49,Glenlivet,2,3,1,0,0,2,2,2,1,2,2,3
50,Glenlossie,1,2,1,0,0,1,2,0,1,2,2,2
51,Glenmorangie,2,2,1,1,0,1,2,0,2,1,2,2
52,Glenrothes,2,3,1,0,0,1,1,2,1,2,2,0
53,Glenturret,2,3,1,0,0,2,2,2,2,2,1,2
54,Highland Park,2,2,3,1,0,2,1,1,1,2,1,1
55,Inchgower,1,3,1,1,0,2,2,0,1,2,1,2
56,Isle of Jura,2,1,2,2,0,1,1,0,2,1,1,1
57,Knochando,2,3,1,0,0,2,2,1,2,1,2,2
58,Lagavulin,4,1,4,4,1,0,1,2,1,1,1,0
59,Laphroaig,4,2,4,4,1,0,0,1,1,1,0,0   
60,Linkwood,2,3,1,0,0,1,1,2,0,1,3,2
61,Loch Lomond,1,1,1,1,0,1,1,0,1,2,1,2
62,Longmorn,3,2,1,0,0,1,1,1,3,3,2,3
63,Macallan,4,3,1,0,0,2,1,4,2,2,3,1
64,Mannochmore,2,1,1,0,0,1,1,1,2,1,2,2
65,Miltonduff,2,4,1,0,0,1,0,0,2,1,1,2
66,Mortlach,3,2,2,0,0,2,3,3,2,1,2,2
67,Oban,2,2,2,2,0,0,2,0,2,2,2,0
68,OldFettercairn,1,2,2,0,1,2,2,1,2,3,1,1
69,OldPulteney,2,1,2,2,1,0,1,1,2,2,2,2
70,RoyalBrackla,2,3,2,1,1,1,2,1,0,2,3,2
71,RoyalLochnagar,3,2,2,0,0,2,2,2,2,2,3,1
72,Scapa,2,2,1,1,0,2,1,1,2,2,2,2
73,Speyburn,2,4,1,0,0,2,1,0,0,2,1,2
74,Speyside,2,2,1,0,0,1,0,1,2,2,2,2
75,Springbank,2,2,2,2,0,2,2,1,2,1,0,1
76,Strathisla,2,2,1,0,0,2,2,2,3,3,3,2
77,Strathmill,2,3,1,0,0,0,2,0,2,1,3,2
78,Talisker,4,2,3,3,0,1,3,0,1,2,2,0
79,Tamdhu,1,2,1,0,0,2,0,1,1,2,2,2
80,Tamnavulin,1,3,2,0,0,0,2,0,2,1,2,3
81,Teaninich,2,2,2,1,0,0,2,0,0,0,2,2
82,Tobermory,1,1,1,0,0,1,0,0,1,2,2,2
83,Tomatin,2,3,2,0,0,2,2,1,1,2,0,1
84,Tomintoul,0,3,1,0,0,2,2,1,1,2,1,2
85,Tomore,2,2,1,0,0,1,0,1,2,1,0,0
86,Tullibardine,2,3,0,0,1,0,2,1,1,2,2,1
;
 
/* For later analysis, identify some of the most popular single-malt 
   Scotch whiskies (eg, Glenlivet, Glenfiddich, Macallan,...) and 
   some that have extreme scores in a PCA of the data.
*/
data Whisky / view=Whisky;
length ID $20;
set Whisky_Orig;
selected = 0;
BestSeller = 0;
if Distillery in (
   'Glenlivet' 'Glenfiddich' 'Macallan' 'Glenmorangie' 'Balvenie' 
   'Laphroaig' 'Aberlour' 'Lagavulin' 'Ardbeg' 'Talisker' 
   ) then BestSeller=1;
/* also select a few whiskies that have extreme values in a PC */
if BestSeller | RowID in (07 11 32 33 35 43 60 62 82 85) then 
   selected = 1;
if selected then 
   ID = Distillery;
else 
   ID = put(RowID, Z2.);
run;
 
%let varNames = Tobacco Medicinal Smoky Body Spicy Winey Nutty Honey Malty Fruity Sweetness Floral;
proc print data=Whisky;
   where selected;
   ID ID;
   var &VarNames;
run;

Appendix 2: Robert Burns' "Scotch Drink" poem

Robert Burns was one of Scotland's most famous poets. One of his poems, "Scotch Drink" (1785), is an ode to whisky. It contains more than 20 verses, but here is one amusing verse:

O thou, my muse! guid auld Scotch drink!
Whether thro' wimplin worms thou jink,
Or, richly brown, ream owre the brink,
In glorious faem,
Inspire me, till I lisp an' wink,
To sing thy name!

If the 18th century spelling is confusing, here is a modern translation that uses conventional spelling. The phrase "wimplin worms thou jink" refers to the distilling process, which uses winding coils of tubing. "To jink" means to change direction.

O you, my muse! Good old Scotch drink!
Whether through winding worms you jink,
Or, richly brown, foam over the brink,
In glorious foam,
Inspire me, until I lisp and wink,
To sing your name!

The post Comparing flavor characteristics of Scotch whiskies: A principal component analysis appeared first on The DO Loop.

• The Pascal matrix is a symmetric positive definite matrix with a known Cholesky decomposition and inverse.
• A Toeplitz matrix is symmetric (and banded) and can be constructed so that it is positive definite.
• A special Toeplitz matrix, called the first-order autoregressive (AR(1)) correlation matrix, has explicit eigenvalues.

I recently learned about another symmetric and positive definite matrix, called the Min Matrix. The Min Matrix of size n is the matrix defined by M[i,j] = min(i,j) for 1 ≤ i,j ≤ n. The general form of the Min Matrix is shown below:

The Min Matrix has some wonderful properties. It is easy to generate. It is symmetric and positive definite (SPD). The inverse and Cholesky roots are also integer valued. The eigenvalues have explicit formulas.

Creating the Min Matrix in SAS IML

In the SAS IML language, you can create the Min Matrix without using any loops. By using the built-in ROW and COL functions in conjunction with the elementwise minimum operator (><), you can generate the matrix in a single line of code:

proc iml;
/* return the Min Matrix of size n. M[i,j] = min(i,j) */
start MinMatrix(n);
   M = row(I(n)) >< col(I(n));
   return M;
finish;
 
n = 5;
M = MinMatrix(n);
print M[L='Min Matrix'];

The Cholesky root and the determinant

A property of the Min Matrix is that its Cholesky factor is known and integer-valued. The Cholesky root of an SPD matrix, M, is an upper-triangular matrix, U, such that M = U^TU. (In PROC IML, the ROOT function returns the upper triangular factor.) For the Min Matrix, the Cholesky root contains all 1s on and above the main diagonal. Not only is this a cool result, but the root immediately tells us that the determinant of M is always 1. This follows because det(M) = det(U^T)det(U) = 1.

U = root(M);
print U[L='Cholesky Root'];

The tridiagonal inverse

The Min Matrix is a dense matrix, but its inverse is a tridiagonal matrix. The inverse matrix has 2s on the main diagonal (except for the very last element, which is 1) and -1s on the superdiagonal and subdiagonal. This inverse matrix is closely related to the matrix for a discrete second-difference operator, except for the (n,n)th element.

M_Inv = inv(M);
print M_Inv[L='Inverse Matrix'];

Eigenvalues

The structure of the inverse matrix enables use to find a formula for its eigenvalues. The k_th eigenvalue for M^-1 is $\mu_k = 4 \sin^2\left( \frac{(2k-1)\pi}{2(2n+1)} \right)$

For any invertible matrix, the eigenvalues of the matrix and its inverse are reciprocals. Therefore, the eigenvalues of the Min Matrix of size n are given by the following expression: $\lambda_k = \frac{1}{4 \sin^2\left( \frac{(2k-1)\pi}{2(2n+1)} \right)}$

You can verify this formula by comparing the analytical results to the numerical computation from the EIGVAL function:

eigenvals = eigval(M);    /* numerical eigenvalues */
/* compare to the analytical formula */
k = T(1:n);
pi = constant('pi');
check_eig = 1 / (4 * sin( ((2*k-1)*pi) / (2*(2*n+1)) )##2);
Diff = eigenvals - check_eig;
print eigenvals check_eig Diff;

From the analytical expression for the eigenvalues, you can deduce that the Min Matrix is ill-conditioned for large n. As n increases, the largest eigenvalue (k=1) grows like O(n²). In contrast, the smallest eigenvalue (k=n) approaches the value 1/4. The condition number of a matrix is the ratio of the largest to smallest eigenvalue, so the condition number of the Min Matrix grows quadratically with n. This makes it an a good test case for matrix computations.

Summary

The Min Matrix is a good test matrix to know about. It is a covariance matrix because it is symmetric and positive definite. It is easy to construct, it has a simple Cholesky root, it has an integer-valued tridiagonal inverse, and its eigenvalues are known exactly.

The post Properties of the "Min Matrix" appeared first on The DO Loop.

The range is also a measure of dispersion, although it is not robust to outliers. I have previously shown how to use the range to form a consistent estimator of σ for normally distributed data. The scaling factor depends on the expected value of the range for a normal sample of size n, which is called d₂(n).

When the sample size is small (n ≤ 25), quality control engineers use the range to estimate the standard deviation. You can find values of d₂(n) in a table or use the previous article to compute values by using the D2 function in SAS/QC software or by solving an integral by using SAS IML.

I recently saw an approximation to d₂(n) by using an asymptotic (large n) formula. An asymptotic formula approximates the expected value of the range in a normal sample when n is large. I was surprised to discover that the asymptotic approximation is not a good approximation for d₂(n) for moderately large samples, such as n < 1000. This post compares the asymptotic formula to the exact value of d₂(n) and compares both methods with a Monte Carlo simulation.

The Expected Range and Extreme Value Theory

The range of a sample is the difference between the maximum and the minimum values: R = X_(n) - X₍₁₎. For a symmetric distribution like the normal distribution, the expected value of the range is simply twice the expected value of the maximum: E[R] = 2·E[X_(n)].

I have previously shown that as n becomes large, the distribution of the maximum of a normal sample approaches a Gumbel distribution. The graph of the expected value of the maximum is shown to the left for large sample sizes. You can write down formulas that describe the asymptotic behavior of this function. I used formulas in David and Nagaraja (2004, Order Statistics, 3rd Ed.). First, define the following functions:

• $s_n = \sqrt{2 \log n}$
• $b_n = s_n - \frac{\log(\log n) + \log(4\pi)}{2s_n}$
• $a_n = \frac{1}{s_n}$

In these formulas, LOG is the natural logarithm. By using these expressions, you can construct several asymptotic approximations for d₂(n). Each approximation adds one additional term to the previous approximation:

1. 1-term approximation: 2s_n. This is the simplest (and crudest) estimate.
2. 2-term approximation: 2b_n. This a better approximation, but it tends to underestimate d₂(n).
3. 3-term approximation: 2(b_n + γa_n), where γ ≈ 0.5772 is the Euler-Mascheroni constant. This accounts for the fact that the mean of the limiting Gumbel distribution is γ. This is the best approximation of the three.

These are large-n approximations for d₂(n). The following SAS DATA step compares the exact values of d₂(n) to these asymptotic formulas.

/* d2(n) is the scale factor to use to adjust the range statistic so that 
   it is a consistent estimator of the standard deviation for normal data.
   Some references use asymptotic approximations for d2(n).
   Define the formulas 
   sn = sqrt(2*log(n));
   bn = sn - (log(log(n)) + log(4*pi)) / (2*sqrt(2*log(n)));
   an = 1 / sn;
 
   Then the following are 1-term, 2-term, and 3-term asymptotic approximations for d2(n):
 
   Number of Terms    Formula
   1 term             2 * sn
   2 terms            2 * bn
   3 terms            2 * (bn + gamma*an), where gamma~0.57721... is the Euler-Mascheroni constant
*/
data k_n;
/* Compute or look up the value of d2(n) in a table for certain values of n.
   Compare those exact values to several asymptotic (large n) approximations. */
array nA[5]   _temporary_ (2,      50,     100,    500,    1000);
array d_nA[5] _temporary_ (1.1284, 4.4982, 5.0151, 6.0734, 6.4829);
pi = constant('pi');
gamma = constant('Euler');  /* 0.57721... */
do i = 1 to dim(nA);
   n  = nA[i];
   d2_n = d_nA[i];
   sn = sqrt(2*log(n));
   bn = sn - (log(log(n)) + log(4*pi)) / (2*sqrt(2*log(n)));
   an = 1 / sn;
   Asymp1 = 2 * sn;
   Asymp2 = 2 * bn;
   Asymp3 = 2 * (bn + gamma*an);
   output;
end;
drop pi gamma i sn bn an;
run;
 
proc print data=k_n noobs; 
   format _all_ best5.; 
   ID n; 
run;

I have seen researchers use these asymptotic approximations for values of n as small as 50 or 100. But, as you can see, the approximations are not close to d₂(n) for those values of n. In fact, David and Nagaraja (2003) warn that these asymptotic approximations converge quite slowly. Even at n = 1000, the approximations are not very close to the true value. I do not recommend these asymptotic approximations for small or moderate sample sizes.

Testing the Estimators: A Monte Carlo Simulation

You can run a Monte Carlo simulation to estimate d₂(n), which is the expected range of a normal sample of size n. The following SAS DATA step generates 1,000 samples of size n = 100 from a standard normal distribution. For each sample, PROC MEANS calculates the range. In a second DATA step, the range is divided by three different factors: the exact d₂(100), the 2-term approximation, and the 3-term approximation. These "adjusted ranges" are the various estimates for the standard deviation, which is 1 for this simulation. I do not show the 1-term adjusted range because it is so bad.

/* simulate 1000 random samples of size n=100 from N(0,1) */
%let N = 100;
%let numSamples = 1000;
data SimNormal;
call streaminit(12345);
do sampleID = 1 to &numSamples;
   do i = 1 to &n;
      x = rand("Normal");
      output;
   end;
end;
keep sampleID x;
run;
 
/* compute the range for each random sample */
proc means data=SimNormal noprint;
    by sampleID;
    var x;
    output out=RangeOut RANGE=Range;
run;
 
/* adjust the range statistics by using d2(100), Asymp2, or Asymp3 formula */
data AdjustRange;
pi = constant('pi');
gamma = constant('Euler');  /* 0.5772156649 */
set RangeOut(rename=(_FREQ_=n));
d2_n = 5.0151;
sn = sqrt(2*log(n));
bn = sn - (log(log(n)) + log(4*pi)) / (2*sqrt(2*log(n)));
an = 1 / sn;
Asymp2 = 2 * bn;
Asymp3 = 2 * (bn + gamma*an);
 
AdjRange_d2     = 1/d2_n   * Range;
AdjRange_2Term  = 1/Asymp2 * Range;
AdjRange_3Term  = 1/Asymp3 * Range;
keep SampleID Range AdjRange_d2 AdjRange_2Term AdjRange_3Term;
label AdjRange_d2="d2" AdjRange_2Term="2-Term Asymptotic" AdjRange_3Term="3-Term Asymptotic";
run;
 
/* Print the Monte Carlo estimates of the adjusted range, which should be close to sigma=1.
   Note that the 95% CI for the mean includes 1 for d2(100), but not for the asymptotic approx. 
*/
proc means data=AdjustRange Mean CLM ndec=3 nolabels;
   var AdjRange_d2 AdjRange_2Term AdjRange_3Term;
run;

Since the data were generated from N(0,1), a consistent estimator should produce an average value close to 1.0. The results show that the Monte Carlo estimate is consistent if you use the d₂(100) value to scale the data. For the asymptotic formulas, you get biased estimates. The two-term asymptotic approximation overestimates the standard deviation whereas the two-term asymptotic approximation underestimates it. These results are reproducible: they do not change if you use additional Monte Carlo trials.

Summary

The expected range of normal data is a quantity known as d₂(n). This article shows several asymptotic approximations. The asymptotic approximations provide a simple expression that you can use to predict the range in a very large random sample of normal data. However, these estimates converge quite slowly. In quality control, the range is used to estimate the standard deviation of normal data when n is small. You should not use the asymptotic approximations of d₂(n) for that purpose.

The post The asymptotic expected value of the range for normal data appeared first on The DO Loop.

Years ago, I wrote an article about the "trap and cap" programming technique. The idea is that programmers should "trap" inputs to functions (like SQRT, LOG, and QUANTILE functions) to avoid domain errors. In addition, when visualizing a function's range, you should "cap" the output to improve graphs of functions (such as LOG, 1/X, or TAN) that contain singularities.

Each of these is an example of what computer scientists call clipping or clamping. Mathematicians call it truncation to an interval. An example is shown to the right. In the example, input values are in the interval [-2, 2]. The "clip" function ensures that the values are clipped into the interval [-1, 1]. For references, the identity function is shown as a dashed line in the background.

This article shows how to write SAS logic for clipping values in the DATA step and in SAS IML.

Clipping in the DATA step: Nested MAX and MIN functions

In the DATA step, the most straightforward way to truncate a value x to the interval [a, b] is to use nested calls to the MIN and MAX functions. First, take the minimum between the value and the upper bound, b. This value is less than or equal to b. Then take the maximum of the lower bound, a, and the previous result. The result is greater than or equal to a.

You can define a simple macro to encapsulate this logic:

/* clip a value x into the interval [low, high] */
%macro CLIP(low, value, high);
   (max(&low, min(&value, &high)))
%mend;
data ClipVals;
input x @@;
clipped = %CLIP(-1, x, 1);
datalines;
-2 -1.5 -0.5 0 0.5 1.5 2
;
 
proc print; run;

In this program, the value for x is clipped into the interval [-1, 1]. All values less than -1 become -1; all values greater than 1 become 1.

Clipping in PROC IML: Elementwise operators

If you are working with vectors or matrices in PROC IML, the previous code is probably not what you want to use. The MIN and MAX functions return a scalar value, which is the minimum or maximum element of the vector or matrix argument. Instead, you probably want to truncate every element of a vector by using the elementwise minimum and maximum operators. If x and y are conformal quantities (for example, vectors of the same size or a scalar and a vector), then you can use the elementwise operators, as follows:

• The elementwise maximum is x <> y.
• The elementwise minimum is x >< y.

You can encapsulate this logic by defining a SAS IML function that clips elements of a vector, x, into an interval [a,b]. The interval is represented as a two-element vector.

/* use the elementwise minimum (><) and elementwise maximum (<>) operators to 
   truncate elements of a vector into an interval. */
proc iml;
start Clip(x, ab);
   a = ab[1]; b = ab[2];
   return( a <> (x >< b) );
finish;
 
x = {-2, -1.5, -0.5, 0, 0.5, 1.5, 2};
clipped = Clip(x, {-1, 1});   /* truncate elements of x into [-1, 1] */
print x clipped;

The result is the same as for the DATA step example. The elementwise operators prevent you from having to loop over elements of the vector.

What happens with missing values?

A missing value is not a mathematical number, so the clipping operation is not defined for missing values. SAS obeys certain conventions when a missing value is used in a MIN or MAX function, or is part of an elementwise operation in IML. As written, the DATA step function will return the lower limit if the value is missing. The IML function will return the upper limit. With a little additional logic, you can determine what you want to happen by using the IFN function in the DATA step and the CHOOSE function in IML. For example, if you decide that you want the clipping function to output a missing value when the input is missing, you can rewrite the body of the macro and the function as follows:

• DATA step: ifn(missing(&value), ., max(&low, min(&value, &high)))
• PROC IML: choose(x=., ., a (x >< b))

Summary

This article shows how to clip (or truncate) values into an interval. This technique is useful for avoiding domain errors in arguments to functions that have a restricted domain. Even for functions that are defined everywhere, such as the EXP function, you can use the clipping trick to avoid numerical underflow and overflow in your programs.

The post Clip values into an interval appeared first on The DO Loop.