# 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 Aug 21 '14 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. – Glen_b -Reinstate Monica Aug 21 '14 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). – Glen_b -Reinstate Monica Aug 21 '14 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. – Glen_b -Reinstate Monica Aug 21 '14 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))