A SAS customer recently asked me about how to implement the modified Bessel function of the second kind in SAS. This Bessel function, which is often denoted by K_ν(x), appears in the formulas for the generalized inverse Gaussian (GIG) distribution. The modified Bessel function of the second kind is also used in survival analysis, spatial statistics, financial modeling, and to model variance parameters in Bayesian models. Although SAS does not have a built-in "KBESSEL" function, you can define the K_ν(x) function in SAS IML.

The definition of the modified Bessel function of the second kind

For the purpose of this article, the modified Bessel function of the second kind, K_ν(x), where the parameter ν is a real number and the domain of the function is x > 0. A good way to compute the function is to use the integral definition from Equation 10.32.10 of the NIST Digital Library of Mathematical Functions:
$K_\nu(x) = \frac{1}{2} \left( \frac{x}{2} \right)^\nu \int_0^\infty t^{-\nu-1} \exp\left( -t - \frac{x^2}{4t} \right) \, dt$

In this formula, x and ν are parameters to the integrand. The dummy variable, t, is integrated from 0 to infinity. The integrand has a peak that can be computed analytically by using calculus. The peak occurs at
$t_{peak} = \frac{\sqrt{(\nu + 1)^2 + x^2} - (\nu + 1)}{2}$

In the SAS IML Language, the QUAD subroutine is used to compute integrals. Knowing the location of the peak is helpful because the QUAD subroutine in IML contains a PEAK= option that you can specify to assist the numerical integration routine. (Alternatively, you can perform a change of variables u = t / t_peak so that the peak is always at u=1. I leave that as an exercise.)

The optional exponential scaling factor

The documentation for the the IBESSEL function, which is the modified Bessel function of the first kind, contains an optional binary argument. By default, the call IBESSEL(nu, x) computes the function I_ν(x). However, the function decays rapidly to zero as x increases, which makes the numerical computation difficult when x is moderately large (like x ≥ 50). Accordingly, a numerical trick is to compute the expression exp(x)*I_ν(x), which you can do by using the syntax IBESSEL(nu, x, 1).

A similar situation occurs for the Bessel function of the second kind, K_ν(x), which also decays quickly to zero for moderately large x. (Equation 10.40.2 in the NIST Digital Library shows that the function decays asymptotically like exp(-x)/sqrt(x).) When x is moderately large, it becomes numerically challenging to evaluate the function by using a standard numerical quadrature routine like CALL QUAD. The integration routine can prematurely terminate because the value of the integral is less than the algorithm's absolute error tolerance.

Accordingly, it is common in numerical software for the syntax to support an optional binary flag that causes the function to return a "scaled value," which is exp(x)*K_ν(x). The exponential scaling factor is evaluated as part of the integrand. This keeps the peak of the integrand near 1 and permits the quadrature routine to evaluate the integral accurately.

Not only does this trick enable you to compute the scaled value exp(x)*K_ν(x) to high precision, but it also enables you to find the value of log(K_ν(x)) accurately. Taking the natural log of the scaled value gives log(exp(x)*K_ν(x)) = x + log(K_ν(x)), which you can solve for log(K_ν(x)). This computation is shown in the next section. The log of the Bessel function is used in applications such as loglikelihood computations.

Implement modified Bessel function of the second kind in SAS

Because SAS already supports the JBESSEL and IBESSEL functions, I will reserve the name KBESSEL for a future built-in SAS implementation. Instead, I will define a SAS IML function named BesselK, which takes three parameters: the parameter ν, the variable x, and an optional binary flag. If the binary flag is 0 (the default), the function will return K_ν(x), otherwise the function returns the scaled value exp(x)*K_ν(x).

proc iml;
/* Define the integrand function for the unscaled modified 
   Bessel function of the second kind, K_{nu}(x)
 */
start BesselK_integ(t) global(g_BessK_nu, g_BessK_x);
   nu = g_BessK_nu; x = g_BessK_x;
   return( t**(-nu - 1) * exp(   -t - (x**2)/(4*t)) );
finish;
 
/* Define the integrand function that includes the EXP term:
   exp(x)*K_{nu}(x). Use this scaled function when x is moderately large or you want log(K_{nu}(x)).
*/
start BesselK_Scale_integ(t) global(g_BessK_nu, g_BessK_x);
   nu = g_BessK_nu; x = g_BessK_x;
   return( t**(-nu - 1) * exp(x -t - (x**2)/(4*t)) );
finish;
 
/* Define the wrapper function */
start BesselK(nu, x, scale=0) global(g_BessK_nu, g_BessK_x);
   g_BessK_nu = nu;
   if scale = 0 then 
      funcName = "BesselK_integ";
   else 
      funcName = "BesselK_Scale_integ";
   n = nrow(x);
   K = j(n, 1, .); /* the result vector K_nu(x). Use missing for domain error such as x <= 0 */
   tau = nu+1;
   do i = 1 to n;
      if x[i] > 0 then do;
         g_BessK_x = x[i];
         /* Calculate the exact peak of the integrand */
         pk = (sqrt(tau**2 + x[i]**2) - tau) / 2;
         call QUAD(integral_at_x, funcName, {0 .I}) PEAK=pk;
         K[i] = integral_at_x;
      end;
   end;
   /* Multiply the integral by the leading constants */
   K = 0.5 * (x / 2)##nu # K;
   return(K);
finish;
store module=(BesselK_integ BesselK_Scale_integ BesselK);
QUIT;

Let's test the function. First, let's visualize the function K_ν(x) for ν=1 and x in (0, 3].

proc iml;
load module=_all_;
 
x = T(do(0,3,0.05));
nu = 1;
 
title "Modified Bessel Function of the Second Kind (nu=1)";
K = BesselK(nu, x);
call series(x, K) grid={x y};

The graph shows that the K_ν(x) decays to 0 rapidly. Let's turn on the binary flag and visualize the function exp(x)*K_ν(x) on the same domain:

title "Scaled Modified Bessel Function exp(x)*K(x) (nu=1)";
expx_K = BesselK(nu, x, 1);
call series(x, expx_K) grid={x y};

The graph shows that the exponentially scaled function still decays, but not as rapidly. For example, whereas K₁(30) ~ 2.5E-14, the scaled function at x=30 is much larger: exp(30)*K₁(30) ~ 0.23.

By using the scaled version, you can compute log(K_ν(x)), which is useful for various applications:

/* What good is the exp scaling? 
   1. Computations are more accurate for large values of x. 
   2. You can easily compute log(K_{nu}(x)) */
x = 50;
expx_K = BesselK(nu, x, 1);
logK = log(expx_K) - x;
print x expx_K logK;

If you plan to compute log(K_ν(x)) regularly, you could encapsulate this logic into a self-contained module.

Summary

This article shows how to implement the modified Bessel function of the second kind, K_ν(x), in SAS. The method involves integrating a complicated expression over an infinite domain. In addition, the article shows how you can compute the exponentially scaled expression exp(x)*K_ν(x), which enables you to evaluate log(K_ν(x)).

In a subsequent article, I will show how you can use the K_ν Bessel function to evaluate the PDF and log-PDF of the generalized inverse Gaussian distribution in SAS.

The post The modified Bessel function of the second kind in SAS appeared first on The DO Loop.

This article provides the mathematical background for a discrete convolution operation. It discusses a SAS IML function, DFCONV, that you can use to compute a discrete convolution. The convolution operation is illustrated by using two classical statistical applications. The first is the sum of two independent discrete random variables. The second is the computation of a moving average for a time series.

The formula for a discrete convolution

The convolution operation comes in two forms: a continuous convolution that operates on functions, and a discrete convolution that operates on sequences. This article discusses only the discrete convolution operator. The discrete convolution operator is a binary operation on two sequences, u and v. It produces a new sequence, denoted u*v, by using multiplication and addition of the terms in the original sequences. Each element of u*v is a dot product for a subsequence of elements in u and v.

In math, the two sequences are usually indexed starting at 0. That is, u = {u₀, u₁, ..., u_M-1} and v = {v₀, v₁, ..., v_N-1}. The sequence u has M elements; the sequence v has N elements. The full convolution u*v is a new sequence that has M+N-1 elements (indexed from 0 to M+N-2). The k-th element of u*v, denoted as (u*v)(k), is defined as:
$(u*v)(k) = \sum_{j=\max(0, k-N+1)}^{\min(k, M-1)} u(j)v(k-j)$

If you adopt the convention that the sequences are infinite and that the unspecified terms of the sequences are zero, this formula can be written more simply:
$(u*v)(k) = \sum_{j=-\infty}^{\infty} u(j)v(k-j)$
This version makes it look more like the continuous convolution operator, which replaces the summation with an integral. The convolution is commutative, so you can also use u(k-j)v(j) in the formula. I use this fact in the section about time-series smoothing.

Computing a convolution: The DFCONV function in SAS IML

In Appendix A, I show how to "manually" compute a discrete convolution between two vectors in the SAS IML language. (The SAS IML language uses one-based indexing, so the indices in the implementation are slightly different than in the formula that I showed in the previous section.) But, you don't need to write any code to use a convolution. The SAS IML language provides a built-in function, the DFCONV function, which automatically computes the convolution and a few useful variations. The documentation shows the formula for vectors that use one-based indexing. The 'DF' prefix stands for 'Digital Filter', which is the name that electrical engineers use for operations that process digital signals.

Let's look at a simple numeric example. Suppose that u = {1,2,3} and v = {4,5,6,7,8} are vectors of length M=3 and N=5, respectively. You can use the DFCONV function to compute the convolution y = u*v, as follows:

proc iml;
u = {1,2,3};
v = {4,5,6,7,8};
y = dfconv(u,v);
print y;

Appendix A provides the step-by-step manual computations that verify that the output is correct. Notice that the output length is N + M - 1 = 3 + 5 - 1 = 7.

Example 1: The sum of two independent random variables

An example of a discrete convolution in probability theory is the distribution of the sum of two independent rolls of dice. When you sum two independent random variables, the distribution of the sum is the convolution of their individual probability distributions.

A classic example is two fair six-sided dice. Let X and Y be the outcomes of the two dice. Both follow a discrete uniform distribution. The probability mass function (PMF) for a single die can be represented as a vector of probabilities for the outcomes 1 through 6. Thus, u = {1/6, 1/6, 1/6, 1/6, 1/6, 1/6} and v = {1/6, 1/6, 1/6, 1/6, 1/6, 1/6} are the PMFs.

Let Z = X + Y be the sum of the dice faces. The random variable Z takes on values in the range 2, 3, 4, ..., 12. The probability that Z=2 is relatively small because there is only one way ("snake eyes") for that sum to occur. In contrast, there are six equally likely ways for the sum to be 7: 1+6, 2+5, 3+4, 4+3, 5+2, and 6+1. In general, if the sum of the two dice is s, and if X=k, then we know that second die has the value s-k. Thus, to find the probability that Z=s, we must sum the probabilities that X=k and Y=s-k for k=1..6.

Because X and Y are independent, the probability that X=k and Y=s-k is the product P(X=k)*P(Y=s-k). To find the complete probability that Z=s, we sum over all the ways that the sum s can occur. This leads to the formula
$P(Z = s) = \sum_{k} P(X = k) P(Y = s - k)$
The right side of the formula is the k_th element of a convolution. Thus, to find the probability distribution of the sum of two dice, you can convolve the PMFs for each die.

You can use the DFCONV function in IML to find the distribution of the sum of two six-sided dice, as follows:

/* Example: Discrete convoution of rolling two fair six-sided dice */
proc iml;
prob_d6 = {1,1,1,1,1,1} / 6;            /* PMF of six-side fair die */
prob_2dice = dfconv(prob_d6, prob_d6);  /* PMF of sum of two dice */
print prob_2dice[r=(2:12) F=FRACT12.];

The PMF for each die has six elements. Thus, the convolution is a vector of length 6 + 6 - 1 = 11, which corresponds to the sums 2-12. The PMF of the sum is the familiar triangular distribution, which peak at the sum of 7 with a probability of 6/36 or 1/6.

If you want the sum of three dice, you can convolve the distribution of the sum with the distribution of another fair die, as follows:

/* Want three dice? Convolve the probability distribution for a six-sided 
   die with the distribution for two independent dice. */
prob_3dice = dfconv(prob_d6, prob_2dice);
print prob_3dice[r=(3:18) F=FRACT12.];

You can continue this process as many times as you like. You can compare the result to the classical formula for rolling n six-sided dice. If the dice are not fair, you can replace the vectors u and v with other vectors that specify the probability for each face.

Example 2: Moving averages in time series

A second application of convolution is smoothing a sequence of data points by using a moving window analysis. When you apply a moving average, you convolve a short sequence (the filter or kernel) with a longer data sequence. I have previously discussed how to compute a moving average and how to compute a moving average in SAS by using PROC EXPAND or SAS IML functions.

In the previous post, I use PROC EXPAND to find two different (backward) moving averages of the monthly price of IBM stock in the late 1990s and early 2000s. The following SAS statements smooth 120 prices and create a graph that overlays the scatter plot, a moving average, and a weighted moving average. See the previous article for an explanation and discussion.

/* Example from https://blogs.sas.com/content/iml/2016/01/27/moving-average-in-sas.html */
title "Monthly IBM Stock Price";
proc sort data=sashelp.stocks(where=(STOCK='IBM') rename=(Date=t Close=y)) 
          out=Series(keep=t y);
  where ('01JAN1996'd <= t <= '31DEC2005'd);
  by t;
run;
 
/* create three moving average curves. Use the TRIM option to remove 
   edge effects and compute only the "middle" of the smoother. */
proc expand data=Series out=out method=none;
   id t;
   convert y = MA   / transout=(movave 5 trim 5);          /* backward MA(5) */
   convert y = WMA  / transout=(movave(1 2 3 4 5) trim 5); /* backward WMA   */
run;
 
proc sgplot data=out cycleattrs;
   scatter x=t y=y;
   series x=t y=MA   / name='MA'   legendlabel="MA(5)";
   series x=t y=WMA  / name='WMA'  legendlabel="WMA(1,2,3,4,5)";
   keylegend 'MA' 'WMA';
   xaxis display=(nolabel) grid;
   yaxis label="Closing Price" grid;
run;

PROC EXPAND has various options to deal with averages at the beginning and end of the sequence, where there are not enough points to perform a full convolution. By default, it adjusts the width of the moving window. To avoid a discussion of edge effects, I used the TRIM option to output missing values for the first and last five averaged values. Thus, the graph shows smoothers that do not predict values for the first five and the last five observations. I will do the same after convolving the window and the data. In this way, I avoid a distracting discussion about why the default behavior of PROC EXPAND deviates from the mathematical definition of a convolution.

Let's perform the moving average computation in SAS IML by using a convolution. We want the filter elements to sum to unity, A simple backward moving average with a window width of 5 corresponds to a filter with values u={1,1,1,1,1}/5. The backward weighted moving average corresponds to u={5,4,3,2,1}/15. Notice that the order of the elements in the weighted kernel are reversed from the PROC EXPAND order. Why? Well, the mathematical definition of a convolution flips the kernel shape (that is, it uses u(k-j) to compute the k_th smoothed point). We therefore need to manually reverse the elements in u before calling the DFCONV function so that the output agrees with the backward moving average from PROC EXPAND.

/* perform a similar computation by using convolution */
proc iml;
/* read the data into v */
use Series;  read all var "y" into v;  close;
 
u = {1,1,1,1,1} / 5;   /* MA(5) */
MA = dfconv(u, v);
/* reverse the filter elements to match the PROC EXPAND output for a backward MA */   
u = {5,4,3,2,1} / 15;  /* WMA(1,2,3,4,5) */
WMA = dfconv(u, v);
 
/* We want to compare with PROC EXPAND. The PROC
   1) outputs N elements.
   2) has special handling for the first and last N obs.
*/
M = nrow(u);
N = nrow(v);
MA  = MA[1:N];    /* keep only N elements */
WMA = WMA[1:N];
 
firstObs = 1:M;
lastObs  = (N-M+1):N; 
MA [firstObs||lastObs] = .;  /* set the first/last obs to missing */
WMA[firstObs||lastObs] = .;
 
create IML_MA var {'MA' 'WMA'};  append;  close;
QUIT;
 
proc compare base=out(keep=MA WMA) compare=IML_MA 
             method=absolute criterion=1E-12;
run;

                                       Observation Summary                                        
 
                                  Observation      Base  Compare                                  
 
                                  First Obs           1        1                                  
                                  Last  Obs         120      120                                  
 
                 Number of Observations in Common: 120.                                           
                 Total Number of Observations Read from WORK.OUT: 120.                            
                 Total Number of Observations Read from WORK.IML_MA: 120.                         
 
                 Number of Observations with Some Compared Variables Unequal: 0.                  
                 Number of Observations with All Compared Variables Equal: 120.

The IML values (computed by using a convolution) agree with the PROC EXPAND values to within 1E-12. (In fact, the largest difference for this example is 5.68E-14.) If we had not used the TRIM option in PROC EXPAND, then there would be differences in the first few observations because PROC EXPAND uses a varying-width kernel near the beginning and end of the data series.

Summary

This article discusses the mathematical definition of a discrete convolution operator. In SAS, you can use the DFCONV function in SAS IML to convolve two vectors of arbitrary lengths. The result is a new longer vector. The convolution operator is illustrated by using two examples: the distribution of the sum of two six-sided dice, and the smoothing of a time series by using a moving average.

Appendix A: Manual computation of a discrete convolution

This appendix shows a simple manual computation of discrete convolution for two vectors. It is presented so that you can compare the output of the DFCONV function with a function that directly implements the convolution formula.

Suppose you define a vector, u={1,2,3}, which has M=3 elements, and a vector, v={4,5,6,7,8}, which has N=5 elements. If you want to compute their convolution, y=u*v, you can use the following manual computations. These computations assume zero-based indexing. For this problem, the convolution has M+N-1 = 7 elements:

y[0] = u[0]v[0] = 1*4 = 4
y[1] = u[0]v[1] + u[1]v[0] = (1*5) + (2*4) = 13
y[2] = u[0]v[2] + u[1]v[1] + u[2]v[0] = (1*6) + (2*5) + (3*4) = 28
y[3] = u[0]v[3] + u[1]v[2] + u[2]v[1] = (1*7) + (2*6) + (3*5) = 34
y[4] = u[0]v[4] + u[1]v[3] + u[2]v[2] = (1*8) + (2*7) + (3*6) = 40
y[5] = u[1]v[4] + u[2]v[3] = (2*8) + (3*7) = 37
y[6] = u[2]v[4] = 3*8 = 24

Thus, the convolution y=u*v is a new vector, y = {4,13,28,34,40,37,24}.

There are two issues when turning this into a computation. First, the IML language uses one-based indexing, not zero-based. Second, a naive implementation would use a pair of nested loops and perform scalar multiplication and addition, which is inefficient in a vectorized language. A better implementation is to use a dot product to perform the sums of the products of subsequences. After the first few elements of the vector v, the dot products will always use subsequences of length M, which is the length of u.

The following computation implements the definition of the discrete convolution formula.

proc iml;
/* y = conv(u,v) convolves column vectors u and v. 
   Algebraically, convolution is the same operation as 
   multiplying the polynomials whose coefficients are the 
   elements of u and v. 
*/
start conv_manual(u,v);
   m = nrow(u);
   n = nrow(v);
   ny = m + n - 1;
   y = j(ny, 1, 0);
 
   do k = 1 to ny;
      idx = max(1, k+1-n):min(k,m); 
      y[k] = u[idx]` * v[k + 1 - idx];
   end;
   return (y);
finish;
 
/* test the computation */
u = {1,2,3};
v = {4,5,6,7,8};
y1 = conv_manual(u,v);
/* compare to the built-in IML function */
y2 = dfconv(u,v);
print y1 y2;

Notice the statement, idx = max(1, k+1-n):min(k,m), which selects the indices to use for the dot product. This statement corresponds to the summation limits in the mathematical formula earlier. In the formula, the lower limit of summation is max(0, k-N+1) and the upper limit is min(k, M-1). The IML statement adjusts these limits to account for 1-based indexing.

Appendix B: Convolution as matrix-vector multiplication

Incidentally, the first time I saw a convolution operator was in a course in Linear Algebra where it was presented as an example of a bilinear operator between two vector spaces. Accordingly, you can represent a convolution as a matrix-vector multiplication where the matrix T=T(v) is a Toeplitz-type matrix of the appropriate dimension and padded with zeros. In honor of my now-deceased Cornell math professor, Keith Dennis, I'll show how to construct the matrix and use it to compute a convolution in IML:

ny = nrow(v) + nrow(u) - 1;
Tv = j(ny, nrow(u), 0);
do i = 1 to nrow(u);
   Tv[i:(i+nrow(v)-1),i] = v;    /* create a Toeplitz-type matrix by copying v into columns */
end;
y3 = Tv*u;                       /* matrix-vector multiplication */
print Tv, y3;

In practice, there is no advantage to construct a large matrix and perform matrix-vector multiplication. You can save RAM memory and compute a series of dot products instead. Nevertheless, it's a cool result!

The post What is convolution and why is it important in statistics? appeared first on The DO Loop.

Owen's T function is defined as a one-dimensional integral, so this article uses the QUAD subroutine in SAS IML to evaluate it. However, it is worth mentioning that Patefield, M. and Tandy, D. (2000, JSS) provide a faster and more accurate algorithm. Their algorithm is more complex. It evaluates the function by selecting from among six different numerical methods, depending on the specified values of (x, α).

The definition of Owen's T function

The definition of Owen's T function, T(x, α), uses α as the upper limit of an integral and X as a parameter to the integrand:
$T(x, \alpha) = \frac{1}{2\pi} \int_0^\alpha \frac{e^{-\frac{x^2}{2}(1+t^2)}}{1+t^2} dt$

From this definition, it is clear that T is an even function in x and an odd function in α:

• T is symmetric in x: T(-x, α) = T(x, α)
• T is anti-symmetric in α: T(x, -α) = -T(x, α)

It is less obvious that Owen's T function simplifies if x=0 or α=1:

• T(0,α) = 1/(2 π) tan^-1(α)
• T(x,1) = 1/2 * Φ(x) (1 - Φ(x)), where Φ(x) is the standard univariate normal CDF.

Define Owen's T function in SAS

The QUAD subroutine in SAS IML is a general-purpose integration routine. The domain of integration in the formula is [0, α]. However, α can be negative, so it might be necessary to reverse the limits of integration. To handle that possibility, I specify the domain of integration as [0, |α|] and negate the answer if α < 0. Here is one way to define Owen's T function:

proc iml;
/* Define the Integrand for Owen's T function:
   The QUAD subroutine passes a scalar value, t, during the integration process.
   Use the GLOBAL statement to pass the parameter value, x (via 'g_x').
*/
start OwenT_Integrand(t) global(g_x);
   t2 = 1 + t##2;
   return ( exp(-0.5 * g_x##2 * t2) / t2 );
finish;
 
/* Define Owen's T Function by integrating OwenT_Integrand on [0, alpha].
   Loop over the elements of x to allow vector inputs. alpha is a scalar.
*/
start OwenT(x, alpha) global(g_x);
   pi = constant("pi");
   n = nrow(x);
   y = j(nrow(x), 1, .);
   /* make sure limits of integration are increasing. If alpha < 0, reverse sign of answer */
   limits = 0 || abs(alpha);
   do i = 1 to n;
      g_x = x[i];   /* set parameter, 'x' */
      call quad(val, "OwenT_Integrand", limits);
      y[i] = val;
   end;
   return( choose(alpha>=0, y/(2*pi), -y/(2*pi)) ); /* negate answer if alpha < 0 */
finish;
 
/* test evaluating the function for x in [-3,3] */
x = T(do(-3, 3, 6/50));
alpha = 1;
T = OwenT(x, alpha);
title "Graph of Owen's T Function: alpha=1";
call series(x, T) grid={x y};

The plot shows the graph of Owen's T function when the parameter α is fixed at α=1 and the parameter x is in the interval [-3, 3]. You can see that the graph looks somewhat bell-shaped.

The following panel of graphs shows the T(x, α) function for various values of α. You can see that the shape of the function is similar for all α values, but the amplitude of the function depends on α. You can show analytically that the magnitude of the T function is greatest at x=0. The height of the function at x=0 is 1/(2 π) tan^-1(α).

If you like 3-D graphs, the Boost C++ library contains a graph of T(x, α) over a rectangular domain.

Summary

Owen's T function is defined as a one-dimensional integral with two parameters. This article uses the QUAD subroutine in SAS IML to evaluate the integral for any parameter values. A subsequent article will discuss applications of Owen's T function in probability and statistics.

The post What is Owen's T function? appeared first on The DO Loop.

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.