# How to find/estimate probability density function from density function in R

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


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.

• 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. – Glen_b -Reinstate Monica Dec 5 '13 at 23:08
• @Glen_b I want to use the estimated density for finding the probability of other values in the distribution. For instance, I have a vector of data ranging from 0 to 10.This data set contains only 70 unique values between 0 and 10. I can plot the density. Now suppose that I'm interested on finding the probability of having X=7.5, which is not in observed data, in a random sampling.How can I get it? I know that ecdf(X) gives me the equivalent percentile of 7.5 but it's not what I'm looking for. – Amin Dec 5 '13 at 23:59
• "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). – Glen_b -Reinstate Monica Dec 6 '13 at 6:17
• 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. – Glen_b -Reinstate Monica Dec 6 '13 at 7:03
• 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. – Glen_b -Reinstate Monica Dec 6 '13 at 7:14

?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)


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.

• 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? – Amin Dec 6 '13 at 16:33
• 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'. – Glen_b -Reinstate Monica Dec 6 '13 at 17:03
• 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! – Amin Dec 6 '13 at 18:04
• 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. – Glen_b -Reinstate Monica Dec 7 '13 at 0:01