This R worksheet does not include assessed questions.You should complete this worksheet in Week 8. If you have difficulty with this worksheet, you can get help at the office hours drop-in sessions.

On the last worksheet, we looked at dealing with famous discrete random variables: the binomial, geometric, and Poisson distributions. In this worksheet, we will look at dealing with arbitrary discrete random variables, as defined by a probability mass function (PMF)

Probability mass function

The key to dealing with a discrete random variable is to set up two things: the range of values the random variable can take, and the PMF itself.

Consider, for example, the following PMF for a random variable \(X\):

\(x\) \(0\) \(0.5\) \(1\) \(2\) \(2.5\) \(4\)
\(p_X(x)\) \(0.1\) \(0.2\) \(0.2\) \(0.3\) \(0.1\) \(0.1\)

We want to store the \(x\) values in one vector, and the \(p_X(x)\) values in another vector. We can do that like this:

x <- c(0, 0.5, 1, 2, 2.5, 4)
pmf_x <- c(0.1, 0.2, 0.2, 0.3, 0.1, 0.1)

It’s good practice to check that the PMF does indeed sum to 1:

sum(pmf_x)
## [1] 1

We can find the probability of a given value of \(x\) by pulling out the relevant value of the PMF. We are used to using square brackets to pull out the elements of a vector. But more complicate expressions can be put inside the brackets, like this:

Exercise 8.1. Let \(Y\) be the sum of fair dice rolls. So \(Y\) takes the values \((2, 3, \dots, 12)\) with probabilities \(\frac{1}{36}(1,2,3,4,5,6,5,4,3,2,1)\).
(a) Enter this random variable into y and pmf_y. Do so as efficiently as possible, by taking advantage of the colon notation.
(b) What is \(\mathbb P(Y = 9)\)?
(c) What is \(\mathbb P(Y \leq 6)\)?

Plotting a PMF

If we want to plot a picture of a PMF (like those in the Section 6 of the notes) we can use plot() together with the optional argument type = "h". Here, "h" stands for “histogram-like vertical lines”, which are most appropriate for displaying PMFs. I also like to make the lines somewhat thicker, with the lwd = ... argument, and make sure that the y-axis starts from 0.

So I would draw a graph for our previous PMF like this:

plot(x, pmf_x, type = "h", lwd = 4, ylim = c(0, max(pmf_x)))

Exercise 8.2. Draw a plot of the PMF for the dice random variable \(Y\) from Exercise 8.1. Draw the lines in blue. Make sure to label the axes appropriately.

Cumulative distribution function

We know that a CDF \(F(x) = \mathbb P(X = x)\) can be calculated by summing up the PDF; that is, \[ F(x) = \sum_{y \leq x} p(x) . \] To calculate this in R, it is helpful to use the cumulative sum function cumsum(). The cumulative sum of a vector \[ (a, \quad b, \quad c, \quad d, \quad e) \] is the vector of partial sums \[ (a, \quad a+b, \quad a+b+c, \quad a+b+c+d, \quad a+b+c+d+e) . \]

So taking the cumulative sum of the PMF,

cumsum(pmf_x)
## [1] 0.1 0.3 0.5 0.8 0.9 1.0

gives a vector of the values of the CDF at the values of \(x\).

However, the CDF also should be defined in between the values of \(x\) in the range of \(X\), where it should behave like a “step function”. We can create such a function using the stepfun() function. This stepfun requires two arguments:

So we form the CDF function with

cdf_x <- stepfun(x, c(0, cumsum(pmf_x)))

We can test that this CDF works between the values in the range of \(X\) by testing it:

cdf_x(1.62)
## [1] 0.5

We can also plot the CDF with the plot() function:

plot(cdf_x)

with other extra graphical commands, as required. The “blobs” are to let the reader know that at the jumps, the function takes the higher value.

Exercise 8.3. Form the CDF for the dice random variable from Exercise 8.1.
(a) What is \(F_Y(4.6)\)?
(b) What is \(\mathbb P(Y > 8.2)\)?
(c) Plot the CDF. Make sure to label your axes appropriately.

Expectation

Recall that the definition of the expectation of a discrete random variable is \[ \mathbb EX = \sum_x x\, p(x) . \]

So to calculate this, we can simply use R’s sum() function. We just want

exp_x <- sum(x * pmf_x)
exp_x
## [1] 1.55

That’s all there is to it!

Exercise 8.4. With the dice random variable again:
(a) Calculate the expectation \(\mathbb EY\).
(b) Calculate \(\mathbb EY^2\).
(c) Using parts (a) and (b), find the variance of \(Y\).

Sampling from a distribution

Finally, we may wish to sample from a discrete random variable. To do this, we will use the sample() function, like this:

sample(x, 20, replace = TRUE, prob = pmf_x)
##  [1] 2.0 2.0 0.0 2.0 2.0 0.5 0.0 2.0 0.5 2.0 2.0 1.0 0.0 2.0 2.0 0.5 2.5 1.0 2.0 1.0

Here, the arguments are:

  1. The first argument is the set of values to be sampled, which for us is the range of the random variable.
  2. The second argument is the number of samples we wish to take. In the above example, we took 20 samples.
  3. replace = TRUE tells us we want to be allowed “repeats” in our sample. If we want IID samples of a random variable, we always want sampling with replacement. (For some sampling-without-replacement probabiltiy problems, it might be approriate to take replace = FALSE, which is R’s default for the sample() function.)
  4. prob = ... tells R with what probabilities each value of x should be sampled. If this argument is omitted, R takes each value to be equally likely.

Exercise 8.5. Sample 10,000 simulated dice rolls from the random variable \(Y\). Calculate the mean of these simulations. Compare your answer with that from Exercise 8.4(a), and explain why you got the answer you did, using results from the module.