Logistic function macro for SAS

Used in logistic regression, neural networks, and countless other applications, the logistic function is simple to understand yet can be quirky to calculate in SAS because of how SAS treats extreme values. First, let’s consider the value -1000:

R calculates this as expected:

> 1/(1+exp(-(-1000)))
[1] 0
# note: the [1] is not part of the value, so the value is 0

However, SAS gives a missing value:

1676  data _null_;
1677      logistic = 1/(1+exp(-(-1000)));
1678      put logistic=;
1679  run;

NOTE: Invalid argument to function EXP at line 1677 column 21.
logistic=.
logistic=. _ERROR_=1 _N_=1
NOTE: Missing values were generated as a result of performing an operation on
      missing values.
      Each place is given by: (Number of times) at (Line):(Column).
      1 at 1677:20
NOTE: Mathematical operations could not be performed at the following places.
      The results of the operations have been set to missing values.
      Each place is given by: (Number of times) at (Line):(Column).
      1 at 1677:21

The problem is that quickly approaches infinity. Computers have limited floating float precision, and in case of this error, SAS chooses to return a missing value.

One workaround is this macro which places bounds on the extreme values:

/* Logistic macro which is robust to extreme values.*/
/* Depending on your platform and application, you may want to tune the magic number 500. */
%macro logistic(z);
	1/(1+exp(min(max(-(&z),-500),500)))
%mend;

This is code to test the macro with extreme and normal values:

1719  /* demonstrate */
1720  data _null_;
1721      g00=%logistic(-1000); /* should be close to 0.0 */
MPRINT(LOGISTIC):  1/(1+exp(min(max(-(-1000),-500),500)))
1722      g02=%logistic(-4); /* should be close to 0.02 */
MPRINT(LOGISTIC):  1/(1+exp(min(max(-(-4),-500),500)))
1723      g05=%logistic(0); /* should be 0.5 */
MPRINT(LOGISTIC):  1/(1+exp(min(max(-(0),-500),500)))
1724      g09=%logistic(5); /* should be close to 0.9 */
MPRINT(LOGISTIC):  1/(1+exp(min(max(-(5),-500),500)))
1725      g10=%logistic(1000); /* should be close to 1.0 */
MPRINT(LOGISTIC):  1/(1+exp(min(max(-(1000),-500),500)))
1726      put _all_;
1727  run;

g00=7.12458E-218 g02=0.01798621 g05=0.5 g09=0.9933071491 g10=1 _ERROR_=0 _N_=1

If you paid the big bucks for the SAS/GRAPH license, draw a scatterplot to check the non-extreme values:

data logistic;
	do x = -6 to 6 by 0.01;
		y = %logistic(x);
		output;
	end;
run;

title 'Logistic function';
proc gplot data=logistic;
	plot y*x;
run;

This code was tested with (a) Base SAS 9.1.3 32-bit on Windows 7 and (b) Enterprise Guide 4.3 on Windows 7 connected to Base SAS 9.2 64-bit on Server 2008.

Advertisements