A Monte Carlo (MC) simulation provides both a statistical estimate for a parameter and a confidence interval for the parameter. (The confidence interval assumes that the estimate follows a normal distribution.) If you run two independent MC simulations to estimate the same parameter, how close will the two estimates be? This is illustrated by the diagram to the right. There are two MC estimates (q₁ and q₂) and two confidence intervals for an unknown parameter, μ. How close are the estimates to each other? Intuitively, you might expect to bound the distance between the estimates by the width of the larger confidence interval, but is there a better way to compare the distance between two estimates?

Yes, there is a better way. You can compute the variance of the difference between the two estimates and construct a 95% confidence interval for that difference. This article shows how to perform the computation in SAS by using the SAS IML language. We will also verify the empirical 95% coverage for an example.

Generate independent Monte Carlo estimates

To make the ideas concrete, let's look at a specific problem. Suppose you want to estimate the expected value of the Exponential distribution with scale parameter 10. Call that quantity μ. You can prove analytically that μ = 10 is the true expected value.

The following SAS IML function performs a Monte Carlo simulation to estimate μ. The function generates B independent samples of size N from the Expo(10) distribution. It calculates the mean of each sample. The Monte Carlo estimate is the average of the estimates. The Monte Carlo standard error is the standard deviation of the estimates divided by sqrt(B).

proc iml;
/* use Monte Carlo simulation to estimate the mean of the
   Expo(lambda) distribution. The simulation uses B independent
   samples of size N. Each sample is independently drawn from Expo(lambda) */
start MC_Est_Expo(N, B, lambda);
    x = j(N, B, .);
    call randgen(x, "Expo", lambda);
    sample_means = mean(x);                /* mean of each column */
    MC_est = mean(sample_means`);          /* MC estimate of statistic */
    SE_est = std(sample_means`) / sqrt(B); /* MC standard error */
    return( MC_est || SE_est );
finish;
 
N = 100;
B = 5000;
lambda = 10;
call randseed(79);
est1 = MC_Est_Expo(N, B, Lambda);
est2 = MC_Est_Expo(N, B, Lambda);
print est1[c={'MCest1' 'SE1'}], est2[c={'MCest2' 'SE2'}];

Because the algorithm relies on random number simulation, every call results in a slightly different point estimate and standard error. Here, the estimate from the first simulation is larger than the true parameter (10); the other estimate is smaller. The standard errors are very similar: approximately 0.142.

The variance of the difference

Let q₁ and q₂ be two independent MC estimates, with standard errors SE₁ and SE₂. Since both estimates are computed by the same method, the expected value of their difference is zero: E[q₁ - q₂] = 0. (If we used different methods, and those methods are biased by different amounts, the expected value will not be zero.)

How much variance should we expect for the difference? Because q₁ and q₂ are independent, the variance of their difference is the sum of the individual variances:
Var(q₁ - q₂) = SE₁² + SE₂²

Let SE_diff be the square root of this sum. For this example, both estimates are normally distributed, so the difference also follows a normal distribution: q₁ - q₂ ∼ N(0, SE_diff). You can use this distribution to construct a 95% interval for the difference. Let's do that, then check whether the distance between these estimates |q₁ - q₂| is less than the width of the interval.

/* given (q1, SE1) and (q2, SE2), what is a confidence interval for the 
   difference q1-q2? The expected value E[q1-q2] = 0, but what is the variance? */
q1 = est1[1];  q2 = est2[1]; 
SE1 = est1[2]; SE2 = est2[2];
/* Standard deviation of the DIFFERENCE between two independent estimates */
SE_diff = sqrt(SE1**2 + SE2**2);
print SE1 SE2 SE_diff;
 
/* Can we bound |q1 - q2| with high probability? 
   Use a prediction interval based on the fact that 
   q1-q2 ~ N(0, SE_diff) */
w95 = quantile("Normal", 0.975) * SE_diff;
if abs(q1 - q2) < w95 then 
    print "The difference |q1 - q2| is inside a 95% confidence interval";
else
    print "The difference |q1 - q2| is NOT inside a 95% confidence interval";

For these estimates, the difference is, indeed, less than the width of a 95% confidence interval. Of course, we might have gotten lucky. The next section investigates what happens if we repeat this experiment many times.

Verify the empirical coverage

According to probability theory, the difference between the two estimates should fall inside the calculated interval 95% of the time. However, I always like to verify statistical theory by running a simulation in SAS.

The following program repeats the experiment 1,000 times. For each iteration, the program generates two independent Monte Carlo estimates, computes the 95% confidence interval for their difference, and saves a 0/1 binary value that indicates if the difference falls inside the CI.

/* If we repeat this computation many times, about 95% of the differences
   should be in the 95% CI */
C = 1000;
crit_val = quantile("Normal", 0.975);  /* ~ 1.96 */
inCI95 = j(C, 1, 0);
do i=1 to C;
    est1 = MC_Est_Expo(N, B, Lambda);
    est2 = MC_Est_Expo(N, B, Lambda);
    q1 = est1[1];  q2 = est2[1]; 
    SE1 = est1[2]; SE2 = est2[2];
    SE_diff = sqrt( (SE1**2 + SE2**2) );
    w95 = crit_val * SE_diff;
    inCI95[i] = (abs(q1 - q2) < w95); 
end; 
 
propInCI = mean(inCI95);
print propInCI[c={'Empirical 95% Cov'}];

Ah, success! The empirical coverage for this random number seed is close to the theoretical coverage of 0.95. This confirms that adding the individual variances provides a statistical interval for comparing the distance between two Monte Carlo estimates of the same quantity.

Modifying this technique

You can change the critical value if you want to test the difference with more confidence. For example, you can use three standard errors to obtain a 99.7% confidence interval, which is known as the "three-sigma rule."

You can also compare the distance between two estimates that come from different (unbiased) estimates. For example, if one estimate is computed by using B₁=5,000 random samples and the other is computed by using B₂=3,000 random samples, the standard error for the second estimate will be larger, but you can still compute the standard error of the difference.

If you are testing a routine that produces a random or quasi-random result, you can use this technique to compare the routine's result to a MC validation procedure that you wrote yourself. Suppose the routine returns the value q₁. Your validation procedure gives the value q₂. 95% of the time, q₁ will be close to q₂. But there is a 5% chance that they are not sufficiently close. Does that indicate that the estimates are wrong? No, most likely you got unlucky. One option is to rerun the routine and your verification procedure again. This gives new estimates q₃ and q₄. These new results also have a 5% chance of failure, but the probability that BOTH tests fail sequentially is a very small 0.25%. If that probability is not small enough for you, you can use 99.7% CIs instead of 95% CIs.

Summary

This article shows a technique for comparing the difference between two Monte Carlo estimates of the same quantity. You can compute the variance of the difference as the sum of the two individual variances. In many cases, the expected value of the difference is zero, and the difference is normally distributed. (You should check this assumption.) Therefore, you can derive a formula for the magnitude of the difference. You can use this formula to verify the consistency of a single Monte Carlo simulation, or you can use it to verify that two different MC algorithms provide consistent estimates of the same quantity.

The post How to compare two independent Monte Carlo estimates appeared first on The DO Loop.

A SAS programmer recently asked me questions about a zero-inflated (ZI) probability distribution. Some readers might be familiar with distributions that are related to zero-inflated counts. These distributions are used in some linear regression models. A zero-inflated model is a mixture of two distributions: a discrete probability mass situated at X=0, and a second probability distribution, which describes the probability that a random variable is not zero. In regression, two common zero-inflated distributions are the ZI Poisson (ZIP) and the ZI negative binomial (ZINB) distributions.

However, the SAS programmer wanted to implement a mixture distribution that is the mixture of a discrete point-mass distribution and a continuous distribution. The discrete distribution gives the probability that X=0. The continuous distribution gives the probability that X > 0.

This article shows how to use SAS to create functions for the zero-inflated Beta distribution. The image to the right shows the PDF for a ZI Beta distribution. The techniques in this article can be generalized to other mixture distributions, such as a ZI gamma distribution. This article implements the functions in PROC FCMP and calls them from the DATA step.

What is the zero-inflated Beta distribution?

The Beta distribution is a two-parameter distribution on the interval (0,1). The two positive shape parameters are often called "alpha" and "beta." To reduce confusion, I capitalize the name of the distribution and lowercase the name of the parameter. Depending on the relative magnitudes of the shape parameters, the probability density function can be U-shaped, J-shaped, L-shaped, or mound-shaped (sometimes called "inverted-U shaped").

The equation for the probability density function (PDF) of a random variable that follows a zero-inflated distribution is
     f(x) = π δ(x) + (1-π) f₂(x; α, β) I(x > 0)
where δ(x) is the Dirac delta function, π is the mixing probability, f₂ is the probability density function of the second distribution, and I is the indicator function (also called the Heaviside function). For this article, f₂ is the PDF of a Beta distribution.

This article uses the values α = 2 and β = 3 for the shape parameters. For the probability of mixing, this article uses π = 0.2. In programs, I use the variable p0 instead of π for the mixing probability so that the parameter will not be confused with the mathematical constant 3.14159....

Define the zero-inflated functions in PROC FCMP

This article defines the four main functions that statistical programmers need: The PDF, CDF, quantile, and random functions for the ZI Beta distribution. They are defined in a PROC FCMP library and can be called from the SAS DATA step. To simplify the article, I've put the PROC FCMP code in the Appendix. To reproduce the images in this article, you must run the PROC FCMP code before running the examples! I have stored the definitions in the WORK libref, but you replace WORK with a libref to a permanent location if you want to persist the definitions. If you do, be sure to specify the CMPLIB= option in each new session of SAS. The CMPLIB= option tells SAS where you stored the functions.

Parameter values for the examples

To keep this article short, all examples use the same parameter values for the Beta distribution. The %DefineParms macro defines three parameter values. It sets the mixing probability (p0) to 0.2. It sets the Beta parameters to α = 2 and β = 3. You can modify the macro to explore other ZI Beta shapes.

%macro DefineParms;
   p0    = 0.2;
   alpha = 2;
   beta  = 3;
   length Params $40;
   Params = cats("alpha=",put(alpha,BEST3.), ", beta=", put(beta,BEST3.));
%mend;

The ZI Beta PDF

When I study a new distribution, I start by visualizing the probability density function (PDF). The following DATA step calls the PDF_ZIBeta function, then uses PROC SGPLOT to visualize the PDF curves. Remember to run the code in the Appendix first!

data PDF_Test;
%DefineParms;
dx = 0.01;
do x = 1E-6, dx/2, dx to 1 by dx;
   PDF = PDF_ZIBeta(x, p0, alpha, beta);
   output;
end;
t = 0; output;
keep p0 x PDF Params t;
run;
 
title "PDF of Zero-Inflated Beta";
proc sgplot data=PDF_Test;
   series x=x y=PDF / group=Params lineattrs=(thickness=2);
   dropline x=t y=p0 / dropto=x lineattrs=GraphData1(thickness=5);
   yaxis grid min=0 offsetmin=0;
   xaxis grid;
run;

The graph is shown at the top of this article. I used a DROPLINE statement to display a vertical bar, which indicates the discrete probability mass at x=0. The Beta curve defines the probability density for x > 0.

The ZI Beta CDF

The cumulative distribution function (CDF) is the integral of the PDF. Given a value, x, the CDF at x tells you the probability that a random observation drawn from the ZI Beta distribution will be less than or equal to x. Because of the point-mass at 0, the CDF for the ZI Beta distribution has the value p0 (here, 0.2) at x=0 and increases monotonically as x increases. The following DATA step calls the CDF_ZIBeta function, then uses PROC SGPLOT to visualize the CDF curves:

data CDF_Test;
%DefineParms;
dx = 0.01;
do x = 1E-6, dx/2, dx to 1 by dx;
   CDF = CDF_ZIBeta(x, p0, alpha, beta);
   output;
end;
t = 0; output;
keep p0 x CDF Params t;
run;
 
title "CDF of Zero-Inflated Beta";
proc sgplot data=CDF_Test;
   series x=x y=CDF / group=Params lineattrs=(thickness=2);
   dropline x=t y=p0 / dropto=x lineattrs=GraphData1(thickness=2);
   yaxis grid min=0 offsetmin=0;
   xaxis grid;
run;

For any value of x along the horizontal axis, the graph shows the probability (on the vertical axis) that a random observation is less than or equal to x. For example, the function at x=0.4 shows that there is approximately a 63% chance that a random variate from that distribution will have a value less than or equal to 0.4. Notice that there is a 20% probability that a random observation is 0.

The ZI Beta quantile function

The quantile function is the inverse of the CDF function. Given a probability, p in (0,1), the quantile function at p is the x value such that a random observation drawn from the ZI Beta distribution will be less than or equal to x with probability p. For a ZI distribution, you can use the quantile function of the continuous distribution, but you need to rescale the probability to account for the point-mass at 0. The following DATA step calls the quantile_ZIBeta function, then uses PROC SGPLOT to visualize the quantile curves.

data Quantile_Test;
%DefineParms;
dp = 0.01;
do p = 0 to 1-dp by dp, 0.999, 0.9999, 1;
   x = quantile_ZIBeta(p, p0, alpha, beta);
   output;
end;
keep p0 p x Params;
run;
 
title "Quantile Function of Zero-Inflated Beta";
proc sgplot data=Quantile_Test;
   series x=p y=x / group=Params lineattrs=(thickness=2);
   yaxis grid min=0 offsetmin=0;
   xaxis grid;
run;

Notice that the graph of the quantile function is the "flipped image" of the CDF graph. Notice also that the derivative of the quantile becomes infinite as x → 1. This is because the slope of the CDF function approaches 0 as x → 1.

Random variates from the ZI Beta distribution

For any mixture distribution, generating random variates requires two steps. The first step is to generate a binary random variate. (In SAS, this is done by using the "Bernoulli" distribution.) With probability π, you return a random variate from the first distribution; with probability 1-π, you return a random variate from the second distribution. For the ZI Beta distribution, this means that you return 0 with probability π, and a random variate from the Beta distribution with probability 1-π. If you draw a histogram of the resulting random variates, the bin that includes 0 will display the counts of the random variates near 0 from both distributions.

The following DATA step calls the rand_ZIBeta function 1000 times. Approximately 20% of the random variates are from the point-mass at 0. The remaining variates are from the Beta distribution.

data Rand_Test;
%DefineParms;
call streaminit(1234);
N = 1000;
do j  = 1 to N;
   x = rand_ZIBeta(p0, alpha, beta);
   output;
end;
keep x Params;
run;
 
title "Random Variates for the Zero-Inflated Beta";
proc sgplot data=Rand_Test;
   histogram x;
   yaxis grid;
   xaxis grid;
run;

Because the program generates 1,000 random variates, the shape of the histogram is similar to the shape of the PDF curve. Specifically, note that bar that includes x=0 is slightly more than 20% for this random sample.

Define the functions for the ZI Beta distribution in SAS IML

Some programmers prefer to use the SAS IML language, especially when conducting a simulation study. Appendix 2 presents SAS IML module definitions. You can store the definition, then call them from a SAS IML program.

Summary

This article demonstrates how to implement a zero-inflated distribution in SAS. This article implements the zero-inflated Beta distribution, but you can use the same techniques for other zero-inflated distributions. Appendix A of this article defines four related functions in PROC FCMP in Base SAS: the PDF, CDF, quantile, and random variate functions. SAS IML versions of the functions are defined in Appendix B.

Appendix A: FCMP functions for the zero-inflated Beta distribution

This appendix shows how to use PROC FCMP in SAS to add new DATA step functions. It follows the same technique as a previous article about the generalized extreme-value (GEV) distribution. Run this PROC FCMP step before trying to call the functions in the DATA step.

/* Define the PDF, CDF, QUANTILE, and RAND functions for the 
   ZI Beta distribution in PROC FCMP.  
 
   The support for the ZI Beta distribution is X in [0,1).
   alpha : Shape1 > 0
   beta  : Shape2 > 0
 
   For simplicity, the function assume alpha > 0 and beta > 0.
   The PDF is undefined if x=1 and beta < 1. It is well-defined
   for other parameter values. For simplicity, we always 
   return a missing value with x=1. Thus, the variates for
   the distribution are defined on [0,1). */
proc fcmp outlib=work.funcs.ProbDist;
   function PDF_ZIBeta(x, p0, alpha, beta);
      if x < 0 | x >= 1 then return (.);
      f = ifn(x=0, p0, (1-p0)*pdf('Beta', x, alpha, beta));
      return (f);
   endsub;
 
   function CDF_ZIBeta(x, p0, alpha, beta);
      if x < 0 | x >= 1 then return (.);
      F = ifn(x=0, p0, p0+(1-p0)*cdf('Beta', x, alpha, beta));
      return (F);
   endsub;
 
   function quantile_ZIBeta(p, p0, alpha, beta);
      if p < 0 | p > 1 then return (.);
      if p <= p0 then return (0);
      if p = 1 then return (1);
      z = (p-p0)/(1-p0);  /* rescale the probability on [p0,1] */
      x = quantile('Beta', z, alpha, beta);
      return (x);
   endsub;
 
   function rand_ZIBeta(p0, alpha, beta);
      b = rand("Bernoulli", 1-p0);
      if b=0 then return (0);
      x = rand('Beta', alpha, beta);
      return (x);
   endsub;
quit;
options cmplib=work.funcs;  /* define location of functions so the DATA step can find them */

Appendix B: SAS IML functions for the zero-inflated Beta distribution

This appendix shows how to define the PDF, CDF, QUANTILE, and RAND functions for the ZI Beta distribution in the SAS IML language.

proc iml;
/* PDF of Zero Inflated Beta Distribution
   x can be a vector; p0, alpha, and beta are scalars. */
start PDF_ZIBeta(x, p0, alpha, beta);
   f = j(nrow(x), ncol(x), .);
   idx = loc(x = 0);
   if ncol(idx)>0 then 
      f[idx] = p0;
   idx = loc(0 < x & x < 1);
   if ncol(idx)>0 then 
      f[idx] = (1-p0)*PDF('beta', x[idx], alpha, beta);
   return f;
finish;
 
/* CDF of Zero Inflated Beta Distribution
   x can be a vector; p0, alpha, and beta are scalars. */
start CDF_ZIBeta(x, p0, alpha, beta);
   F = j(nrow(x), ncol(x), .);
   idx = loc(x = 0);
   if ncol(idx)>0 then 
      F[idx] = p0;
   idx = loc(x = 1);
   if ncol(idx)>0 then 
      F[idx] = 1;
   idx = loc(0 < x & x < 1);
   if ncol(idx)>0 then 
      F[idx] = p0 + (1-p0)*CDF('beta', x[idx], alpha, beta);
   return F;
finish;
 
/* Quantile function of Zero Inflated Beta Distribution
   p can be a vector in (0,1); p0, alpha, and beta are scalars. */
start Quantile_ZIBeta(p, p0, alpha, beta);
   x = j(nrow(p), ncol(p), .);
   idx = loc(0 <= p & p <= p0);
   if ncol(idx)>0 then 
      x[idx] = 0;
   idx = loc(p = 1);
   if ncol(idx)>0 then 
      x[idx] = 1;
   idx = loc(p0 < p & p < 1);
   if ncol(idx)>0 then do;
      z = (p[idx]-p0)/(1-p0);  /* rescale the probability on [p0,1] */
      x[idx] = quantile('Beta', z, alpha, beta);
   end;
   return x;
finish;
 
/* Random variate function of Zero Inflated Beta Distribution
   N is the number or random variates; p0, alpha, and beta are scalars. */
start Rand_ZIBeta(N, p0, alpha, beta);
   x = j(N, 1, .);
   b = randfun(N, 'Bernoulli', 1-p0);
   idx = loc(b = 0);
   if ncol(idx)>0 then 
      x[idx] = 0;
   idx = loc(b = 1);
   nBeta = ncol(idx);
   if nBeta>0 then
      x[idx] = randfun(nBeta, 'Beta', alpha, beta);
   return (x);
finish;
store module=(PDF_ZIBeta CDF_ZIBeta quantile_ZIBeta Rand_ZIBeta);
QUIT;
 
/* Test the functions for the ZI Beta distribution */
proc iml;
load module=(PDF_ZIBeta CDF_ZIBeta quantile_ZIBeta Rand_ZIBeta);
 
p0    = 0.2;
alpha = 2;
beta  = 3;
 
dx = 0.01;
x = 0 // 1E-6 // T(do(dx, 1, dx));
PDF = PDF_ZIBeta(x, p0, alpha, beta);
 
title "PDF of Zero-Inflated Beta";
call series(x, PDF) grid={x y};
 
title "CDF of Zero-Inflated Beta";
CDF = CDF_ZIBeta(x, p0, alpha, beta);
call series(x, CDF) grid={x y} other="refline 0.2/axis=y; yaxis grid values=(0 to 1 by 0.1);";
 
title "Quantile Function of Zero-Inflated Beta";
dp = 0.01;
p = T(do(0, p0, dp)) // T(do(p0+dp/10, 1, dp));
q = quantile_ZIBeta(p, p0, alpha, beta);
call series(p, q) grid={x y};
 
title "Random Variates for the Zero-Inflated Beta";
x = Rand_ZIBeta(1000, p0, alpha, beta);
call histogram(x) grid={y};
QUIT;

The post Implement a zero-inflated distribution in SAS appeared first on The DO Loop.

Certain fields, such as quantitative finance and climatology, routinely use models with thicker tails. If you need to analyze data that have "heavy tails," one possible choice is the multivariate t (MVT) distribution. The MVT distribution has a parameter, ν, that controls the tail thickness. For historical reasons, ν is called the degrees-of-freedom parameter. When used for inferential statistics, the parameter is an integer. However, for the MVT distribution, ν is a shape parameter that is not restricted to integer values. For small values of ν, the MVT provides a heavy-tailed alternative to the MVN distribution. As ν → ∞, the MVT distribution approaches the multivariate normal distribution.

I have previously written about how to visualize the multivariate t distribution in SAS, including how to evaluate the multivariate t PDF and how to generate random variates from MVT by using the RANDMVT function in SAS IML software. This article discusses how to compute the cumulative distribution function (CDF) for the multivariate t distribution. The CDF requires computing a high-dimensional integral. The article includes a SAS IML function that evaluates the CDF for the multivariate t distribution in two dimensions.

A general formula for the CDF of the multivariate t distribution

This section presents the advanced math that enables you to compute the CDF for the multivariate t (MVT) distribution. If you prefer programming over math, you can skip to the next sections.

The book Computation of Multivariate Normal and t Probabilities (Genz and Bretz, 2009) shows that you can obtain the CDF of the standard MVT(Σ, ν) distribution as an integral of the CDF of the MVN(0, Σ) distribution. The goal is to evaluate the probability that a k-dimensional random variable X ~ MVT(Σ, ν) is in a hyper-rectangular region: P(a1 < X1 < b1, a2 < X2 < b2, ..., ak < Xk < bk). Here Σ is a k x k positive semidefinite covariance matrix that defines the correlations between the components of X.

The vectors a = (a1, a2, ..., ak) and b = (b1, b2, ..., bk) define the rectangular region. To obtain the probability that a random variate is in this region, you use the components of a for the lower limits of integration, and the components of b for the upper limits. The components of a can be -∞, and the components of b can be +∞.

Although you can define the probability directly as a multivariate integral over the hyper-rectangular region, Genz and Bretz show that you can also compute the probability by integrating a function that includes the CDF for a MVN distribution. The first formula appears on p. 3 and is listed as Eqn 1.3. The formula shows how to compute the MVT CDF by integrating a related CDF for the MVN distribution over the infinite interval (0, ∞). Then, in Chapter 4 (Eqn 4.6, p. 32-33), Genz and Bretz show that you can transform the integral on the infinite interval (0, ∞) into an equivalent integral on (0,1). The formula is:
$T_k(\mathbf{a}, \mathbf{b}; \boldsymbol{\Sigma}, \nu) = \int_{0}^{1} \Phi_k \left( \frac{\chi_\nu^{-1}(t)\mathbf{a}}{\sqrt{\nu}}, \frac{\chi_\nu^{-1}(t)\mathbf{b}}{\sqrt{\nu}}; \boldsymbol{\Sigma} \right) dt$

In this formula:

• T_k(a, b; Σ, ν) is the CDF function for the MVT distribution. It gives the probability that a random variate from MVT(Σ, ν) is in the region defined by the vectors a and b.
• Φ_k is the MVN CDF in dimension k. $\Phi_k (\mathbf{L}, \mathbf{U}; \boldsymbol{\Sigma})$ is the probability that a random variable from the centered MVN distribution is in the hyper-rectangular region defined by the lower-limit vector, L, and the upper-limit vector, U. The vector L = $\frac{\chi_\nu^{-1}(t)\mathbf{a}}{\sqrt{\nu}}$, and U = $\frac{\chi_\nu^{-1}(t)\mathbf{b}}{\sqrt{\nu}}$, where t is the dummy integration variable in (0,1).
• χ_ν^-1 is the quantile function for the chi distribution. If you've never heard of the chi distribution, don't worry. It is closely related to the familiar chi-square distribution. In fact, if Y ~ χ²(ν) is distributed as a chi-squared distribution with ν degrees of freedom, then X = sqrt(Y) is distributed as a chi distribution with ν degrees of freedom. This definition means that you can use the QUANTILE("CHISQ") distribution in Base SAS to compute the quantile function of the chi distribution. You simply take the square-root of the resulting quantile!

This formula enables you to compute the exact CDF of a MVT distribution by integrating a function that involves only the (transformed) probability of the MVN distribution over the interval (0,1). In SAS, you can use the PROBBNRM function to evaluate 2-D MVN CDF. And SAS IML software enables you to integrate arbitrary functions! So, let's write a SAS IML function that uses the formula to evaluate probabilities for the 2-D MVT distribution. The function will enable you to compute exact CDF values in this case.

The CDF of the bivariate t distribution in SAS

The formula for Eqn 4.6 provides the CDF on any rectangular region. However, for ease of presentation, I will present the simpler case of the "left-tailed" CDF in 2-D, which is the probability that a random variable is in the region (-∞, b1)x(-∞, b2). In other words, I will set the vector of lower limits to a = (-∞, -∞). This computation can be used to compute the probability of general rectangular regions.

A left-tailed region is shown in the image to the right. In the image, I generated 1,000 random variates according to a bivariate t distribution with ν=5 degrees of freedom. Of those points, 180 of them fall into the region {(x,y) | x < -1, y < 0}, which means that a Monte Carlo estimate for the probability is 0.18. Let's use the formulas in the previous section to obtain a better estimate of the probability.

You can use the QUAD routine in the SAS IML language to perform numerical integration. It is helpful to use the PEAK= option to tell the QUAD routine to avoid evaluating the integral near the boundary of the integration region. The following IML program defines and store three functions. The BVT_CDF is the top-level function that evaluates the bivariate t CDF. 'BVT_CDF' is short for 'bivariate t CDF'.

/* A useful result in computational statistics is Eqn 1.3 and Eqn 4.6 in 
   Genz and Bretz (2009, Computation of Multivariate Normal and t Probabilities, p. 3, 32-33).
   The equation shows that you can compute any multivariate t CDF by computing an integral
   that uses the multivariate normal (MVN) CDF. 
   Specifically, Eqn 4.6 shows that you can use the quantile of the CHI distribution 
   (CHI, not chi-squared!) to transform Eqn 1.3 into an integral on the interval (0,1).
*/
proc iml;
/* truncate every element of x into the interval [ab[1], ab[2]].
   See https://blogs.sas.com/content/iml/2026/02/04/clip-values.html
*/
start TruncateToRange(x, ab);
   a = ab[1]; b = ab[2];
   return( (a <> x) >< b );
finish;
 
/* Define the integrand from Equation (4.6) in Genz (2009)
   by calling the PROBBNRM function in Base SAS to compute the bivariate normal CDF.
   For numerical reasons, limit the evaluation of PROBBNRM to [-10,10] in each coordinate. 
   In case nu is large, use the exp(log(.)) trick to evaluate the constants.
*/
start BVT_CDF_Integrand(t) global(g_b, g_rho, g_nu);
   b = g_b; rho = g_rho; nu = g_nu;
   if      t<=0 then s=0;
   else if t>=1 then s=1E6;
   else s = sqrt( quantile("ChiSq", t, nu));
   w = s * b / sqrt(nu);
   w = TruncateToRange(w, {-10, 10}); /* protect against floating point underflow; max value is 10 */
   f = probbnrm(w[1], w[2], rho);
   return( f );
finish;
 
/* wrapper function to set global variables and call the integration */
start BVT_CDF(b, rho, nu) global(g_b, g_rho, g_nu);
   g_b = b; g_rho = rho; g_nu = nu;
   call quad(BVT_prob, "BVT_CDF_Integrand", {0 1}) peak=0.1;
   return( BVT_prob );
finish;
store module=(TruncateToRange BVT_CDF_Integrand BVT_CDF);
QUIT;

Having defined the functions, you can now call them to compute the CDF for the bivariate t distribution, as follows:

proc iml;
load module=_all_;
 
/* call the BVT_CDF with different values of b and nu */
b = {-1, 0};     /* upper limits for X1 and X2 */
rho = sqrt(2)/2; /* correlation */
nu = 5;          /* DOF */
/* first, compute the probability for the MVN distribution */
BVN_prob = probbnrm(b[1], b[2], rho);
print (b`)[c={'b1' 'b2'}] rho BVN_prob[L="Bivariate Normal CDF"];   
 
/* compare to the MVT probability (nu=5) */
BVT_prob = BVT_CDF(b, rho, nu);
print nu (b`)[c={'b1' 'b2'}] rho BVT_prob[L="Bivariate t CDF"];

For these calculations, we assume that two random variables have the correlation 0.7071. The first output uses the built-in PROBBNRM function in SAS to compute the probability P(Z1 < -1, Z2 < 0) for Z ~ MVN(0, R), where R is the correlation matrix with off-diagonal elements ρ. The output shows the bivariate normal CDF at (-1,0) is 0.146. In contrast, the CDF for the bivariate t distribution with 5 DoF is 0.162 because a t distribution has heavier tails. This is the probability P(T1 < -1, T2 < 0) for T ~ MVT(0; R, ν). This value is much more accurate than the Monte Carlo estimate shown in the previous section. In addition, the Monte Carlo estimate suffers from variability: If you choose a different set of 1,000 random points, you will get a different Monte Carlo estimate. In contrast, the integration method is deterministic and highly accurate.

If you reverse the correlation of the variables, you get two negatively correlated variables. You can compute the CDF probability for the MVT distribution with ρ = -0.7071 by using the following call:

rho = -sqrt(2)/2; /* correlation */
BVT_prob2 = BVT_CDF(b, rho, nu);
print nu (b`)[c={'b1' 'b2'}] rho BVT_prob2[label="Bivariate t CDF"];

The probability in the region is about 0.019, which is smaller than for the positively correlated variables. You can understand this smaller number by looking back at the scatter plot in the previous section. That scatter plot has a southwest-to-northeast orientation. When you reverse the correlation, the analogous scatter plot will have a southeast-to-northwest orientation. (Mentally flip the scatter plot about the vertical line x=0.) The new scatter plot will have only a small proportion of points (about 19 per thousand) that fall into the specified region.

Summary

This article shows a direct method to evaluate the CDF for the bivariate t distribution in SAS. The computation relies on a mathematical result in Genz and Bretz (2009), which shows that you can evaluate the probability that a MVT variable (in any dimension) is in a rectangular region by integrating the CDF for the MVN distribution over the interval (0, 1). To simplify the presentation, I showed the special case of computing the CDF for the bivariate t distribution on left-tailed regions {(T1,T2) | T1 < b1, T2 < b2}, where T = (T1,T2) is a random variable distributed as MVT(Σ, ν). You can use this left-tail probability to compute the probability in any other 2-D rectangular region.

It is difficult to compute the exact CDF of the MVN distribution in k dimensions when k > 2. However, you can use the formula and numerical methods such as quasi-Monte Carlo computations to compute the CDF of the MVT distribution to high precision.

The post How to evaluate the multivariate t CDF in SAS appeared first on The DO Loop.

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

A simple transformation that changes correlation

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

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

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

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

A SAS IML program that changes correlation

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

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

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

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

Summary

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

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

Correlation matrices that have a common spectrum

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

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

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

The Marsaglia and Olkin method for random correlation matrices

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

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

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

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

Generate random correlation matrices in SAS

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

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

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

Summary

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

Appendix: The Marsaglia and Olkin Algorithm in IML

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

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

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

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

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

The computation of MVN probabilities

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

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

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

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

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

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

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

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

An overview of the Genz transformation

Genz's method consists of three main steps:

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

A 2-D Implementation in SAS IML

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

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

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

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

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

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

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

The Genz integrand for uncorrelated variables

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

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

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

The Genz integrand for small probabilities

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

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

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

Summary

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

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

Appendix: The traditional Monte Carlo estimate

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

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

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

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

The Gram-Schmidt orthogonalization of a matrix

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

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

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

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

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

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

Computing orthogonal vectors to a linear subspace: The manual way

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

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

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

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

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

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

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

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

Computing orthogonal vectors to a linear subspace: The easier way

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

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

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

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

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

Summary

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

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

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

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

An unconstrained feature set for Scotch whiskies

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

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

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

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

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

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

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

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

A sparseness-constrained NMF algorithm

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

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

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

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

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

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

The tradeoffs between sparseness and goodness of fit

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

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

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

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

Summary

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

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

Appendix: Sparseness and goodness of fit

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

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

Sparseness versus sparsity

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

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

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

Hoyer's sparseness measure

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

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

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

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

The range of the Hoyer measure

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

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

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

Examples of Hoyer's sparseness

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

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

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

Visualize Hoyer's sparseness for 2-D vectors

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

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

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

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

Visualize Hoyer's sparseness for 3-D vectors

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

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

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

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

A SAS function to compute Hoyer's sparseness

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

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

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

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

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

Summary

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

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

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

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

The Log-KDE method

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

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

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

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

Implement the method in SAS

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

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

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

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

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

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

Visualize the back-transformed KDE

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

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

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

Peculiarities of the log-KDE method

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

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

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

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

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

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

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

Summary

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

Further Reading

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

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