Manipulation of Discrete Random Variables with discreteRV

A prominent issue in statistics education is the sometimes large disparity between the theoretical and the computational coursework. discreteRV is an R package for manipulation of discrete random variables which uses clean and familiar syntax similar to the mathematical notation in introductory probability courses. The package offers functions that are simple enough for users with little experience with statistical programming, but has more advanced features which are suitable for a large number of more complex applications. In this paper, we introduce and motivate discreteRV, describe its functionality, and provide reproducible examples illustrating its use. Introduction One of the primary hurdles in teaching probability courses in an undergraduate setting is to bridge the gap between theoretical notation from textbooks and lectures, and the statements used in statistical software required in more and more classes. Depending on the background of the student, this missing link can manifest itself differently: some students master theoretical concepts and notation, but struggle with the computing environment, while others are very comfortable with statistical programming, but find it difficult to translate their knowledge back to the theoretical setting of the classroom. discreteRV (Buja et al., 2015) is an approach to help with bringing software commands closer to the theoretical notation. The package provides a comprehensive set of functions to create, manipulate, and simulate from discrete random variables. It is designed for introductory probability courses. discreteRV uses syntax that closely matches the notation of standard probability textbooks to allow for a more seamless connection between a probability classroom setting and the use of statistical software. discreteRV is available for download on the Comprehensive R Archive Network (CRAN). discreteRV was derived from a script written by Dr. Andreas Buja for an introductory statistics class (Buja). The package rv2 (Buja and Wickham, 2014), available on GitHub, provides a useful example of using devtools (Wickham and Chang, 2015) to begin basic package development, and also uses Dr. Buja’s code as a starting point. The goal of rv2 seems more focused on learning package development, while the goal of discreteRV is to be a useful statistics education and probability learning tool. The functions of discreteRV are organized into two logical areas, termed probabilities and simulations. This document will illustrate the use of both sets of functions. All code used in this document is available in a vignette, accessible by loading discreteRV and calling vignette("discreteRV").


Introduction
One of the primary hurdles in teaching probability courses in an undergraduate setting is to bridge the gap between theoretical notation from textbooks and lectures, and the statements used in statistical software required in more and more classes.Depending on the background of the student, this missing link can manifest itself differently: some students master theoretical concepts and notation, but struggle with the computing environment, while others are very comfortable with statistical programming, but find it difficult to translate their knowledge back to the theoretical setting of the classroom.
discreteRV (Buja et al., 2015) is an approach to help with bringing software commands closer to the theoretical notation.The package provides a comprehensive set of functions to create, manipulate, and simulate from discrete random variables.It is designed for introductory probability courses.discreteRV uses syntax that closely matches the notation of standard probability textbooks to allow for a more seamless connection between a probability classroom setting and the use of statistical software.discreteRV is available for download on the Comprehensive R Archive Network (CRAN).discreteRV was derived from a script written by Dr. Andreas Buja for an introductory statistics class (Buja).The package rv2 (Buja and Wickham, 2014), available on GitHub, provides a useful example of using devtools (Wickham and Chang, 2015) to begin basic package development, and also uses Dr. Buja's code as a starting point.The goal of rv2 seems more focused on learning package development, while the goal of discreteRV is to be a useful statistics education and probability learning tool.
The functions of discreteRV are organized into two logical areas, termed probabilities and simulations.This document will illustrate the use of both sets of functions.All code used in this document is available in a vignette, accessible by loading discreteRV and calling vignette("discreteRV").

Probabilities
discreteRV includes a suite of functions to create, manipulate, and compute distributional quantities for discrete random variables.A list of these functions and brief discriptions of their functionality is available in Table 1.

Creating random variables
The centerpiece of discreteRV is a set of functions to create and manipulate discrete random variables.A random variable X is defined as a theoretical construct representing the value of an outcome of a random experiment (see e.g.Wild and Seber, 1999).A discrete random variable is a special case that can only take on a countable set of values.Discrete random variables are associated with probability mass functions, which map the set of possible values of the random variable to probabilities.Probability mass functions must therefore define probabilities which are between zero and one, and must sum to one.
Throughout this document, we will work with two random variables, a simple example of a discrete random variable representing the value of a roll of a fair die, and one representing a realization of a Poisson random variable with mean parameter equal to two.Formally, we can define such random variables and their probability mass functions as follows: Let X be a random variable representing a single roll of a fair die; i.e., the sample space Ω = {1, 2, 3, 4, 5, 6} and X is the identity, mapping the upturned face of a die roll to the corresponding number of dots visible.
Let Y be a random variable distributed according to a Poisson distribution with mean parameter λ.In this case, Y takes on values in the non-negative integers {0, 1, 2, . ..}.Then, In discreteRV, a discrete random variable is defined through the use of the RV function.RV accepts a vector of outcomes, a vector of probabilities, and returns an 'RV' object.The code to create X, a random variable representing the roll of a fair die, is as follows: > (X <-RV(outcomes = 1:6, probs = 1/6))

Displaying first 12 outcomes
Several common distributions are natively supported so that the functions need not be defined manually.For instance, an equivalent method of defining Y is: Random variable with outcomes from 0 to Inf Outcomes 0 1 2 3 4 5 6 7 8 9 10 11 Probs 0.135 0.271 0.271 0.180 0.090 0.036 0.012 0.003 0.001 0.000 0.000 0.000

Displaying first 12 outcomes
The RV function also allows the definition of a random variable in terms of odds.We construct a loaded die in which a roll of one is four times as likely as any other roll as: Random variable with 6 outcomes Outcomes 1 2 3 4 5 6 Odds 4:5 1:8 1:8 1:8 1:8 1:8

Structure of an 'RV' object
The syntactic structure of the included functions lends itself both to a natural presentation of elementary probabilities and properties of probability mass functions in an introductory probability course, as well as more advanced modeling of discrete random variables.In Table 2, we provide an overview of the notational similarities between discreteRV and the commonly used probability textbook by Casella and Berger (2001).
The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 A random variable object is constructed by defining a standard R vector to be the possible outcomes that the random variable can take (the sample space Ω).It is preferred, though not required, that these be encoded as numeric values, since this allows the computation of expected values, variances, and other distributional properties.This vector of outcomes then stores attributes which include the probability of each outcome.By default, the print method for a random variable will display the probabilities as fractions in most cases, aiding in readability.The probabilities can be retrieved as a numeric vector by using the probs function: > probs(X) 1 2 3 4 5 6 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

Probability-based calculations
By storing the outcomes as the principle component of the object X, we can make a number of probability statements in R. For instance, we can calculate the probability of obtaining a roll greater than 1 by using the code P(X > 1).R will check which of the values in the vector X are greater than 1.In this case, these are the outcomes 2, 3, 4, 5, and 6.Hence, R will return TRUE for these elements of X, and we compute the probability of this occurrence in the function P by simply summing over the probability values stored in the names of these particular outcomes.Likewise, we can make slightly more complicated probability statements such as P(X > 5 ∪ X = 1), using the %OR% and %AND% operators.Consider our Poisson random variable Y, and suppose we want to obtain the probability that Y is within a distance δ of its mean parameter λ = 2: > delta <-3; lambda <-2 > P((Y >= lambda -delta) %AND% (Y <= lambda + delta)) [1] 0.9834364 Alternatively, we could have also used the slightly more complicated looking expression: > P((Y -lambda)^2 <= delta^2) [1] 0.9834364 Conditional probabilities often provide a massive hurdle for students of introductory probability classes.These types of questions often make it necessary to first translate the problem from everyday language into the mathematical concept of conditional probability, e.g., what is the probability that you will not need an umbrella when the weather forecast said it was not going to rain?Similarly, what is the probability that a die shows a one, if we already know that the roll is no more than 3? The mathematical solution is, of course, P(X = 1 | X ≤ 3).In discreteRV this translates directly to a solution of P(X == 1 | X <= 3).The use of the pipe operator may be less intuitive to the seasoned R programmer, but overcomes a major notational issue in that conditional probabilities are most commonly specified with the pipe.Using the pipe for conditional probablity, we had to create alternative %OR% and %AND% operators, as specified previously.
We can compute several other distributional quantities, including the expected value and the variance of a random variable.In notation from probability courses, expected values can be found with the E function.To compute the expected value for a single roll of a fair die, we run the code E(X).
The expected value for a Poisson random variable is its mean, and hence E(Y) in our example will return the value two.The function V(X) computes the variance of random variable X.Alternatively, we can also work from first principles and assure ourselves that the expression E((X -E(X))ˆ2) provides the same result:

Joint distributions
Aside from moments and probability statements, discreteRV includes a powerful set of functions used to create joint probability distributions.Once again letting X be a random variable representing a The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859  single die roll, we can use the iid function to compute the probability mass function of n trials of X. Table 3 gives the first eight outcomes for n = 2, and Table 4 gives the first eight outcomes for n = 3.Notice again that the probabilities have been coerced into fractions for readability.Notice also that the outcomes of the joint distribution are encoded by the outcomes on each trial separated by a comma.Outcome 1,1,1  1,1,2  1,1,3  1,1,4  1,1,5  1,1,6  1,2,1  1,2,  The * operator has been overloaded in order to allow a more seamless syntax for defining joint distributions.Suppose we wish to compute the joint distribution of X, our toss of a fair coin, and a coin flip.After defining the coin flip variable, the joint distribution can be defined as follows: > Z <-RV(0:1) > X * Z
Joint distributions need not be the product of iid random variables.Joint distributions in which the marginal distributions are dependent can also be defined.Consider the probability distribution defined in Table 5.Note that A and B are dependent, as the product of the marginal distributions does not equal the joint.We can define such a random variable in discreteRV by using the jointRV function, which is a wrapper for RV: > (AandB <-jointRV(outcomes = list(1:3, 0:2), probs = 1:9 / sum(1:9)))
The SofI function computes the random variable representing the sum of two independent, but not necessarily identically distributed, random variables.The + operator is overloaded to make this computation even more syntactically friendly.Note, however, that similar limitations apply as in the joint distribution case:

Plotting
discreteRV includes a plot method for random variable objects so that visualizing outcomes and their probabilities is as simple as calling plot(X).Figure 1 on the left shows a visual representation of the probability mass function (pmf) of a fair die.The x axis includes all outcomes, and the y axis includes the probabilities of each particular outcome.Figure 1 on the right shows the pmf of the sum of two independent rolls of a fair die.The pmf of a sum of 20 independent rolls of a die is given in Figure 2 on the left.
In addition to a plotting method, there is also a method for qqnorm to allow assessment of normality for random variable objects, as displayed in Figure 2 on the right.While very close to a normal, the sum of 20 independent rolls of a fair die still shows a slight S curve in the Q-Q plot.

Simulation
discreteRV also includes a set of functions to simulate trials from a random variable.A list of these functions and brief discriptions of their functionality is available in Table 7.

Creation
Creating a simulated random vector is done by using the rsim function.rsim accepts a parameter X representing the random variable to simulate from, and a parameter n representing the number of independent trials to simulate.For example, suppose we'd like to simulate ten trials from a fair die.We have already created a random variable object X, so we simply call rsim as follows: The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 Name Description plot Plot method for a simulated random vector, i.e., an object of class 'RVsim' Prop Proportion of an event observed in a vector of simulated trials props Proportions of observed outcomes in one or more vectors of simulated trials rsim Simulate n independent trials from a random variable X skewSim Skew of the empirical distribution of simulated data > (X.sim <-rsim(X, 10)) Simulated Vector: 4 3 2 3 4 2 2 3 6 5 Random variable with 6 outcomes Outcomes 1 2 3 4 5 6 Probs 1/6 1/6 1/6 1/6 1/6 1/6 The object returned is a vector of simulated values, with an attribute containing the random variable that was used for the simulation.If we would like to retrieve only the simulated values and exclude the attached probabilities, we can coerce the object into a vector using R's built-in as.vector function.
> as.vector(X.sim) [1] 4 3 2 3 4 2 2 3 6 5 It is also possible to retrieve some quantities from the simulation.We can retrieve the empirical distribution of simulated values with the props function.This will return the outcomes from the original random variable object, and the observed proportion of simulated values for each of the outcomes.We can also compute observed proportions of events by using the Prop function.Similar to the P function for probability computations on random variable objects, Prop accepts a variety of logical statements.

Extended example: playing Craps
Craps is a common dice game played in casinos.The game begins with what is called the "Come Out" roll, in which two fair dice are rolled.If a sum of seven or eleven is obtained, the player wins.If a sum of two, three, or twelve is obtained, the player loses.In all other cases, the roll obtained is declared the "Point" and the player rolls again in an attempt to obtain this same point value.If the player rolls the Point, they win, but if they roll a seven, they lose.Rolls continue until one of these two outcomes is achieved.
discreteRV allows for a seamless analysis and simulation of the probabilities associated with different events in Craps.Let us begin by asking "What is the probability that the game ends after the first roll?"To answer this question we construct two random variables.We note that calling RV(1:6) returns a random variable for a single roll of a fair die, and then we use the overloaded + operator to sum over two rolls to obtain the random variable Roll.> (Roll <-RV(1:6) + RV(1:6)) The Recall that the game ends after the first roll if and only if a seven or eleven is obtained (resulting in a win), or a two, three, or twelve is obtained (resulting in a loss).Hence, we calculate the probability that the game ends after the first roll as follows: > P(Roll %in% c(7, 11, 2, 3, 12)) [1] 0.3333333 Now suppose we would like to condition on the game having ended after the first roll.Using the conditional probability operator in discreteRV, we can obtain the probabilities of winning and losing given that the game ended after the first roll: > P(Roll %in% c(7, 11) | Roll %in% c(7, 11, 2, 3, 12)) [1] 0.6666667 > P(Roll %in% c(2, 3, 12) | Roll %in% c(7, 11, 2, 3, 12)) [1] 0.3333333 Now, let's turn our attention to calculating the probability of winning a game in two rolls.Recall that we can use the iid function to generate joint distributions of independent and identically distributed random variables.In this case, we would like to generate the joint distribution for two independent rolls of two dice.Now, we will have 11 2 possible outcomes, and our job is to determine which outcomes result in a win.We know that any time the first roll is a seven or eleven, we will have won.
We also know that if the roll is between four and ten inclusive, then we will get to roll again.To win within two rolls given that we've received a four through ten requires that the second roll matches the first.We can enumerate the various possibilities to calculate the probability of winning in two rolls, which is approximately 30%.
> TwoRolls <-iid(Roll, 2) > First <-marginal(TwoRolls, 1) > Second <-marginal(TwoRolls, 2) > P(First %in% c(7, 11) %OR% (First %in% 4:10 %AND% (First == Second))) [1] 0.2993827 Finally, suppose we are interested in the empirical probability of winning a game of Craps.Using the simulation functions in discreteRV, we can write a routine to simulate playing Craps.Using the rsim function, we simulate a single game of Craps by rolling from our random variable Roll, which represents the sum of two dice.We then perform this simulation 100000 times.

Figure 2 :
Figure 2: Left: plot method called on a sum of 20 fair die random variables; right: qqnorm method called on a sum of 20 fair die random variables.

Table 1 :
Overview of functions provided in discreteRV ordered by topics.

Table 2 :
Probability functions in discreteRV and their corresponding syntax in introductory statistics courses.

Table 3 :
First eight outcomes and their associated probabilities for a variable representing two independent rolls of a die.

Table 4 :
First eight outcomes and their associated for a variable representing three independent rolls of a die.

Table 6 :
Outcomes and their associated probabilities for a variable representing the sum of two independent rolls of a die.

Table 7 :
List of the simulation functions contained in discreteRV.
R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859 The results indicate that the player wins a game of craps approximately 49% of the time.