Kernel density estimates for positive data

The kernel density estimate (KDE) is a powerful tool for estimating the density of univariate data. The KDE is a flexible model. For example, it can fit data distributions that are multimodal or have long tails. It is nonparametric, which means that you do not need to assume any form for the distribution of the data. In contrast, parametric estimates (such as normal or lognormal distributions) impose constraints on the shape of the data distribution.

However, the KDE's flexibility comes at a cost. If the underlying distribution is known to have certain properties, the KDE does not necessarily preserve those properties. A property that arises often in practice is positivity. Many real-world measurements are strictly positive, such as mass, length, and many health-related clinical lab values. If you have a distribution of lengths, and you fit a KDE to the data, it is disconcerting (and wrong!) to display a density estimate that predicts positive density for negative lengths!

There are two ways to deal with this problem. This article shows how to use SAS to truncate the density estimate by setting the density to zero when the variable is negative. A subsequent article shows how to log-transform the data, fit a KDE on the log scale, and back-transform the results in SAS.

Fit a KDE to positive data

The following data set shows the position (in millimeters) of a rivet relative to a bulkhead. The distance is a positive quantity. For convenience, I assign the name of the data sets and the variables to macro variables. The remainder of the examples use only the macro variables, which makes it easy for you to run the code on your own data.

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

You can use either PROC UNIVARIATE or PROC KDE to fit a kernel density estimate. I will use PROC KDE because, as we will see in the next section, it supports options to truncate the KDE. The following call fits the KDE to the data and visualizes the result. The UNISTATS option displays descriptive statistics for the data and shows the bandwidth chosen for the fit:

proc kde data=&DSName;
   univar &varName / unistats;
run;

I have highlighted a few entries in the output tables. Although the range of the data is [1.5, 21.28], the Controls table shows that the range of the KDE is [-4.99, 27.77]. The graph of the KDE curve confirms this fact. The KDE extends about three bandwidths past the extremes of the data, and the kernel bandwidth for these data is 2.13. However, it does not make sense to estimate the density to be nonzero for Position < 0 because Position must be a positive quantity.

PROC KDE provides options that enable you to truncate the density estimate in two ways:

Truncate the KDE to the range of the data by using the TRUNCATE option.
Truncate the KDE to an arbitrary range by using the GRIDL= option to specify a lower limit and/or the GRIDU= option to specify an upper limit.

It is somewhat unfortunate that PROC KDE does not have an option to display the histogram and KDE on the percentage scale. Most graphs in this article are shown on the frequency scale. You can rescale the graphs to display them on the percent or density scale. The last graph in this article is shown on the percentage scale.

Truncate the density into the data range

The following call uses the TRUNCATE option to restrict the KDE to the range of the data:

/* truncate KDE to data range */
proc kde data=&DSName;
   univar &varName / truncate out=KDE_truncData; /* save KDE to data set for later analysis */
run;

Notice that the KDE is truncated on the interval [1.5, 21.28], which is the range of the data.

Truncate the density onto the positive half-line

The Position variable cannot be negative, but there is no intrinsic reason why it cannot be smaller than 1.5. There is also no obvious reason to bound the density from above. An alternative way to truncate the density is to use the GRIDL=0 option to truncate the density to strictly positive values.

proc kde data=&DSName;
   univar &varName / gridl=0 out=KDE_trunc0; /* save KDE to data set for later analysis */
run;

This KDE is restricted to positive values. It is unbounded above, but the density is exactly zero for negative values of the Position variable. If you have an upper limit, you can set it by using the GRIDU= option. For example, if you are plotting the density of student scores on an exam that has 100 as a perfect score, you could use GRIDL=0 GRIDU=100 to bound the KDE to the range [0, 100].

Limitations of the truncation method

There are two problems with truncating a KDE:

On a probability scale, the area under a true density must integrate to 1. If you tell PROC KDE to truncate the KDE, it does so by throwing away some area. Thus, area is less than 1 for the restricted KDE.
Even if the estimate is rescaled, the KDE underestimates the expected value of distribution near 0 (Silverman, 1986, Section 2.10).

You can rescale to the density to correct for the first problem. In the previous sections, I wrote the values of the KDE to the data sets KDE_truncData and KDE_trunc0, respectively. The following macro computes the area under these curves by using a trapezoidal approximation to the area.

/* the area under a density estimate should be 1. See
   https://blogs.sas.com/content/iml/2011/07/08/the-area-under-a-density-estimate-curve-nonparametric-estimates.html
   The following macro computes the area under a KDE estimate that is created by PROC KDE and written
   to the 'DSName' data set by using the OUT= option on the UNIVAR statement.
   In the data set the x variable is named 'value', and the 'density' variable is f(x).
   */
%macro AreaUnderKDE(DS);
data _null_;
   set &DS end=EOF;
   dx = dif(value);
   meanY = (density + lag(density)) / 2;
   Area + dx*meanY;
   if EOF then do;
      call symputx('Area', Area, "G"); /* create global macro */
   end;
run;
%mend;
 
%AreaUnderKDE(KDE_truncData);
%put Trunc Data &=Area;
%AreaUnderKDE(KDE_trunc0);
%put Trunc to 0 &=Area;

Trunc Data AREA=0.9031115838
Trunc to 0 AREA=0.97522813

The log shows that the area is less than 1 for both KDE curves. This can be a problem if you intend to use the KDE to estimate probabilities. However, it is simple enough to rescale the truncated density by dividing it by its area:

/* If 'density' is a variable that does not integrates to 1,
   you can renormalize by dividing by the Area under the density curve. */
%macro Renormalize(DS, Area);
   data &DS;
      set &DS;
      density = density / &Area;
   run;
%mend;
 
%Renormalize(KDE_trunc0, &Area);

The following graph overlays the renormalized truncated KDE on a histogram and fringe plot of the data. The KDE is truncated at 0 because the measured quantity must be positive. The graph was created by using the %EmulateHistogram macro and the process described in the article, "Overlay multiple custom density curves on a histogram in SAS."

The second issue (underestimating the KDE near the boundary) does not have a simple solution. Some researchers suggest reflecting the data across x=0 before fitting the KDE. I don't want to say more about that method here, but I have used it previously to create smoothers of periodic data.

Summary

If you naively compute a kernel density estimate for data that are inherently positive, the KDE might estimate density for negative values of the variable. PROC KDE in SAS provides two methods for correcting this problem:

The TRUNCATE option on the UNIVAR statement truncates the KDE to the range of the data.
The GRIDL=L and GRIDU=U options truncate the KDE to the range [L, U]. For example, you can use GRIDL=0 to truncate the KDE onto the positive axis.

Technically, if you truncate the KDE, you should renormalize it so that the area under the curve is unity. This article includes two macros (%AreaUnderKDE and %Renormalize) that you can use for that purpose. I do not deal with the fact that the truncated KDE underestimates the expected density near the boundary. In a subsequent article, I present an alternative method: log-transforming the data.

Blogs