5
$\begingroup$
set.seed(42)
a <- rnorm(200) # generate random data
hist(a, probability = T) # histogram
d <- density(a) # calculate density
polygon(d, density = 10, col = "deeppink") # plot density
abline(v=0, col="blue") # draw boundaries
abline(v=0.1, col="blue")

If I want to know the probability of the variable between 2 points (since it is considered continuous), e.g. 0 and 0.1, the way to calculate it is:

integrate(approxfun(d), lower = 0, upper = 0.1)

Since the probability density function is the derivative of the cumulative distribution function, when I apply integrate(approxfun(d), lower = 0, upper = 0.1), I am reversing the derivation. Does this mean that the result I get is cumulative probability?

$\endgroup$
0

2 Answers 2

6
$\begingroup$

For easier reading, I have combined three extensive Comments (now deleted) into an Answer:

You don't have the true PDF $f(x)$ from density in R. From the code, we know $X$ is standard normal, so the exact value of $p=P(0<X<1)$ could be found from in R as $0.3413447.$

diff(pnorm(c(0,1)))
[1] 0.3413447

However, I suppose you want to get $p$ from your $n=200$ observations a. The most direct way to do that is to find the proportion of values of a in (0,1):

set.seed(42)
a = rnorm(200) # generate random data

mean((a > 0) & (a < 1)) 
[1] 0.33

Alternatively, if you know data are normal, then you could estimate $μ,σ$ from data and use R's normal pdf function to get $0.3429.$

mu = mean(a);  sd = sd(a)
diff(pnorm(c(0,1), mu, sd)) 
[1] 0.3428732

I wouldn't expect a density estimator to do very much better.

The output of density in R is a sequence of 512 x-values and 512 y-values that can be used to plot the estimated PDF (enclosing unit area).

density(a)

Call:
        density.default(x = a)

Data: a (200 obs.);     Bandwidth 'bw' = 0.2895

       x                 y            
 Min.   :-3.8616   Min.   :0.0000781  
 1st Qu.:-2.0036   1st Qu.:0.0124148  
 Median :-0.1456   Median :0.0685148  
 Mean   :-0.1456   Mean   :0.1344194  
 3rd Qu.: 1.7124   3rd Qu.:0.2425435  
 Max.   : 3.5704   Max.   :0.4123541 

hist(a, prob=T, col="skyblue2")
 lines(density(a), col="brown", lwd=2)
 rug(a)

The figure below shows a histogram of a along with the density estimator. Tick marks along the horizontal axis show locations of the $n=200$ observations. [Sometimes density estimators are informally called 'smoothed histograms', but they are based on individual data points without reference to the binning of any histogram. The density estimator used here is the default estimator from density in R; variations are available via parameters not used here.]

enter image description here

You might try to use this output to estimate $p,$ as follows, to get $p \approx 0.337867.$

xx = density(a)$x; yy = density(a)$y
sum(yy[xx > 0 & xx < 1])/sum(yy)
[1] 0.337867

This method does have the advantage of not needing to know the population family of distributions (e.g., normal).


Addendum, showing results for a much larger sample: $n=10\,000.$

set.seed(2022)

# counting points
A = rnorm(10000)
mean((A > 0) & (A < 1))
[1] 0.3396

# assuming normality 
mu = mean(A);  sd = sd(A)
diff(pnorm(c(0,1), mu, sd))
[1] 0.3407714

# density estimation 
xx = density(A)$x;  yy = density(A)$y
sum(yy[xx > 0 & xx < 1])/sum(yy)
[1] 0.3373269
$\endgroup$
1
  • $\begingroup$ BruceET, thanks for the concise and extensive response - I greatly appreciate that you took the time to write this response. $\endgroup$
    – Reader 123
    May 10 at 19:05
7
$\begingroup$

The cumulative density function (cdf, $F(x)$) and the probability density function (pdf, $f(x)$) are related by the equation:

$F(x) = P(X \leq x) = \int_{-\infty}^{x} f(x)dx$

By the fundamental theorem of calculus, we have that $\frac{\partial}{\partial x}F(x) = f(x)$

(We can replace $-\infty$ with the lower limit of the distribution if that distribution doesn't have support extending to $-\infty$.)

For continuous distributions, we can show that $P(X = a) = 0$ for all points $a$. So we can be a little flexible with our uses of $<$ vs $\leq$.

The probability that $X$ takes a value in a range is defined as

$P(a < X < b) = \int_{a}^{b}f(x)dx$

Which is also equal to

$P(a < X \leq b) = F(b) - F(a)$

So when we take an integral $\int_{a}^{b}f(x)dx$, we are not getting a cumulative probability, except in the case that $a = -\infty$ (or the lower limit of the distribution).

$\endgroup$
3
  • $\begingroup$ Thanks a lot for taking the trouble to answer, @David Luke Thiessen, I appreciate it. I had seen the formulas while I was researching it - I guess the part I find difficult is that I assumed that the difference of two cdfs (F(b) - F(a) must be a cdf itself... $\endgroup$
    – Reader 123
    May 10 at 19:08
  • 1
    $\begingroup$ @Reader123 One of the required properties of a CDF is that as x goes to infinity, F(x) goes to 1. So unfortunately, transformations of the CDF usually result in something that isn't a CDF. Instead it's easier to make transformations inside the P(X <= x) term and find the new CDF from that. $\endgroup$ May 12 at 22:47
  • $\begingroup$ Thank you, @David Luke Thiessen - very useful, I didn't think of this! $\endgroup$
    – Reader 123
    May 15 at 12:15

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.