32
$\begingroup$

Suppose that I have a variable like X with unknown distribution. In Mathematica, by using SmoothKernelDensity function we can have an estimated density function.This estimated density function can be used alongside with PDF function to calculate probability density function of a value like X in the form of PDF[density,X] assuming that "density" is the result of SmoothKernelDensity. It would be good if there is such feature in R.This is how it works in Mathematica

http://reference.wolfram.com/mathematica/ref/SmoothKernelDistribution.html

As an example (based on Mathematica functions):

data = RandomVariate[NormalDistribution[], 100]; #generates 100 values from N(0,1)

density= SmoothKernelDistribution[data]; #estimated density

PDF[density, 2.345] returns 0.0588784 

Here you can find more information about PDF:

http://reference.wolfram.com/mathematica/ref/PDF.html

I know that I can plot its density function using density(X) in R and by using ecdf(X) I can obtain its empirical cumulative distribution function.Is it possible to do the same thing in R based on what I described about Mathematica?

Any help and idea is appreciated.

$\endgroup$
19
  • $\begingroup$ density(x) gives an estimate of the pdf, as you already noted, but its suitability depends on the purpose for which you want to have the density. Note, for example, that the variance is biased up (in performing convolution, you add the variance of the kernel to the variance of the data, itself an unbiased estimate) - such bias-variance tradeoffs are ubiquitous. There are other alternatives, such as log-spline density estimation, for example -- but again, its suitability partly depends on what you want to do with it. $\endgroup$
    – Glen_b
    Commented Dec 5, 2013 at 23:08
  • 1
    $\begingroup$ "finding the probability of having X=7.5" -- there's your problem! Either you have a continuous distribution (in which case the actual answer is "0"), or you don't (in which case you shouldn't be using density estimation, because you don't have a density). $\endgroup$
    – Glen_b
    Commented Dec 6, 2013 at 6:17
  • 2
    $\begingroup$ Note the definition of the ecdf (or the cdf more generally); ecdf(b)-ecdf(a) would estimate $P(\text{a}<X\leq \text{b})$. Of course with continuous variables the distinction between $<$ and $\leq$ is unimportant. If $X$ is discrete, then you can estimate $P(X=7.5)$ by computing the sample proportion of values that are 0.75. $\endgroup$
    – Glen_b
    Commented Dec 6, 2013 at 7:03
  • 1
    $\begingroup$ Sorry, that was an error. I mean the sample proportion of values that are 7.5; my son distracted me as I was typing the last couple of words. Your sample estimate of the probability of an unobserved event is zero. Did you want to apply a prior? Did you want a confidence interval for the proportion instead of a point estimate? Your actual issue is not yet an R issue, your issue is in correctly explaining what it is you actually want. You should probably edit your question, or post a new one. $\endgroup$
    – Glen_b
    Commented Dec 6, 2013 at 7:14
  • 1
    $\begingroup$ @Glen_b regardless of our discussion, my question is clear. I want to know if it's possible to find probability based on estimated density (assume for a continuous variable) in R as what we can do in Mathematica. $\endgroup$
    – Amin
    Commented Dec 6, 2013 at 15:50

2 Answers 2

46
$\begingroup$

?density points out that it uses approx to do linear interpolation already; ?approx points out that approxfun generates a suitable function:

x <- log(rgamma(150,5))
df <- approxfun(density(x))
plot(density(x))
xnew <- c(0.45,1.84,2.3)
points(xnew,df(xnew),col=2)

enter image description here

By use of integrate starting from an appropriate distance below the minimum in the sample (a multiple - say 4 or 5, perhaps - of the bandwidth used in df would generally do for an appropriate distance), one can obtain a good approximation of the cdf corresponding to df.

$\endgroup$
5
  • $\begingroup$ this is interesting. It seems that df(2.3) gives the value of estimated density function at x=2.3 but what PDF does in Mathematica is giving the area under curve above x=2.3. I'm not quite sure about this.This is just my guess.Can you re-produce what I did in Mathematica? $\endgroup$
    – Amin
    Commented Dec 6, 2013 at 16:33
  • $\begingroup$ My function above demonstrably gives a kernel-based estimate of a "probability density function" ... "evaluated at x". Either you want that, or you don't. If you don't, you have to explain what you do want - in statistical terms, not just as 'reproduce this behavior'. $\endgroup$
    – Glen_b
    Commented Dec 6, 2013 at 17:03
  • $\begingroup$ I think that I mistakenly and unintentionally have promoted that density is probability which is not. I didn't mean to be misleading.If you think that PDF in Mathematica does what you described in your answer (i.e. finding the value of density function for given X value) then I think that I got my answer. Just there are many confusion on using words! $\endgroup$
    – Amin
    Commented Dec 6, 2013 at 18:04
  • 2
    $\begingroup$ From what the PDF page says it does, it returns the same kind of thing I do, but the methods it uses in its calculation in this case are likely to be somewhat more accurate (for such a purpose additional accuracy has little value, however). For some discussion of the probability/density distinction, see here and here. $\endgroup$
    – Glen_b
    Commented Dec 7, 2013 at 0:01
  • 1
    $\begingroup$ Because density samples equally spaced values, cumsum is a more expedient way to construct points along the CDF than integration might be. Linear interpolation then finishes the job of obtaining a function for the CDF. $\endgroup$
    – whuber
    Commented Nov 24, 2021 at 19:34
1
$\begingroup$

spatstat.core::CDF() can be used to to create a cumulative density function from a given output from density().

set.seed(123)
x <- rnorm(10000000)

x_density <- density(x, n = 10000)

x_cdf <- spatstat.core::CDF(x_density)

sds <- c(-2, -1, 0, 1, 2)
names(sds) <- sds

# check cdf at different values
setNames(
  x_cdf(sds), 
  sds)
#>         -2         -1          0          1          2 
#> 0.02285086 0.15889356 0.50009332 0.84134448 0.97717762

# compare against theoretical
pnorm(sds)
#>         -2         -1          0          1          2 
#> 0.02275013 0.15865525 0.50000000 0.84134475 0.97724987

Created on 2021-11-22 by the reprex package (v2.0.0)

Update

A previous version of this answer copied code from the deprecated spatstat:::CDF() which was broken up (in ?2020?) into several other packages. If anyone knows a lighter weight package where this CDF function currently exists would love to hear about it in the comments!

$\endgroup$
3
  • 1
    $\begingroup$ When z references the object returned by density, a CDF corresponding to z can be created easily via f <- with(z, splinefun(x, cumsum(y) / sum(y), method="monoH.FC")). That's pretty lightweight ;-). $\endgroup$
    – whuber
    Commented Nov 24, 2021 at 15:11
  • 1
    $\begingroup$ @whuber that looks pretty great. I was also pointed to the {logspline} package. I did a quick benchmark and seems that the splinefun or approxfun approaches for interpolation are both way faster: gist.github.com/brshallo/ea2e04347e14fae7ff969a54e2266359 $\endgroup$ Commented Nov 24, 2021 at 18:26
  • $\begingroup$ I wrote a toy package here: github.com/brshallo/densdpqr that uses whuber's approach for doing this and other density --> distribution operations. $\endgroup$ Commented Jan 5, 2022 at 23:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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