1
$\begingroup$

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

histogram of simulated data and line of the kernel

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

$\endgroup$
  • 2
    $\begingroup$ It's not correct; solutions to the problem are given at stats.stackexchange.com/questions/65866/…. $\endgroup$ – whuber Aug 21 '14 at 21:08
  • $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Aug 21 '14 at 22:55
  • $\begingroup$ 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). $\endgroup$ – Glen_b -Reinstate Monica Aug 21 '14 at 23:06
  • $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Aug 21 '14 at 23:37
2
$\begingroup$

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))
curve(kde,-6,2,add=T,col="red")
# Transformed KDE
lkde <- function(x) kde(log(x))/x
# Fit
hist(Gamma, prob=T,ylim=c(0,2), xlim=c(-1,2.5))
curve(lkde,0,3,add=T,col="red")
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – whuber Aug 21 '14 at 21:09
  • $\begingroup$ @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. $\endgroup$ – Munchausen Aug 21 '14 at 21:12
  • 1
    $\begingroup$ I think the one presented by Gavin Simpson is extremely promising; I expect it to perform better than the approach I had posted. $\endgroup$ – whuber Aug 21 '14 at 21:14
  • $\begingroup$ @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. $\endgroup$ – Munchausen Aug 21 '14 at 21:21
  • 1
    $\begingroup$ @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. $\endgroup$ – Munchausen Aug 21 '14 at 21:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.