SAS博客列表

I previously discussed how to use the PUTLOG statement to write a message from the DATA step to the log in SAS. The PUTLOG statement is commonly used to write notes, warnings, and errors to the log. This article shows how to use the PRINTTOLOG subroutine in SAS IML software to write messages to the log. The PRINTTOLOG subroutine was released in 2020 in SAS IML on SAS Viya. It runs in PROC IML and in the iml action on Viya. The PRINTTOLOG subroutine is not part of SAS 9.4, but I show how to construct a SAS/IML module for SAS 9.4 that has the same syntax and functionality.

The PrintToLog subroutine in SAS IML on Viya

The syntax of the subroutine is
    call PrintToLog(msg, flag);
The flag parameter is optional. If you omit that argument, the message is displayed in the log. If you use the flag parameter, the valid values are 0, 1, and 2:

• If flag is 0, the string "NOTE: " is prepended to the message. For example, call PrintToLog("This is a note.", 0);.
• If flag is 1, the string "WARNING: " is prepended to the message. For example, call PrintToLog("This is a warning.", 1);.
• If flag is 2, the string "ERROR: " is prepended to the message. For example, call PrintToLog("This is an error.", 2);.

I will illustrate the routine by using an IML version of the program that I used in my previous post about the PUTLOG statement. The input to the program is a three-column matrix. Each row of the matrix contains the coefficients (a, b, c) for the quadratic equation a x² + b x + 3 = 0. The program is vectorized. It uses the discriminant of the quadratic equation to determine whether the quadratic equation has any real roots. If so, the program uses the quadratic formula to find the roots. To demonstrate the PrintToLog function, the program does the following:

• Display an ERROR if any coefficient of the quadratic term is 0.
• Display a WARNING if any quadratic equation does not have any real roots.
• Display a NOTE if any quadratic equation has a repeated root.

proc iml;
/* NOTE: The PrintToLog subroutine is available in SAS IML in Viya. 
   It is not available as a built-in subroutine for SAS/IML in SAS 9.4. */
Q = { 1 -1 -6,    /* each row is (a,b,c) */
      1  7  6,
      0  7  6,
      2  7  7,
     -2  5 12,
     -2  4 -2,
      5  4 10};
a = Q[,1]; b = Q[,2]; c = Q[,3];
Root1 = j(nrow(a), 1, .); /* initialize roots to missing */
Root2 = j(nrow(a), 1, .); /* initialize roots to missing */
 
if any(a=0) then
   call PrintToLog("The coefficient of the quadratic term is 0 for at least one observation", 2); /* display ERROR */
 
discrim = b##2 - 4*a#c;
if any(discrim < 0) then 
   call PrintToLog("At least one equation does not have real roots.", 1); /* display WARNING */
if any(discrim=0) then 
   call PrintToLog("At least one equation has repeated roots.", 0); /* display NOTE */
 
idx = loc(discrim >= 0 && a^=0);
if ncol(idx)>0 then do;
	Root1[idx] = (-b[idx] - sqrt(discrim[idx])) / (2*a[idx]);
	Root2[idx] = (-b[idx] + sqrt(discrim[idx])) / (2*a[idx]);
end;
print a b c discrim Root1 Root2;

The SAS log contains the following message, which is color-coded according to whether the message is an error (red color), a warning (green color), or a note (blue color).

For completeness, the output of the program is shown below.

A PrintToLog module for SAS 9.4

Although the PrintToLog subroutine is not a built-in subroutine in SAS 9.4, you can define a module that has the same functionality. The definition uses two SAS/IML techniques and one macro trick:

• You can use the %PUT macro statement to write a message to the log from a SAS/IML program.
• You can use string concatenation and the EXECUTE subroutine in SAS/IML to run a global SAS statement that requires a literal string (not a SAS/IML variable).
• You can use the SYSVER system macro to determine whether you are executing a program in SAS 9.4 or SAS Viya. In Viya, the value of the SYSVER macro is V.04.00. In SAS 9.4, the value is 9.4. Thus, you can use conditional logic to create a macro that defines a SAS/IML module in SAS 9.4 but does nothing in SAS Viya.

The following SAS macro defines a SAS/IML module in SAS 9.4, but not in SAS Viya. Thus, you can run this macro in any PROC IML program. If the program runs in Viya (where the PrintToLog subroutine is already defined), the macro does nothing. If the program runs in SAS 9.4, the macro defines a PrintToLog module that you can use to write messages to the log.

/* Check the SYSVER macro to see if SAS 9.4 is running.
   In SAS Viya, the macro is empty and does nothing.
   In SAS 9.4, the macro defines a function that emulates the PrintToLog call.
   The syntax is as follows:
   call PrintToLog("This is a log message.");
   call PrintToLog("This is a note.", 0);
   call PrintToLog("This is a warning.", 1);
   call PrintToLog("This is an error.", 2);
*/
%macro DefinePrintToLog;
%if %sysevalf(&sysver = 9.4) %then %do;
start PrintToLog(msg,errCode=-1);
   if      errCode=0 then prefix = "NOTE: ";
   else if errCode=1 then prefix = "WARNING: ";
   else if errCode=2 then prefix = "ERROR: ";
   else prefix = "";
   stmt = '%put ' + prefix + msg + ';';
   call execute(stmt);
finish;
%end;
%mend;
 
/* this program runs in SAS 9.4 or in SAS Viya */
proc iml;
%DefinePrintToLog;
call PrintToLog("This is a log message.");
call PrintToLog("This is a note.", 0);
call PrintToLog("This is a warning.", 1);
call PrintToLog("This is an error.", 2);

The output is shown for SAS 9.4, where the %DefinePrintToLog macro defines the SAS/IML module.

Summary

The PRINTTOLOG statement is way for SAS IML programmers to write messages to the log. In SAS Viya, this subroutine is built-in and can be used from PROC IML or from the iml action. This article shows how to define a module for PROC IML in SAS 9.4 that has similar functionality.

The post Write to the log from SAS IML programs appeared first on The DO Loop.

Recently, I needed to know "how much" of a piecewise linear curve is below the X axis. The coordinates of the curve were given as a set of ordered pairs (x1,y1), (x2,y2), ..., (xn, yn). The question is vague, so the first step is to define the question better. Should I count the number of points on the curve for which the Y value is negative? Should I use a weighted sum and add up all the negative Y values? Ultimately, I decided that the best measure of "negativeness" for my application is to compute the area that lies below the line Y=0 and above the curve. In calculus, this would be the "negative area" of the curve. Because the curve is piecewise linear, you can compute the area exactly by using the trapezoid rule of integration.

An example is shown to the right. The curve is defined by 12 ordered pairs. The goal of this article is to compute the area shaded in blue. This is the "negative area" with respect to the line Y=0. With very little additional effort, you can generalize the computation to find the area below any horizontal line and above the curve.

Area defined by linear segments

The algorithm for computing the shaded area is simple. For each line segment along the curve, let [a,b] be the interval defined by the left and right abscissas (X values). Let f(a) and f(b) be the corresponding ordinate values. Then there are four possible cases for the positions of f(a) and f(b) relative to the horizontal reference line, Y=0:

• Both f(a) and f(b) are above the reference line. In this case, the area between the line segment and the reference line is positive. We are not interested in this case for this article.
• Both f(a) and f(b) are below the reference line. In this case, the "negative area" can be computed as the area of a trapezoid: $A = 0.5 (b - a) (f(b) + f(a))$.
• The value f(a) is below the reference line, but f(b) is above the line. In this case, the "negative area" can be computed as the area of a triangle. You first solve for the location, c, at which the line segment intersects the reference line. The negative area is then $A = 0.5 (c - a) f(a)$.
• The value f(a) is above the reference line and f(b) is below the line. Again, the relevant area is a triangle. Solve for the intersection location, c, and compute the negative area as $A = 0.5 (b - c) f(b)$.

The three cases for negative area are shown in the next figure:

You can easily generalize these formulas if you want the above the curve and below the line Y=t. In every formula that includes f(a), replace that value with (f(a) – t). Similarly, replace f(b) with (f(b) – t).

Compute the negative area

The simplest computation for the negative area is to loop over all n points on the line. For the i_th point (1 ≤ i < n), let [a,b] be the interval [x[i], x[i+1]] and apply the formulas in the previous section. Since we skip any intervals for which f(a) and f(b) are both positive, we can exclude the point (x[i], y[i]) if y[i-1], y[i], and y[i+1] are all positive. This is implemented in the following SAS/IML function. By default, the function returns the area below the line Y=0 and the curve. You can use an optional argument to change the value of the horizontal reference line.

proc iml;
/* compute the area below the line y=y0 for a piecewise
   linear function with vertices given by (x[i],y[i]) */
start AreaBelow(x, y, y0=0);
   n = nrow(x);
   idx = loc(y<y0);                /* find indices for which y[i] < 0 */
   if ncol(idx)=0 then return(0);
 
   k = unique(idx-1, idx, idx+1);  /* we need indices before and after */
   jdx = loc(k > 0 & k <  n);      /* restrict to indices in [1, n-1] */
   v = k[jdx];                     /* a vector of the relevant vertices */
 
   NegArea = 0;
   do j = 1 to nrow(v);            /* loop over intervals where f(a) or f(b) negative */
      i = v[j];                    /* get j_th index in the vector v */
      fa = y[i]-y0; fb = y[i+1]-y0;/* signed distance from cutoff line */
      if fa > 0 & fb > 0 then 
         ;                         /* segment is above cutoff; do nothing */
      else do;
         a  = x[i];  b = x[i+1];
         if fa < 0 & fb < 0 then do;  /* same sign, use trapezoid rule */
            Area = 0.5*(b - a) * (fb + fa);
         end;
         /* different sign, f(a) < 0, find root and use triangle area */
         else if fa < 0 then do;  
            c = a - fa * (b - a) / (fb - fa);
            Area = 0.5*(c - a)*fa;
         end;
         /* different sign, f(b) < 0, find root and use triangle area */
         else do;
            c = a - fa * (b - a) / (fb - fa);
            Area = 0.5*(b - c)*fb;
         end;
         NegArea = NegArea + Area;
      end;
   end;
   return( NegArea );
finish;
 
/* points along a piecewise linear curve */
x = {   1,    2, 3.5,   4,5,     6, 6.5,   7,   8,  10, 12,  15};
y = {-0.5, -0.1, 0.2, 0.7,0.8,-0.2, 0.3, 0.6, 0.3, 0.1,-0.4,-0.6};
 
/* compute area under the line Y=0 and above curve (="negative area") */
NegArea = AreaBelow(x,y);
print NegArea;

The program defines the AreaBelow function and calls the function for the piecewise linear curve that is shown at the top of this article. The output shows that the area of the shaded regions is -2.185.

The area between a curve and a reference line

I wrote the program in the previous section so that the default behavior is to compute the area below the reference line Y=0 and above a curve. However, you can use the third argument to specify the value of the reference line. For example, the following statements compute the areas below the lines Y=0.1 and Y=0.2, respectively:

Area01 = AreaBelow(x,y, 0.1); /* reference line Y=0.1 */ Area02 = AreaBelow(x,y, 0.2); /* reference line Y=0.2 */ print Area01 Area02;

As you would expect, the magnitude of the area increases as the reference line moves up. In fact, you can visualize the area below the line Y=t and the curve as a function of t. Simply, call the AreaBelow function in a loop for a sequence of increasing t values and plot the results:

t = do(-0.7, 1.0, 0.05); AreaUnder = j(1,ncol(t)); do i = 1 to ncol(t); AreaUnder[i] = AreaBelow(x, y, t[i]); end; title "Area Under Reference Line"; call series(t, AreaUnder) label={'Reference Level' 'Area'} grid={x y} xvalues=do(-0.7, 1.0, 0.1);

The graph shows that the area under the reference line is a monotonic function. If the reference line is below the minimum value of the curve, the area is zero. As you increase t, the area below the line Y=t and the curve increases in magnitude. After the reference line reaches the maximum value along the curve (Y=0.8 for this example), the magnitude of the area increases linearly. It is difficult to see in the graph above, but the curve is actually linear for t ≥ 0.8.

Summary

You can use numerical integration to determine "how much" of a function is negative. If the function is piecewise linear, the integral over the negative intervals can be computed by using the trapezoid rule. This article shows how to compute the area between a reference line Y=t and a piecewise linear curve. When t=0, this is the "negative area" of the curve.

Incidentally, this article is motivated by the same project that inspired me to write about how to test whether a function is monotonically increasing. If a function is monotonically increasing, then its derivative is strictly positive. Therefore, another way to test a function for monotonic increasing is to test whether the derivative is never negative. A way to measure how far a function deviates from being monotonic is to compute the "negative area" for the derivative.

The post The area under a piecewise linear curve appeared first on The DO Loop.

A previous article discusses the definitions of three kinds of moments for a continuous probability distribution: raw moments, central moments, and standardized moments. These are defined in terms of integrals over the support of the distribution. Moments are connected to the familiar shape features of a distribution: the mean, variance, skewness, and kurtosis. This article shows how to evaluate integrals in SAS to evaluate the moments of a continuous distribution, and how to compute the mean, variance, skewness, and kurtosis of the distribution.

Notice that this article is not about how to compute the sample estimates. You can easily use PROC MEANS or PROC UNIVARIATE to compute the sample mean, sample variance, and so forth, from data. Rather, this article is about how to compute the quantities for the probability distribution.

Definitions of moments of a continuous distribution

Evaluating a moment requires performing an integral. The domain of integration is the support of the distribution. This is often infinite (for example, the normal distribution on (-∞, ∞)) or semi-bounded (for example, the lognormal distribution on [0,∞)), but sometimes it is finite (for example, the beta distribution on [0,1]). The QUAD subroutine in SAS/IML software can compute integrals on finite, semi-infinite, and infinite domains. You need to specify a SAS/IML module that evaluates the integrand at an arbitrary point in the domain. To compute a raw moment, the integrand is computed as the product of x and the probability density function (PDF). To compute a central moment, the integrand is the product the PDF and a power of (x - μ), where μ is the mean of the distribution. Let f(x) be the PDF. The moments relate to the mean, variance, skewness, and kurtosis as follows:

• The mean is the first raw moment: $\mu = \int_{-\infty}^{\infty} x f(x)\, dx$
• The variance is the second central moment: $\sigma^2 = \int_{-\infty}^{\infty} (x - \mu)^2 f(x)\, dx$
• The skewness is the third standardized moment: $\mu_3 / \sigma^3$ where $\mu_3 = \int_{-\infty}^{\infty} (x - \mu)^3 f(x)\, dx$ is the third central moment.
• The full kurtosis is the fourth standardized moment: $\mu_4 / \sigma^4$ where $\mu_4 = \int_{-\infty}^{\infty} (x - \mu)^4 f(x)\, dx$ is the fourth central moment. For consistency with the sample estimates, you can compute the excess kurtosis by subtracting 3 from this quantity.

Compute moments in SAS

To demonstrate how to compute the moments of a continuous distribution, let's use the gamma distribution with shape parameter α=4. The PDF of the Gamma(4) distribution is shown to the right. The domain of the PDF is [0, ∞). The mean, variance, skewness, and excess kurtosis for the gamma distribution can be computed analytically. The formulas are given in the Wikipedia article about the gamma distribution. The following SAS/IML program evaluates the formulas for α=4.

/* compute moments (raw, central, and standardized) for a probability distribution */
proc iml;
/* Formulas for the mean, variance, skewness, and excess kurtosis, 
   of the gamma(alpha=4) distribution from https://en.wikipedia.org/wiki/Gamma_distribution
*/
alpha = 4;                 /* shape parameter */
mean  = alpha;             /* mean */
var   = alpha;             /* variance */
skew  = 2 / sqrt(alpha);   /* skewness */
exKurt= 6 / alpha;         /* excess kurtosis */
colNames = {'Mean' 'Var' 'Skew' 'Ex Kurt'};
print (mean||var||skew||exKurt)[c=colNames L="Formulas for Gamma"];

Let's evaluate these same quantities by performing numerical integration of the integrals that define the moments of the distribution:

/* return the PDF of a specified probability distribution */
start f(x) global(alpha);
   return  pdf("Gamma", x, alpha);  /* return the PDF evaluated at x */
finish;
 
/* the first raw moment is the mean: \int x*f(x) dx */
start Raw1(x);
   return( x # f(x) );              /* integrand for first raw moment */
finish;
 
/* the second central moment is the variance: \int (x-mu)^2*f(x) dx */
start Central2(x) global(mu);
   return( (x-mu)##2 # f(x) );      /* integrand for second central moment */
finish;
 
/* the third central moment is the UNSTANDARDIZE skewness: \int (x-mu)^3*f(x) dx */
start Central3(x) global(mu);
   return( (x-mu)##3 # f(x) );      /* integrand for third central moment */
finish;
 
/* the fourth central moment is the UNSTANDARDIZE full kurtosis: \int (x-mu)^4*f(x) dx */
start Central4(x) global(mu);
   return( (x-mu)##4 # f(x) );      /* integrand for fourth central moment */
finish;
 
Interval = {0, .I};                 /* support of gamma distribution is [0, infinity) */
call quad(mu, "Raw1", Interval);    /* define mu, which is used for central moments */
call quad(var, "Central2", Interval);
call quad(cm3, "Central3", Interval);
call quad(cm4, "Central4", Interval);
 
/* do the numerical calculations match the formulas? */
skew = cm3 / sqrt(var)##3;          /* standardize 3rd central moment */
exKurt = cm4 / sqrt(var)##4 - 3;    /* standardize 4th standard moment; subtract 3 to get excess kurtosis */
print (mu||var||skew||exKurt)[c=colNames L="Numerical Computations for Gamma"];

For this distribution, the numerical integration produces numbers that exactly match the formulas. For other distributions, the numerical integration might differ from the exact values by a tiny amount.

Checking the moments against the sample estimates

Not all distributions have analytical formulas for the mean, variance, and so forth. When I compute a new quantity, I like to verify the computation by comparing it to another computation. In this case, you can compare the distributional values to the corresponding sample estimates for a large random sample. In other words, if you generate a large random sample from the Gamma(4) distribution and compute the sample mean, sample, variance, and so forth, you should obtain sample estimates that are close to the corresponding values for the distribution. You can perform this simulation in the DATA step or in the SAS/IML language. The following statements use the DATA step and PROC MEANS to generate a random sample of size N=10,000 from the gamma distribution:

data Gamma(keep = x);
alpha = 4;
call streaminit(1234);
do i = 1 to 10000;
   x = rand("Gamma", alpha);
   output;
end;
run;
 
proc means data=Gamma mean var skew kurt NDEC=3;
run;

As expected, the sample estimates are close to the values of the distribution. The sample skewness and kurtosis are biased statistics and have relatively large standard errors, so the estimates for those parameters might not be very close to the distributional values, even for a large sample.

Summary

This article shows how to use the QUAD subroutine in SAS/IML software to compute raw and central moments of a probability distribution. This article uses the gamma distribution with shape parameter α=4, but the computation generalizes to other distributions. The important step is to write a user-defined function that evaluates the PDF of the distribution of interest. For many common distributions, you can use the PDF function in Base SAS to evaluate the density function. You also must specify the support of the distribution, which is often an infinite or semi-infinite domain of integration. When the moments exist, you can compute the mean, variance, skewness, and kurtosis of any distribution.

The post Compute moments of probability distributions in SAS appeared first on The DO Loop.