Working with Joint Probability Tables in R

For this assignment we will create a joint probability table and use it to compute marginal and conditional probabilities, expectations and conditional expectations, variances, and pmf’s and CDF’s.

We have two random variables, \(X\) and \(Y\), with their respective sample spaces: \[\Omega_x = \{0,5,10\}\] \[\Omega_y = \{0,5,10,15\}\]

The joint probability table for these random variables is given by:

	\(P(Y=0)\)	\(P(Y=5)\)	\(P(Y=10)\)	\(P(Y=15)\)
\(P(X=0)\)	0.02	0.06	0.02	0.10
\(P(X=5)\)	0.04	0.15	0.20	0.10
\(P(X=10)\)	0.01	0.15	0.14	0.01

In R, we can generate this matrix, which we call p, by writing:

p <- matrix(c(.02,.04,.01,.06,.15,.15,.02,.20,.14,.10,.10,.01),ncol=4) ## this line creates matrix p
p                                                                      ## this line displays matrix p

##      [,1] [,2] [,3] [,4]
## [1,] 0.02 0.06 0.02 0.10
## [2,] 0.04 0.15 0.20 0.10
## [3,] 0.01 0.15 0.14 0.01

A few things to note here

We assigned the values to a variable called p.
The assignment operator is “<-”. This means that if I want to create an object (like a vector or a matrix) that I want to store in the “global environment” (top right hand side), I will write: “name of object” <- “write what the object is”. So if I wanted to store the variable \(xx=2+5\) I would write “xx <- 2+5”.
The values in p do not automatically display. We can display p by typing its name in the console.

SIDE EXERCISE to understand what c() and matrix() do.

It is important to understand this since data is stored in R as a vector (or matrix). This means that R keeps track of the order that the data is entered in. In particular there is a first element, a second element up to a last element. This means that the vector c(1,2,3) is not the same as c(2,1,3).

For this part you have to submit both the code (the line that displayed the result) and the result that you got.

Write

test1 <- c(1,2,3,4)  ## this line creates the object called "test1"
test1                ## this line displays object "test1"

## [1] 1 2 3 4

Test1 is a row vector with 4 elements. One thing to note is that when you display test1, R gives you an [1] in front of 1,2,3,4. This [1] means that the object you displayed is a vector.

Write

test2 <- matrix(c(1,2,3,4), ncol=1)  ## this line creates the object called "test2"

Then display test2. Test2 is a column vector with 4 rows, right?

Now write

test3 <- matrix(c(1,2,3,4),ncol=2)   ## this line creates the object called "test3"

and display test3. You should now have a 2x2 matrix where the first column is given by (1,2).

Now try

test4 <- matrix(c(1,2,3,4,5,6), ncol=2)  ## this line creates the object called "test4"

and display test4. You should see a 3x2 matrix.

Now try

test5 <- matrix(c(1,2,3,4,5,6), ncol=3)  ## this line creates the object called "test5"

and display test5.

Write the matrices below in R using matrix() and c() notation like above. Call the first matrix, matrix1, and the second, matrix2.

  |1  | 4  | 7 |      and    | 1  | 5 | 
  |---|----|---|             |----|---|
  |2  | 5  | 8 |             | 2  | 6 | 
  |--------|---|             |----|---|
  |3  | 6  | 9 |             | 3  | 7 |
                             |----|---|
                             | 4  | 8 |

Ok. Now we go back to our matrix p.

Matrix p stores all joint probabilities. If you wanted to display P(X=5 and Y=10), which is the element on row 2, column 3 of matrix p, you would have to write

p[2,3]   ## in matrix p, look up element on row 2, column 3

## [1] 0.2

Question 2: Display now P(X=0 and Y=5). What is the code for that? (what is the row and the column for this element stored in matrix p?)

To check that matrix p contains the pmf of \(X\) and \(Y\), you know that if you add all the elements of p you should get \(1\). To see that this is the case, type:

sum(p)  ## sum all elements of matrix p

## [1] 1

sum() is one of the many functions that R comes with. You can type mean(p) or median(p) and see what happens.

Computing marginal probabilities.

Remember that, for example, the marginal probability of \(X=0\) is given by: \[ P(X=0) = P(X=0 \land Y=0) + P(X=0 \land Y=5) + P(X=0 \land Y=10) + P(X=0 \land Y=15)\] This is the sum of all the elements in the first row of matrix p, right? Then \(P(X=5)\) is given by the sum of all elements in the second row of matrix p, and \(P(X=10)\) is given by the sum of all elements in the last row of matrix p. You can compute all three marginal probabilities with the following line:

px <- apply(p,1,sum) ## create marginal probabilities for X  
px                   ## display these marginal probabilities

## [1] 0.20 0.49 0.31

The apply function says: Take matrix p, and, row by row (“1” means “by row”), compute the sum of the elements in the row. This function is very powerful since it computes the sums for all rows at once.

To compute the marginal probabilities for Y, we would have to sum over the columns. To do this, we will use again the apply function but we change \(1\) to \(2\), because \(2\) means “column-by-column”. The code below says: Take matrix p, and column by column, compute the sum of the elements in each column.

py <- apply(p,2,sum) ## create marginal probabilities for Y

Question 3: Display \(P(Y=0)\) and \(P(Y=15)\). Check these numbers by hand. Are the numbers you computed by hand the same as those in the vector py computed with the code above?

Computing conditional probabilities

To compute conditional probabilities, we apply the formula that links conditional probabilities to joint and marginal probabilities. For example, to compute \[P(X=5|Y=5)=\frac{0.15}{0.36} = 0.4166667\] we would have to take the element p[2,2] and divide it by py[2]. The code for this is:

p_x5_y5 <- p[2,2]/py[2]  ## computes conditional probability P(X=5|Y=5)
p_x5_y5                  ## to display this conditional probability

## [1] 0.4166667

Question 4: Compute \(P(X=0|Y=5)\) and \(P(X=10|Y=5)\). Call the first conditional probability p_x0_y5 and the second one p_x10_y5.

We now store the conditional probability \(P(X|Y=5)\) as a row vector with 3 elements and call this object p_x_y5:

p_x_y5 <- c(p_x0_y5,p_x5_y5,p_x10_y5)

Question 5: Display the vector for \(P(X|Y=5)\).

Computing expectations and variances

To compute expectations and variances we need to enter the values of \(X\) and \(Y\):

x<- c(0,5,10)
y<- c(0,5,10,15)

To compute the expectation of \(X\) and of \(X^2\), we write:

EX  <- sum(px*x)   ## expectation of X
EX2 <- sum(px*x^2) ## expectation of X^2

Question 6: Display EX and EX2. What are they equal to? Do you get the same numbers when you compute them by hand?

Question 7: Write the code for \(E(Y)\) and \(E(Y^2)\)

To compute \(E(X|Y=5)\) which is given by \[E(X|Y=5) = 0P(X=0|Y=5) + 5P(X=5|Y=5) + 10P(X=10|Y=5)\] we could write:

EX_Y5 <- sum(p_x_y5*x)

Question 8: What is this conditional expectation?

To compute the variance of \(X\), we write down the formula:

VX <- EX2 - (EX)^2
VX

## [1] 12.4475

Question 9: Write down the code for VY, the variance of \(Y\).

Compute cumulative distribution functions for \(X\) and \(Y\)

Remember the formula for CDFs. For example: \[P(X \leq 5) = P(X \leq 0) + P(X \leq 5) = px[1] + px[2]\] and \[P(X \leq 10) = P(X \leq 0) + P(X \leq 5) + P(X \leq 10) = px[1] + px[2] + px[3]\] where e.g. px[1] is the first element of the vector px, which stands for \(P(X \leq 0)\).

To write a one-liner for the CDF of X we use the command “cumsum()”.

CDF_x <- cumsum(px)  ## returns vector whose elements are the cumulative sums of the elements of px

Question 10: Write down the code for the CDF of \(Y\), and call it CDF_y.

Plot the pmfs and the CDFs