Goal

We have empirical data and want to calculate probabilities based on these data. We should be aware that we need many data points, otherwise the probability estimates will be incorrect. The more data, the better the approximation of the true distribution! We will start with some explanations. If you are in a hurry, you can just use the code provided and the end of this document.

Data simulation

We start with simulating some example data using the rnorm function of the R base package. These data represents some random variable that we call uppercase \(X\). The variable \(X\) could for instance be the day temperature of a randomly chosen day in Alaska (Sorry, I didn’t have a better idea. Do you?).

We randomly choose 300 values from a normal distribution with mean = 2 and standard deviation = 3:

set.seed(77)  # guarantees identical random data when 'rnorm' command is repeated
dat = rnorm(300, mean = 2, sd = 3)   
str(dat)                             
##  num [1:300] 0.351 5.273 3.919 5.128 2.509 ...

Short remark: Why do we label the random variable with \(X\)? Hmmm, because it is much easier to write for the probability \(P(X \le a)\) instead of \(P(\text{day temperature} \le a)\). Of course, you can use the latter if you like. We live in a free country. (Still).

Histogram and density estimation

We use the hist function. We store the output of the function in a variable in order to see what hist calculates:

histo = hist(dat, col = "lightblue", breaks = 25, freq = FALSE, xlab = "x")  # this plots already, but we save in histo anyway

str(histo)
## List of 6
##  $ breaks  : num [1:20] -9 -8 -7 -6 -5 -4 -3 -2 -1 0 ...
##  $ counts  : int [1:19] 1 0 1 2 7 4 11 17 29 35 ...
##  $ density : num [1:19] 0.00333 0 0.00333 0.00667 0.02333 ...
##  $ mids    : num [1:19] -8.5 -7.5 -6.5 -5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 ...
##  $ xname   : chr "dat"
##  $ equidist: logi TRUE
##  - attr(*, "class")= chr "histogram"

The setting “freq = FALSE” is important because it normalises the total area of the histogram to 1. You don’t trust it? Okay, let us check this:

width_of_bins = histo$breaks[2] - histo$breaks[1] # all distances are equal, see str(histo) above 
area = sum(histo$density)*width_of_bins           
area
## [1] 1

We get the density function using density in R.

PDF <- density(dat)  # estimated Probability Density Function (PDF)
str(PDF)
## List of 8
##  $ x         : num [1:512] -10.8 -10.7 -10.7 -10.6 -10.6 ...
##  $ y         : num [1:512] 1.75e-05 2.05e-05 2.39e-05 2.77e-05 3.21e-05 ...
##  $ bw        : num 0.847
##  $ n         : int 300
##  $ old.coords: logi FALSE
##  $ call      : language density.default(x = dat)
##  $ data.name : chr "dat"
##  $ has.na    : logi FALSE
##  - attr(*, "class")= chr "density"

We plot the histogram and the density together:

hist(dat, col = "lightblue", breaks = 25, freq = FALSE, xlab = "x")
lines(PDF$x, PDF$y, col = 2, lwd = 2)  # "lines" adds to an existing plot

This looks pretty nice, I guess.

Estimation of the CDF

The Cumulative Distribution Function (CDF) is a key concept in statistics. Often, it is just called distribution function.
The CDF of a random variable \(X\) is designated \(F(x)\) or even \(F_X(x)\). Here, the uppercase \(X\) denotes the variable, and the lowercase \(x\) is the value we want to calculate \(F\) for (this is often confusing, I know …).
\(F_X(x)\) is the probability that the outcome of \(X\) will take a value less than or equal to \(x\), i.e. \(F_X(x) = P(X <= x)\). The value of \(F_X(x)\) is the same as the area under the PDF to the left of \(x\). (Don’t mix up CDF and PDF, here). I tried to visualise the connection between PDF and CDF in this document. Take a look at it!

If we once have \(F_X(x)\), we can calculate the probabilities using the following expressions:

\[ \begin{aligned} P \left( a < X \leq b \right) &= F_X(b) - F_X(a) \\ P \left( X > a \right) &= 1 - F_X(a) \\ P \left( X \le a \right) &= F_X(a) \end{aligned} \]

The formulas denote (in that order):

In R, a CDF can be estimated from empirical data using the function ecdf.

CDF <- ecdf(dat) 
plot(CDF, main = "CDF", xlab = "x", ylab = "F(x)", col = "blue", cex = 0.5, font.main = 1)

That looks pretty okay. The CDF must start at 0 at the extreme left position, and end up at 1 at the extreme right position.

Calculation of probabilities from the CDF

Now, knowing the CDF and the three formulas listed above, we can calculate whatever probability we want. Let’s look at a few examples:

a = 0
b = 2
Prob_a_b = CDF(b) - CDF(a) # probability that X ends up between a and b
round(Prob_a_b,2)
## [1] 0.27
a = 2
Prob_a = CDF(a) # probability that X ends up lower than or equal to a 
Prob_a
## [1] 0.5066667
a = 2
Prob_b = 1 - CDF(a) # probability that X ends up greater than  a 
Prob_b
## [1] 0.4933333

Because we have chosen the same point \(a\) in the last two calculations, the sum of the respective probabilities must be 1 (why?).

Prob_a + Prob_b
## [1] 1

The promised short version

We make it as short as possible but include data simulation so that you have a working example. If you need plots, you can take the code above.

set.seed(77)  # garantees same random data when command is repeated
dat = rnorm(300, mean = 2, sd = 3)
CDF <- ecdf(dat)

We have already all what we need.

Prob_a_b = CDF(2) - CDF(0) # probability that X ends up between 0 and 2
round(Prob_a_b,2)
## [1] 0.27
Prob_a = CDF(2)            # probability that X ends up lower than or equal to 2 
round(Prob_a, 2)
## [1] 0.51
Prob_b = 1 - CDF(2)        # probability that X ends up greater than  2 
round(Prob_b, 2)
## [1] 0.49
Prob_a + Prob_b            # just a quick check 
## [1] 1

Graphics with ggplot

If you prefer ggplot for graphics (I don’t), you can get it from here (R-script).

Note that all code snippets are designed for comprehensibility, not for efficiency!


uwe.menzel@matstat.org