# Why simulated gamma distributed data have negative kernel values?

I know that Gamma distribution does not allow 0 or negative values. I was doing some simulation and when I write this code in R

set.seed(1234)
Gamma<- rgamma(2032,shape=1.77, rate= 4.33)
hist(Gamma, prob=T,ylim=c(0,2), xlim=c(-1,2.5))
lines(density( Gamma), col='blue')


I obtain this

I can see that the kernel line strart before the axes origin. I'm sure it is correct but can someone explain me why?

• It's not correct; solutions to the problem are given at stats.stackexchange.com/questions/65866/….
– whuber
Commented Aug 21, 2014 at 21:08
• Some discussion of the cause is here. Besides the solutions discussed in answers at whuber's link, there's the solution I offered in comments there - take logs, perform KDE, and then transform the the density estimate back (not forgetting the Jacobian). It's very simple and often quite successful. Commented Aug 21, 2014 at 22:55
• I should add, however that sometimes - especially with gamma-like data, if it's jammed up against 0 - the default bandwidth may need to be increased (by a factor of perhaps 3 or so in some cases), or a different transform might work better (such as a cube root). Commented Aug 21, 2014 at 23:06
• On your data: denlGa=density(log(Gamma),bw=.3) (about double the normal bandwidth), then lines(exp(denlGa$x),denlGa$y/exp(denlGa\$x),col=4) works reasonably well. In some cases the transformation approach isn't so suitable, so use your judgement. Commented Aug 21, 2014 at 23:37

The reason for this is that the command density uses by default a Gaussian kernel, which has support on the whole real line:

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/density.html

A (non-optimal, but enough for practical purposes) estimator of the density for the positive data can be obtained by fitting a KDE to the log of the data and then transforming the KDE back as follows:

set.seed(1234)
Gamma<- rgamma(2032,shape=1.77, rate= 4.33)
hist(Gamma, prob=T,ylim=c(0,2), xlim=c(-1,2.5))

# Log data
hist(log(Gamma), prob=T)
# bandwidth for log-data
hT = bw.nrd0(log(Gamma))
# KDE for log-data
kde <- Vectorize(function(x) mean(dnorm((x-log(Gamma))/hT)/hT))
# Transformed KDE
lkde <- function(x) kde(log(x))/x
# Fit
hist(Gamma, prob=T,ylim=c(0,2), xlim=c(-1,2.5))

• Have you seen the solution at stats.stackexchange.com/a/71291? The unique aspect of the present question is to ask for an explanation why the estimator produces positive density for negative values. Although your explanation is correct, it assumes something the OP probably doesn't know--namely how a kernel density is calculated--and so would benefit from a little more elaboration and perhaps an illustration using a small simple dataset.
– whuber
Commented Aug 21, 2014 at 21:09
• @whuber No, I haven't, but thank you for the reference. It would be interesting to see a comparison of the alternative approach you proposed in your answer. Commented Aug 21, 2014 at 21:12
• I think the one presented by Gavin Simpson is extremely promising; I expect it to perform better than the approach I had posted.
– whuber
Commented Aug 21, 2014 at 21:14
• @whuber Apparently he approximates the log-density of the data using splines. His proposal is analogous to what I suggested (approximating the log-density using KDE). Although I am not sure if my proposal if optimal in any sense (and I am not sure if it is worth investing time on investigating it). I think it is OK (and fast) for practical purposes. Commented Aug 21, 2014 at 21:21
• @Filippo OP=Original Poster (urbandictionary.com/define.php?term=op). In this case, you. You may want to start with the wikipedia entry: en.wikipedia.org/wiki/Kernel_density_estimation. Commented Aug 21, 2014 at 21:24