A previous article discusses the problem of fitting a kernel density estimate (KDE) to data that are strictly positive. If you use a standard KDE, the resulting density curve might estimate non-zero probability for negative values. This is unsatisfactory because quantities like lengths and mass cannot be negative. The previous article shows how to use PROC KDE in SAS to truncate the KDE (and optionally renormalize the density).

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

The Log-KDE method

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

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

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

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

Implement the method in SAS

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

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

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

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

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

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

Visualize the back-transformed KDE

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

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

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

Peculiarities of the log-KDE method

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

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

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

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

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

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

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

Summary

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

Further Reading

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