A Finite Math Companion in R Markdown

Foreword

The disciplines of statistics and probability build upon elementary subjects in the field of finite mathematics. This guide compiles many frequently used formulas, visualizations and theorems in this field of study into R programming language and presents these concepts in a written, textbook-style guide. In this guide, each chapter contains a different general subject matter, and contains a variety of formulas that can be called as a function. These functions may be called in a series, or nested into other formulas, to generate useful calculations and visualizations.

Put simply, finite mathematics is a collection of idioms, principles, and theorems concerning generally whole numbers. This includes a variety of subjects “other than calculus”. The subjects of matrix operations, linear programming, and set theory are key aspect of finite maths - and these subjects can act as the basis for many statistical analysis projects. Understanding how to define a sample space and determine expected value underpins all calculations of probability.

Many important financial formulas are also considered germane to the study of finite math, providing formulas to solve problems in everyday finance. Subjects of interest, annuities, sinking funds, amortization, and present value can help answer questions like: How much will I pay in total for a car or home loan?; How much do I need to save to reach a specified retirement goal?

We will also explore some real-world applications of these concepts, such as: determining how the quantity of a set of ingredients needed to have to make a certain quantity of a final recipe; to allocate resources needed to maximize profits, or applying the Leontief input-output model to extend this idea to the entire economy - determining the “true cost” of bringing a quantity of a resource to market, in cases when that resource (like fuel) may require some of itself in the process of producing or transporting of that resource.

This paper is not intended as a comprehensive guide to the field of finite mathematics; rather, it is intended to show step-by-step how the R programming language can be used as a tool to complete finite mathematics problems, and to highlight the capabilities of R for matrix operations, visualization of sets, and calculation of expected value. We will explore how various code constructs can be used to quickly define sample spaces, run Bernoulli trials, and generate randomized and specified elements for set-wise comparison. We will explore R’s ability to analyze sets for intersections, unions, complements and subsets, as well as to test and confirm that the sets conform to mathematical laws and counting formulas.

By carefully considering the problem, and constructing sample spaces where each outcome is equally likely, the expected value can be estimated, and the probability of any given outcome can be calculated using Bayesian statistics and rules about conditional probability under conditions of dependence, independence, or mutual exclusivity.

R and R Markdown provide a powerful platform for general statistical programming. Visualization capabilities in R are enhanced through the “Grammar of Graphics” as implemented through the GGPlot2 and GGForce packages. Data manipulation functions like select, mutate, and summarize are added by the dplyr package, enabling more convenient and readable R language constructs, using the pipe %>% operator. This operator can help your write intuitive syntax in R that resembles common structured query languages like SQL.

This document will not provide basic instruction in R. It is assumed that the reader has some familiarity with the R language, syntax patterns, installation, and configuration. If this is your first introduction to R, please check out the appendix to learn more about the basics of this language and installing the RStudio programming environment. The Rmd files (which this document was written in) may be downloaded from github, enabling the reader to change any portion of this document and customize the word problems to cover different applied subjects.

This guide is intended for finance professionals, statisticians, data scientists, math teachers, students*, or anyone else who needs to generate, visualize, apply, or solve finite math problems.

Note to students: Many of these formulas can be used to solve finite math challenges; but you should avoid using this guide to answer questions on exams or as a shortcut to doing the hard work of understanding this important subject. The only way to the reach the goal of deep understanding, is by doing the work by hand. So please, do yourself a favor, and learn these concepts on pencil and paper before you use this guide. Your instructors will always want to see your work, and this guide won’t help you with that!

Load packages

This Rmarkdown file will use several packages as defined below. This code will load the packages only of they are not already installed. These packages will be loaded into the local scope via the library functions.

if(!require(pracma)) install.packages("pracma", repos = "http://cran.us.r-project.org")
if(!require(dplyr)) install.packages("dplyr", repos = "http://cran.us.r-project.org")
if(!require(ggplot2)) install.packages("ggplot2", repos = "http://cran.us.r-project.org")
if(!require(ggforce)) install.packages("ggforce", repos = "http://cran.us.r-project.org")
if(!require(gtools)) install.packages("gtools", repos = "http://cran.us.r-project.org")
if(!require(matlib)) install.packages("matlib", repos = "http://cran.us.r-project.org")
if(!require(MASS)) install.packages("MASS", repos = "http://cran.us.r-project.org")
if(!require(kableExtra)) install.packages("kableExtra", repos = "http://cran.us.r-project.org")
if(!require(RcppAlgos)) install.packages("RcppAlgos", repos = "http://cran.us.r-project.org")
if(!require(latex2exp)) install.packages("latex2exp", repos = "http://cran.us.r-project.org")

library(pracma)
library(dplyr)
library(ggplot2)
library(ggforce)
library(gtools)
library(matlib)
library(MASS)
library(kableExtra)
library(RcppAlgos)
library(latex2exp)

options(scipen=999)

Linear Equations

This section demonstrates some ways to perform basic visualizations that we will build upon later. Many applied problems in finite math involve drawing lines on a plane. GGplot makes this type of visualization fast, easy and convenient. When you pass data frames into a ggplot object and set the x and y variables, ggplot automatically determines the range of each axis and generates visual output to the graphics device in R. These plots can be rendered to the screen or saved in a variable that you add additional layers to.

Graphing 2 Points on a Plane

Create two points and build a simple plot with just those 2 points at {0,0} and {10,10}. We will make them size 3 to ensure they are easy to see.

x_points <- c(0,10)
y_points <- c(0,10)
coords <- data.frame(x=x_points,y=y_points)

simple_plot <- ggplot(data=coords, aes(x=x, y=y)) + 
  geom_point(size=3) 
simple_plot

Note how GGPlot has automatically made the plot area appropriate to the location of the two points.

Graphing a Line between 2 Points

Create two points and build draw a line with just 2 points at {0,0} and {10,10}

To add a line between these two points, simply add a geom_line object. In the GGPlot syntax - the aes() argument defines where the data is coming from (coords) and which fields in the dataset should correspond to x/y points drawn on the chart.

x_points <- c(0,10)
y_points <- c(0,10)
coords <- data.frame(x=x_points,y=y_points)

simple_plot <- ggplot(data=coords, aes(x=x, y=y)) + 
  geom_point(size=3) + 
  geom_line()
simple_plot

Slope of a Line

Our first equation will be to calculate the slope of a line. For each formal equation we create, we will create a function that encapsulates that function.

Calculate the slope of a line and display it on a chart as an annotation:

#x/y points
x_points <- c(0,8)
y_points <- c(0,6)
coords1 <- data.frame(x=x_points,y=y_points)

#slope of a line equation
slope_line <- function(coords1){
  m <- (coords1$y[2] - coords1$y[1]) / (coords1$x[2] - coords1$x[1])
  m
} 

#run the function and set the output to variable m
m <- slope_line(coords1)

#plot
simple_plot <- ggplot(data=coords1, aes(x=x, y=y)) + 
  geom_point(size=3) + 
  geom_line() + 
  annotate("text", x=3,y=3,label = paste("Slope:",m))
simple_plot

Graphing a Line with geom_smooth

This is a variant of the above line graphs. the geom_smooth produces a linear regression line between n points. This will draw a line, based on your points, using a model of your choosing. We will start with a linear model with only 2 observations - this will produce a straight line betwee n these two points.

This object will be important in later chapters as we explore plots that contain many points. In this section we are introducing the coord_cartesian object. This object will “zoom in” on a portion of the area, without dropping any points. This is similar to the scale object which sets the scale of the area. The difference is that scale removes observations outside the visible scale limits.

Adding the fullrange=TRUE tells GGplot to extend the line into an infinite space, filling the chart area from edge to edge. This is appropriate for linear programming formulas since the line is not finite, and we are interested in knowing where 2 or more lines intersect, without defining a start or end point of the lines. The starting and ending points of the line don’t exist - all that exists is a set of observations. A line should follow the slope of the known observations without terminating at any specified location in an x/y coordinate space.

#x/y points
x_points <- c(2,6)
y_points <- c(2,8)
coords1 <- data.frame(x=x_points,y=y_points)

#set bounds of focus area
upper_bound <- 10
lower_bound <- -10
left_bound <- -10
right_bound <- 10

#set limits to extend beyond the bounds
lower_bound_f <- ifelse(lower_bound>0,0,lower_bound)
left_bound_f <- ifelse(left_bound>0,0,left_bound)
upper_bound_f <- ifelse(upper_bound<0,0,upper_bound)
right_bound_f <- ifelse(right_bound<0,0,right_bound)

#plot
simple_plot <- ggplot() + 
  geom_point(data=coords1, size=3, aes(x=x, y=y)) +
  xlim(left_bound_f*2, right_bound_f*2) +
  ylim(lower_bound_f*2, upper_bound_f*2) +
  geom_smooth(data=coords1, aes(x=x, y=y),method="lm",fullrange=TRUE,color="black") +
  geom_hline(yintercept=0) + 
  geom_vline(xintercept=0) + 
  coord_cartesian(xlim = c(left_bound, right_bound),
                  ylim = c(lower_bound, upper_bound), 
                  expand = FALSE)
simple_plot

Graphing a Linear Equation with 3 or more Points with geom_smooth

This geom_smooth produces a linear regression line between 3 points. The shaded area represents the standard error of this regression line using a 95% confidence interval by default.

#x/y points
x_points <- c(0,3,5,8)
y_points <- c(0,3,3,6)
coords1 <- data.frame(x=x_points,y=y_points)

#plot
simple_plot <- ggplot(data=coords1, aes(x=x, y=y)) + 
  geom_point(size=3) + 
  geom_smooth(method="lm")
simple_plot

The geom_smooth may also be used to draw a line using a loess kernal method to produce a curved regression line between 3 or more points. Drawing a smoothed curve is not strictly considered part of the discipline of finite maths, but it will be important in later chapters on probability and statistics.

#x/y points
x_points <- c(0,3,5,8)
y_points <- c(0,3,3,6)
coords1 <- data.frame(x=x_points,y=y_points)

#plot
simple_plot <- ggplot(data=coords1, aes(x=x, y=y)) + 
  geom_point(size=3) + 
  geom_smooth(method="loess")
simple_plot

Graphing a line with arrows with geom_segment

Another way to represent an infinite line that continues into an endless space, is to use the geom_segment object. This geom_segment object produces a line between 2 points with an arrowhead on the end. You can add two lines that are both covering the same x/y space, but going opposite directions, with points size set to 0, to simulate the appearance of an infinite line as it might appear in a traditional math text. This approach is a bit cumbersome!

x_points <- c(0,10)
y_points <- c(0,10)
coords <- data.frame(x=x_points,y=y_points)

#plot
simple_plot <- ggplot(data=coords, aes(x=x, y=y)) + 
  geom_point(size=0) + 
  geom_segment(x=coords$x[1],
               xend=coords$x[2],
               y=coords$y[1],
               yend=coords$y[2],
               arrow = arrow(length = unit(.2, 'cm'))) + 
  geom_segment(x=coords$x[2],
               xend=coords$x[1],
               y=coords$y[2],
               yend=coords$y[1],
               arrow = arrow(length = unit(.2, 'cm')))
simple_plot

Graphing 2 Lines

Create 2 sets of 2 points and draw a line between them.

In this example, the two equations are being treated as 2 separate data frames. This enables us to pass these into the ggplot objects data parameter for each element we want to draw. We are also adding some text to show the {x,y} coordinates of each point.

Tech Note: The data parameter in this example has moved out of the ggplot() object, and into the individual geom_ elements. This is because this ggplot constructor does not have one “primary” data set, but two. Any settings specified in the ggplot() constructor including data or aes mappings, will be inherited by child objects added to the plot via the + operator. Further, any parameters set in the child object will override any setting in the ggplot constructor. If you are familiar with CSS style sheets - GGplot has similar inheritance structures. You set styles in the GGplot that cascade from a less specific element (GGplot) to a more specific element like geom_line, geom_text, etc.

# points
x_points <- c(0,10)
y_points <- c(0,4)
coords1 <- data.frame(x=x_points,y=y_points)

x_points <- c(0,6)
y_points <- c(8,0)
coords2 <- data.frame(x=x_points,y=y_points)

m1 <- round(slope_line(coords1),2)
m2 <- round(slope_line(coords2),2)

multi_plot <- ggplot() +  
  geom_text(data=coords1, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="black") +
  geom_point(data=coords1, aes(x=x, y=y),color="blue") +
  geom_line(data=coords1, aes(x=x, y=y),color="blue") + 
  annotate("text", x=9,y=3,color="blue",label = paste("Slope:",m1)) +
  geom_text(data=coords2, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="black") +
  geom_point(data=coords2, aes(x=x, y=y),color="darkgreen") +
  geom_line(data=coords2, aes(x=x, y=y),color="darkgreen") + 
  annotate("text", x=2,y=6,color="darkgreen",label = paste("Slope:",m2))
multi_plot

In the chart above, the slope of the blue line is a positive number since x and y axis are increasing together - visually the line is rising from left to right. The slope of the green line is a negative number, since the y axis is decreasing as the x axis is increasing.

Check Intersection of 2 Lines

This function will to determine if 2 lines intersect by checking the slope of each line using our slope_line function on both sets of coordinates. If the slope is not equal - the lines must intersect (eventually). The function will return TRUE if the lines intersect.

Tech Note: This function assumes that the lines continue forever in both directions. It is only checking that the slope of the two lines is not equal, not detecting if any two segments overlap.

#x/y points
x_points <- c(0,10)
y_points <- c(0,4)
coords1 <- data.frame(x=x_points,y=y_points)

x_points <- c(0,6)
y_points <- c(8,0)
coords2 <- data.frame(x=x_points,y=y_points)

#check if the slopes intersect
check_intersect <- function(coords1,coords2){
  m1 <- slope_line(coords1)
  m2 <- slope_line(coords2)
  is_intersecting <- ifelse(m1 == m2,FALSE,TRUE)
  is_intersecting
}

#run intersect check
intersects <- check_intersect(coords1,coords2)

#plot
multi_plot <- ggplot() +  
  geom_point(data=coords1, aes(x=x, y=y),color="blue") +
  geom_line(data=coords1, aes(x=x, y=y),color="blue") +
  geom_point(data=coords2, aes(x=x, y=y),color="green") +
  geom_line(data=coords2, aes(x=x, y=y),color="green") + 
  annotate("text", x=5,y=5,label = paste("Intersecting?:",intersects))
multi_plot

If the lines do not intersect, as shown below, they must be parallel lines.

#x/y points
x_points <- c(0,10)
y_points <- c(10,0)
coords1 <- data.frame(x=x_points,y=y_points)

x_points <- c(0,9)
y_points <- c(9,0)
coords2 <- data.frame(x=x_points,y=y_points)

#run intersect check
intersects <- check_intersect(coords1,coords2)

#plot
multi_plot <- ggplot() +  
  geom_point(data=coords1, aes(x=x, y=y),color="blue") +
  geom_line(data=coords1, aes(x=x, y=y),color="blue") +
  geom_point(data=coords2, aes(x=x, y=y),color="green") +
  geom_line(data=coords2, aes(x=x, y=y),color="green") + 
  annotate("text", x=6,y=6,label = paste("Intersecting?:",intersects))
multi_plot

### Shading an Area above a Line

This chart includes a polygon shading of the area above the line. This could be useful for linear programming and optimization problems. The build_polygon function defines the bounds of a shading polygon.

build_polygon <- function(xr, yr, slope = 1, intercept = 0, above = TRUE){
    #Assumes ggplot default of expand = c(0.05,0)
    x_bool <- xr + 0.05*diff(xr)*c(-1,1)
    y_bool <- yr + 0.05*diff(yr)*c(-1,1)

    #Find where the line crosses the plot edges
    y_sect <- (y_bool - intercept) / slope
    x_sect <- (slope * x_bool) + intercept

    #Build polygon by cases
    if (above & (slope >= 0)){
        area_rows <- data.frame(x=-Inf,y=Inf)
        if (x_sect[1] < y_bool[1]){
            area_rows <- rbind(area_rows,c(-Inf,-Inf),c(y_sect[1],-Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(-Inf,x_sect[1]))
        }
        if (x_sect[2] < y_bool[2]){
            area_rows <- rbind(area_rows,c(Inf,x_sect[2]),c(Inf,Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(y_sect[2],Inf))
        }
    }
    if (!above & (slope >= 0)){
        area_rows <- data.frame(x= Inf,y= -Inf)
        if (x_sect[1] > y_bool[1]){
            area_rows <- rbind(area_rows,c(-Inf,-Inf),c(-Inf,x_sect[1]))
        }
        else{
            area_rows <- rbind(area_rows,c(y_sect[1],-Inf))
        }
        if (x_sect[2] > y_bool[2]){
            area_rows <- rbind(area_rows,c(y_sect[2],Inf),c(Inf,Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(Inf,x_sect[2]))
        }
    }
    if (above & (slope < 0)){
        area_rows <- data.frame(x=Inf,y=Inf)
        if (x_sect[1] < y_bool[2]){
            area_rows <- rbind(area_rows,c(-Inf,Inf),c(-Inf,x_sect[1]))
        }
        else{
            area_rows <- rbind(area_rows,c(y_sect[2],Inf))
        }
        if (x_sect[2] < y_bool[1]){
            area_rows <- rbind(area_rows,c(y_sect[1],-Inf),c(Inf,-Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(Inf,x_sect[2]))
        }
    }
    if (!above & (slope < 0)){
        area_rows <- data.frame(x= -Inf,y= -Inf)
        if (x_sect[1] > y_bool[2]){
            area_rows <- rbind(area_rows,c(-Inf,Inf),c(y_sect[2],Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(-Inf,x_sect[1]))
        }
        if (x_sect[2] > y_bool[1]){
            area_rows <- rbind(area_rows,c(Inf,x_sect[2]),c(Inf,-Inf))
        }
        else{
            area_rows <- rbind(area_rows,c(y_sect[1],-Inf))
        }
    }

    return(area_rows)
}

Now we will apply the build_polygon to a linear graph and shade above a portion of the graph area.

#x/y points
x_points <- c(-5,5)
y_points <- c(4,2)
coords1 <- data.frame(x=x_points,y=y_points)

#set bounds of focus area
upper_bound <- 10
lower_bound <- -5
left_bound <- -5
right_bound <- 5

#set limits to extend beyond the bounds
lower_bound_f <- ifelse(lower_bound>0,0,lower_bound)
left_bound_f <- ifelse(left_bound>0,0,left_bound)
upper_bound_f <- ifelse(upper_bound<0,0,upper_bound)
right_bound_f <- ifelse(right_bound<0,0,right_bound)

sl <- diff(coords1$y) / diff(coords1$x)
int <- coords1$y[1] - (sl*coords1$x[1])

#Build the polygon
polygon_area <- build_polygon(range(coords1$y),range(coords1$x),
            slope=sl,intercept=int,above=FALSE)

min_x <- min(coords1$x)
max_x <- max(coords1$x)

df_poly_above <- coords1 %>%
  add_row(x = c(max_x, min_x),
                  y = c(Inf, Inf))

#plot
simple_plot <- ggplot() + 
  geom_point(data=coords1, aes(x=x, y=y)) +
  xlim(left_bound_f*2, right_bound_f*2) +
  ylim(lower_bound_f*2, upper_bound_f*2) +
  geom_smooth(data=coords1, aes(x=x, y=y),method="lm",fullrange=TRUE) +
  geom_polygon(data = df_poly_above, 
               aes(x = x, y = y),
               fill = "red",
               alpha = 1/5) +
  geom_hline(yintercept=0) + 
  geom_vline(xintercept=0) + 
  coord_cartesian(xlim = c(left_bound, right_bound),
                  ylim = c(lower_bound, upper_bound), 
                  expand = FALSE)
simple_plot

Application: Rating a Ski Slope

In this example, we must decide whether a ski slope should be rated. This requires some basic trigonometry. We will use the slope measurement to determine the average angle of the slope, and apply the Pythagorean theorem to determine the actual distance of the run based on the distance from the top and bottom, when measured on a map. We will also use the inverse tangent (atan) function to convert from the slope, to the angular measurement (in radians), and from that number, convert radians to degrees.

To decide the rating, we will use the rule that any ski slope with less than a 25 degrees average angle is a Green Circle (Easy), 25-40 degrees is a Blue Square (Intermediate), 40 - 50 degrees is Black Diamond (Advanced), and over 50 degrees is rated Double Black Diamond (Expert).

# get rating
get_rating <- function(height_ft, length_ft){
  #get the slope as rise / run
  slope <- height_ft / length_ft
  #convert the slope to radians
  angle_rad <- atan(slope)
  #convert radians to degrees
  angle_deg <- angle_rad * (180/pi)
  #round the angle to 2 decimal places
  angle <- round(angle_deg,2)
 
  #determine the correct rating
  rating <- 'Double Black Diamond'
  rating <- ifelse(angle < 50,'Black Diamond',rating)
  rating <- ifelse(angle < 40,'Blue Square',rating)
  rating <- ifelse(angle < 25,'Green Circle',rating)
  
  #add a color scale for each rating
  color <- 'orange'
  color <- ifelse(angle < 50,'black',color)
  color <- ifelse(angle < 40,'blue',color)
  color <- ifelse(angle < 25,'green',color)
  
  # determine the actual length of the hypotenuse
  # Pythagorean theorem states that A squared + B squared = C squared
  distance <- sqrt(height_ft^2 +  length_ft^2)
  distance <- round(distance)
  
  #return the values object
  return_obj <- list(
    height = height_ft,
    length = length_ft,
    rating = rating,
    color = color,
    angle = angle,
    distance = distance
  )

return_obj
}

# height from top
height_ft <- 5000

# distance from top
length_ft <- 6000

# return details
ski_slope  <- get_rating(height_ft,length_ft)

#set returned variables into a local scope
rating <- ski_slope$rating
distance <- ski_slope$distance
angle <- ski_slope$angle

result <- paste("A straight ski slope that is ", height_ft," feet tall with the peak a length of ", length_ft," feet away from the bottom (when measured from above on a map), is a total of about ",distance," feet long, has an average angle of ",angle," degrees and should be rated: ",rating, sep="")

Result
A straight ski slope that is 5000 feet tall with the peak a length of 6000 feet away from the bottom (when measured from above on a map), is a total of about 7810 feet long, has an average angle of 39.81 degrees and should be rated: Blue Square

Systems of Linear Equations

In this chapter, we consider how one might solve a system of equations in R using matrix methods.

Single Linear Expression

This code demonstrates how you might set up a linear equation. In the example below - you are trying to determine the values of x1 and x2 for each set, such that the eval_linear variable would evaluate to TRUE.

First we will establish the basic logic, then, we will build some code to simulate the steps required to solve a system of linear equations.

#8x + 6y = 3580

#define coefficients
a1 <- 8
a2 <- 6

#define right hand side
b <- 3580

#set variables to 0
x1 <- 0
x2 <- 0

#--------
set1 <- data.frame(a1=a1,a2=a2,b=b,x1=x1,x2=x2)

#a linear equation with 2 variables has the form
eval_linear1 <- (set1$a1 * set1$x1) + (set1$a2 * set1$x2) == set1$b
eval_linear1

## [1] FALSE

Given the code above - its clear our equation is not solved yet - because $8 * 0$ + $6 * 0$ is not 3580. So the challenge here is to figure out what is a correct value for x1 and x2, to make this evaluate to TRUE?

One way to solve this would be to plug in a number for y, then solve for x. Lets make y = 100 and see what we get for x.

#8x + 6y = 3580

#define coefficients
a1 <- 8
a2 <- 6

#define right hand side
b <- 3580

#set variables to 0
x1 <- 0
x2 <- 100

#figure out what is 3580 - 600
c <- b - (a2 * x2)

#divide that number by a1 to get x1
x1 <- c / a1 

#--------
set1 <- data.frame(a1=a1,a2=a2,b=b,x1=x1,x2=x2)

#a linear equation with 2 variables has the form
eval_linear1 <- (set1$a1 * set1$x1) + (set1$a2 * set1$x2) == set1$b
eval_linear1

## [1] TRUE

OK, we have figured out one feasible value for x and y. But let’s face it - this is pretty easy - you could choose any value for Y and get a value for X that satisfies just one linear equation.

result <- paste("If ", a1,"x + ",a2,"y = ",b," then x could be ",x1," and y could be ",x2, sep="")

Result
If 8x + 6y = 3580 then x could be 372.5 and y could be 100

Set of Linear Expressions

Expanding on the section above - this code represents the setup of a set of 2 linear equations each having 2 variables. In this case, we need to find a value for x and y that satisfies not just one, but both of the systems. This is much tougher - since you can’t just pick a random number for y and solve for x - you need your linear equations to evaluate to true for BOTH of the systems using the same X and the same Y.

#8x + 6y = 3580

#define coefficients
a1 <- 8
a2 <- 6

#define right hand side
b <- 3580

#set variables to 0
x1 <- 0
x2 <- 0

set1 <- data.frame(a1=a1,a2=a2,b=b,x1=x1,x2=x2)


#--------

#define coefficients
a1 <- 1
a2 <- 1

#define right hand side
b <- 525

#set variables to 0
x1 <- 0
x2 <- 0

set2 <- data.frame(a1=a1,a2=a2,b=b,x1=x1,x2=x2)

#--------


#a linear equation with 2 variables has the form
eval_linear1 <- (set1$a1 * set1$x1) + (set1$a2 * set1$x2) == set1$b
eval_linear2 <- (set2$a1 * set2$x1) + (set2$a2 * set2$x2) == set2$b

eval_linear1

## [1] FALSE

eval_linear2

## [1] FALSE

Now that we have selected some values - we need to figure out how to get an X and Y that work in both systems. So how will we do this? Picking a random number is not very efficient - we would have to pick a million random numbers and maybe by some dumb luck we would find one that works. But this is where matrix operations come in super handy - Lets take a look.

Matrix Operations

Matrix operations are important in solving systems of equations. The matrix construct in R enables a simple vector to be converted into an $m \times n$ matrix. Standard addition, multiplication, and other arithmetic operations may then be performed on the matrices. R’s natural handling of matrix algebra makes it a great tool for mathematical analysis projects. IN this section we’ll introcude perhaps the most important construct in the R language: the c(). A c() construct how R defines and instantiates a generic vector, in other words, a list of values. A c() can contain as many values as you need to define, and acts as the basis of many additional constructs that can define a list based on rules, ranges, or randomization of elements. WHats important about the c() construct, is that it’s so basic. All other constructs like a matrix, ultimately “feed off of” the c().

#define a list with 9 elements
list_elements <- c(1, 2, 3, 1, 3, 2, 3, 2, 1)

# generate 2 3x3 matrix objects A and B using the same list of elements
A <- matrix(list_elements, 3, 3, byrow = TRUE)
B <- matrix(list_elements, 3, 3, byrow = TRUE)

Lets take a look at the matrix A and matrix B - since they both have the same elements, and they were both created to have the same dimensions, and they both were created by row, they are both - you guessed it - exactly the same.

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    3    2
## [3,]    3    2    1

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    3    2
## [3,]    3    2    1

Matrix Equality

You can check if all of the elements in two matrices are the same by using the double-equals sign equality operator.

# check matrix equality
A_equals_B <- A == B
A_equals_B

##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

You may also use the identical function to determine if two matrices are the same across all elements.

# check matrix equality
A_ident_B <- identical(A,B)
A_ident_B

## [1] TRUE

Now, taking the same list elements we already defined above - lets try one minor alteration - in this example we will construct matrix B “by column” instead of “by row”

# generate 2 3x3 matrix objects A and B using the same list of elements
A <- matrix(list_elements, 3, 3, byrow = TRUE)
B <- matrix(list_elements, 3, 3, byrow = FALSE)

OK, now looking at the two matrices, they are no longer the same.

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    3    2
## [3,]    3    2    1

##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2    3    2
## [3,]    3    2    1

Instead the second matrix was created using the same list of elements but running down column 1, then column 2, then column 3. So now of we try the identical operator - we can see they are not identical - they were built from the same list of values, but they are not the same.

# check matrix equality
A_ident_B <- identical(A,B)
A_ident_B

## [1] FALSE

Now that we got these basics out of the way -lets consider how we can work with our matrices and do various mathematical operations on them.

Arithmetic Operations on Matrices

Arithmetic operations - addition, subtraction, multiplication, and division, can be performed on matrices. But unlike typical arithmetic- matrix arithmetic has different rules, and is not always performed in the way you might expect. It is possible to perform element by element (“element-wise”) arithmetic - which is the “default” in R using standard operators like + , - , * and / but these are not always the correct operators to use for mathematical analysis of matrices.

In this section we will explore the various ways that matrices can be operated upon using booth the standard arithmetic operators, and operators that are specialized to these constructs. IN addition, When performing arithmetic operations on a matrix, one must be very careful to consider several factors before determining what operators to use, as these factor can affect your choices regarding which are the right tools for the job.

WHen thinking about doing arithmetic on a matrix - you must first consider: * what is the type of arithmetic operation * what are the dimensions of the matrices (rows x columns) * what is the order of the matrices (which one is on the left and which one is on the right, etc)

Element-Wise Matrix Addition

Two same-dimension matrices can be added together by simply using the addition operator. This is equivalent to adding each element in matrix A, to the same element in matrix B.

This is called element wise, because each element is performed upon, by running the chosen operations against the same element in the other matrix. This is an important distinction for matrices, because it is actually NOT the “typical” way you would do some of the operations.

# get the additive matrix sums
A_plus_B <- A + B
A_plus_B

##      [,1] [,2] [,3]
## [1,]    2    3    6
## [2,]    3    6    4
## [3,]    6    4    2

Element-Wise Matrix Subtraction

Two matrices can be subtracted by simply using the subtraction operator. This is equivalent to subtracting each element in matrix B, from the same element in matrix A.

# get the additive matrix sums
A_minus_B <- A - B
A_minus_B

##      [,1] [,2] [,3]
## [1,]    0    1    0
## [2,]   -1    0    0
## [3,]    0    0    0

Element-Wise Matrix Multiplication

Two matrices can be multiplied together by simply using the standard multiplication operator. This is equivalent to multiplying each element in matrix A, to the same element in matrix B.

Its important to understand the difference between “element-wise” Matrix multiplication, and traditional Matrix Multiplication. So we will start with the way you might expect a matrix to be multiplied - element by element.

# get the multiplicative matrix product
A_times_B <- A * B
A_times_B

##      [,1] [,2] [,3]
## [1,]    1    2    9
## [2,]    2    9    4
## [3,]    9    4    1

Element-Wise Matrix Division

Two matrices can be divided by simply using the division operator. This is equivalent to dividing each element in matrix A, by the same element in matrix B.

# get the multiplicative matrix product
A_dividedby_B <- A / B
A_dividedby_B

##      [,1] [,2] [,3]
## [1,]  1.0    2    1
## [2,]  0.5    1    1
## [3,]  1.0    1    1

Scalar Matrix Multiplication

To multiply a matrix by a scalar - a numeric constant N - you simply multiply each element by the same scalar. This will produce a matrix with the same number of elements. Each element in your matrix is multiplied by the same scalar, resulting in a matrix where every element is N times your original matrix.

scalar <- 2
A_times_scalar <- A * scalar
A_times_scalar

##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    2    6    4
## [3,]    6    4    2

The Dot Product

So what is the dot product? Hint: it’s not a stain remover for when your pocket protector fails!

The dot product is a reference to a method of multiplication when you have two lists of the same length. The dot product is a result of multiplication of these two lists together to result in a final, singular value. To get the dot product you simply multiply two lists, element by element, and add each of these products together.

$(a1 * b1) + (a2 * b2) + ... + (aN * bN)$

The dot product is perhaps the simplest form demonstrating the concept of an inner product in linear algebra, this is a concept which we will cover in greater depth, in a later installment of this guide.

Matrix Multiplication

Matrix multiplication, is when you calculate the dot product that results from multiplying each row in matrix A by each column in matrix B. This requires that the two matrices are conformable.

To multiply 2 conformable matrices in R, you must use the %*% operator. This is the only way to perform correct matrix multiplication- which differs from element-wise matrix multiplication. Matrix multiplication multiplies each row by each column in a this specific way.

This results in an output matrix where each element is the sum of the multiplied elements from each row of matrix A (from left to right) multiplied by the element in the same column of matrix B (from top to bottom).

As an example, you would take row 1 column 1 of A, times column 1 row 1 of B plus, row 1 column 2 of A times column 1 row 2 of B, plus row 1 column 3 of A, times column 1 row 3 of B, etc. All these values, added together would result in a single numeric outcome, and be placed into row 1 column 1 of your output matrix.

To put it mathematically -

$Cij = (Ai_1*B_1j)+(Ai_2*B_2j) + ... + (Ain * Bnj)$

By considering each pair of rows in a and columns in B as a list, and getting the dot product, you will find the output of that element in your output matrix. It sounds confusing at first, but after you have done a few of these operations, it all begins to feel very natural.

Lets reload our identical matrices:

# generate 2 3x3 matrix objects A and B using the same list of elements
A <- matrix(list_elements, 3, 3, byrow = TRUE)
B <- matrix(list_elements, 3, 3, byrow = TRUE)

Recall that the elements of these 2 matrices are the same, in the same order. But the rows of A do not equal the columns of B - taken this way, the rows and columns of these two identical matrices, are actually quite different.

A_times_B <- A %*% B
A_times_B

##      [,1] [,2] [,3]
## [1,]   12   14   10
## [2,]   10   15   11
## [3,]    8   14   14

Individual elements, or entire rows or columns of a matrix can be retrieved using the base R square brace syntax.

This syntax returns an element or vector from your matrix in the position specified by [row,column]

If both the row and column are specified, it returns one element, at the specified row and column.

If the column side is left blank, it takes all elements in all the columns in that row.

If the row side is left blank, it takes all elements in all the rows in that column.

In our matrix multiplication operation above, we have taken row 1 of A - just like using the syntax:

row_1_A <- A[1,]
row_1_A

## [1] 1 2 3

And column 1 of B using the syntax:

col_1_B <- B[,1]
col_1_B

## [1] 1 1 3

And we have calculated the dot product of these two lists:

$(1*1) + (2*1) + (3*3) = 12$

The %*% performs this same operations for each combination of rows of A and columns of B.

Matrices with Unequal Dimensions

In matrix operations, you can take two matrices with unequal dimensions, and run matrix algebra on them. This is only possible given if $m \times n$ matrix A is multiplied by an $n \times p$ matrix B.

In other words - the number of columns in the first matrix must equal the number of rows in the second matrix.

This is one way that matrix algebra differs from other forms of algebra - because ORDER MATTERS. Unlike the commutative property in traditional algebra: A * B does not always equal B * A.

# generate a 3x1 matrix objects A
A <- matrix(c(1, 2, 3), 3, 1, byrow = TRUE)

# generate a 1x3 matrix objects B
B <- matrix(c(1, 2, 3), 1, 3, byrow = TRUE)

A

##      [,1]
## [1,]    1
## [2,]    2
## [3,]    3

##      [,1] [,2] [,3]
## [1,]    1    2    3

ONce again we will use the the %*% operator to multiply our two matrices. In this case, we are creating a matrix A which is 3 x 1 and a matrix B which is 1 x 3. These are conformable because matrix the number of columns in A equals the number of rows in B.

A_times_B <- A %*% B
A_times_B

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    4    6
## [3,]    3    6    9

Using the matrix multiplication operator we can see that $A \times B$ produces a 3x3 matrix with 9 individual entries.

B_times_A <- B %*% A
B_times_A

##      [,1]
## [1,]   14

Using the same matrix multiplication operator we can see that $B \times A$ produces a 1x1 matrix with just one value.

dim(A_times_B)

## [1] 3 3

dim(B_times_A)

## [1] 1 1

Since the dimensions are not the same, its clear that $A \times B$ and $B \times A$ are NOT the same, and by consequence we have proven that multiplication of matrices is most definitely not commutative.

In the example above our A and B are transpose - the rows in A equal the columns in B and the columns in A equal the rows in B. However you can just as easily multiply two matrices that are not transpose - as long as they follow the rule that the columns in A equal the rows in B.

# generate a 3x1 matrix objects A
A <- matrix(c(1, 2, 3, 4, 5, 6), 3, 2, byrow = TRUE)

# generate a 1x3 matrix objects B
B <- matrix(c(1, 2, 3, 4), 2, 2, byrow = TRUE)

A

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
## [3,]    5    6

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4

A_times_B <- A %*% B
A_times_B

##      [,1] [,2]
## [1,]    7   10
## [2,]   15   22
## [3,]   23   34

Using the matrix multiplication operator we can see that $A \times B$ produces a 3x2 matrix with 6 individual entries.

If we were to try to multiply $B \times A$ - we would receive an Error noting that these are “Non-Conformable Arrays”.

An easy way to verify if your two matrices are conformable is using the nrow and ncol function:

is_conformable <- function(left,right){
  columns_on_left <- ncol(left)
  rows_on_right <- nrow(right)
  is_conformable <- columns_on_left == rows_on_right
}

Now we will define two matrices, one that is 3x2 and one that is 2x4.

# generate a 3x1 matrix objects A
A <- matrix(c(1, 2, 3, 4, 5, 6), 3, 2, byrow = TRUE)

# generate a 1x3 matrix objects B
B <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8), 2, 4, byrow = TRUE)

Our function above checks if A and B are conformable again, considering that order matters and A always comes before (to the left) of B.

can_A_mult_B <- is_conformable(A,B)
can_A_mult_B

## [1] TRUE

Switching the order we pass in a and B:

can_B_mult_A <- is_conformable(B,A)
can_B_mult_A

## [1] FALSE

So once again, we have seen that order matters - A on the left can be multiplied by B on the right, but not the other way around.

Inverse of a Matrix

The inverse of a matrix can be found by adding the identity matrix to the right hand side of a matrix, and performing row reduction until the identity matrix appears on the left hand side. The values on the right hand side, are the inverse of the original matrix.

# generate a 3x1 matrix objects A
A <- matrix(c(-4,7,1,0,1,2,0,1), 2, 4, byrow = TRUE)
A

##      [,1] [,2] [,3] [,4]
## [1,]   -4    7    1    0
## [2,]    1    2    0    1

#reduced a to its rediuced row echelon form
A_ref <- rref(A)
A_ref

##      [,1] [,2]        [,3]      [,4]
## [1,]    1    0 -0.13333333 0.4666667
## [2,]    0    1  0.06666667 0.2666667

#take the 2 columns on right hand side of A_ref to get the inverse of A
A_1 <- A_ref[,3:4]
A_1

##             [,1]      [,2]
## [1,] -0.13333333 0.4666667
## [2,]  0.06666667 0.2666667

SInce the native form of the returned data is in decimal form, we will use the fractions function in MASS package to convert these to vulgar fractions instead.

fractions(A_1)

##      [,1]  [,2] 
## [1,] -2/15  7/15
## [2,]  1/15  4/15

starting with an original 2x2 matrix, with identity matrix appended,

\[ \left[ \begin{array}{cc|cc} -4 & 7 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ \end{array} \right] \]

The reduced row echelon form of the inverse matrix will look like this:

\[ \left[ \begin{array}{cc|cc} 1 & 0 & -2/15 & 7/15 \\ 0 & 1 & 1/15 & 4/15 \\ \end{array} \right] \]

Application: Cryptography

One popular use for matrix multiplication is in the field of cryptography. A message can be encoded by using a lexicon which associated leters with numbers, and then applying matrix mlultipilication to the message in groups of N.

letters <- c('_','A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z')
lexicon <- c(0,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8,9,-9,10,-10,11,-11,12,-12,13,-13)

lex_df <- data.frame(let=letters,lex=lexicon)
lex_df

We now have a mapping between actual letters and numbers. This could follow any pattern, or no pattern at all; as long as the lexicon and key is known to the receiving party, the message can be decoded.

Now we will write out a message and convert it into a series of matrices.

message <- c('A','T','T','A','C','K','_','A','T','_','F','I','R','S','T','_','L','I','G','H','T')
list_triplets <- list()

count_letters <- 0
count_triplets <- 0

triplet <- c()
for(x in message){
  count_letters <- count_letters + 1
  triplet <- c(triplet,x)
  
  if(count_letters %% 3 == 0){
    count_triplets <- count_triplets + 1
    list_triplets[[count_triplets]] <- triplet
    triplet <- c()
  }
  
}

list_triplets

## [[1]]
## [1] "A" "T" "T"
## 
## [[2]]
## [1] "A" "C" "K"
## 
## [[3]]
## [1] "_" "A" "T"
## 
## [[4]]
## [1] "_" "F" "I"
## 
## [[5]]
## [1] "R" "S" "T"
## 
## [[6]]
## [1] "_" "L" "I"
## 
## [[7]]
## [1] "G" "H" "T"

We now have a list of N triplets. Each of our triplets must now be multiplied by a 3x3 encryption matrix - this matrix is our cryptographic key. It can be any number - but must be invertible.

enc_mat <- matrix(c(2,-1,1,1,4,2,4,9,1), 3, 3, byrow = TRUE)
enc_mat

##      [,1] [,2] [,3]
## [1,]    2   -1    1
## [2,]    1    4    2
## [3,]    4    9    1

We can determine if the matrix is invertible using the det function in matlab. If det != 0 - the matrix is invertible.

#is invertible helper function
is_invertible <- function(mat){
  invertible <- FALSE
  if(det(enc_mat)!=0){
    invertible <- TRUE
  }
  invertible
}
#is the matrix invertible?
invertible <- is_invertible(enc_mat)
invertible

## [1] TRUE

We have verified the matrix can be inverted - Now we must multiply each 3-letter triplet by finding the associated numbers in our lexicon for each letter in the matrix

list_outputs <- list()
count_outputs <- 0

for(t in list_triplets){
  #break the triplet into its 3 parts
  letter_1 <- t[1]
  letter_2 <- t[2]
  letter_3 <- t[3]
  
  #pull the associated numeric representation of the letter from our lexicon
  lex_1 <- lex_df %>% filter(let==letter_1) %>% pull(lex)
  lex_2 <- lex_df %>% filter(let==letter_2) %>% pull(lex)
  lex_3 <- lex_df %>% filter(let==letter_3) %>% pull(lex)
  
  #convert the triplet to a matrix
  mult_mat <- matrix(c(lex_1,lex_2,lex_3), 3, 1, byrow = TRUE)
  
  #multiply the triplet by our crytpo key
  out_mat <- enc_mat %*% mult_mat
  
  #add the matrix to a list of matrices
  count_outputs <- count_outputs + 1
  list_outputs[[count_outputs]] <- out_mat
}

#the list_outputs we are sending is now just a bunch of numbers that have no relation to our lexicon.
#this is only useful if you have the original crypto key, and the lexicon.
list_outputs

## [[1]]
##      [,1]
## [1,]    2
## [2,]  -59
## [3,]  -96
## 
## [[2]]
##      [,1]
## [1,]    6
## [2,]   21
## [3,]   28
## 
## [[3]]
##      [,1]
## [1,]  -11
## [2,]  -16
## [3,]   -1
## 
## [[4]]
##      [,1]
## [1,]    8
## [2,]   -2
## [3,]  -22
## 
## [[5]]
##      [,1]
## [1,]  -38
## [2,]   11
## [3,]   44
## 
## [[6]]
##      [,1]
## [1,]   11
## [2,]  -14
## [3,]  -49
## 
## [[7]]
##      [,1]
## [1,]    2
## [2,]  -32
## [3,]  -30

We now have a fully encrypted message, ready to send to our clandestine partner. Since our partner knows the lexicon we developed and the cryptographic key (a 3x3 matrix) - she can decode the message and learn the big secret. To do so, she must simply do the same thing we did - but essentially in reverse.

For this to work, we must invert the matrix. The matlab package contains a function to get the inverse of the matrix. Assuming our partner has the enc_mat matrix she can determine the inverse very easily.

inv_mat <- inv(enc_mat)
inv_mat

##            [,1]        [,2]        [,3]
## [1,]  0.3333333 -0.23809524  0.14285714
## [2,] -0.1666667  0.04761905  0.07142857
## [3,]  0.1666667  0.52380952 -0.21428571

decrypted_message <- ''

for(m in list_outputs){
  #get the decrypted matrix by multiplying the inverse matrix on the left by the triplet
  #rounding is necessary since the inv function returns a decimal that will be off by a few millionths
  decrypted_mat <-  round(inv_mat %*%  m,0)
  
  #assign each number pulled out to a local variable
  lex_1 <- decrypted_mat[1]
  lex_2 <- decrypted_mat[2]
  lex_3 <- decrypted_mat[3]
  
  #lookup the letter in the lexicon matching this number
  lex_1 <- lex_df %>% filter(lex==lex_1) %>% pull(let)
  lex_2 <- lex_df %>% filter(lex==lex_2) %>% pull(let)
  lex_3 <- lex_df %>% filter(lex==lex_3) %>% pull(let)
  
   #paste the entries to the end of a string
  decrypted_message <- paste(decrypted_message,lex_1,lex_2,lex_3,sep='')
}
#show the string
decrypted_message

## [1] "ATTACK_AT_FIRST_LIGHT"

This method of encryption demonstrates one of the many valuable uses for matrices and their inverses.

Gaussian Elimination

Gaussian elimination is a method by which a successive series of row operations may be performed on a matrix to enable the matrix to appear in row reduced echelon form abbreviated to RREF. A consistent system of equations with a single solution, when fully reduced, will take the form of the identity matrix $I$

Valid row operations are as follows:

Swap any 2 rows
Replace a row by any nonzero multiple of that row
Replace a row by the sum of that row, and a nonzero multiple of another row

Gaussian elimination may be performed step-by-step in R using a method similar to the examples below. This is not a very efficient way to perform this series of tasks, but it does show each step very clearly for demonstration and instruction purposes.

# generate 2 3x3 matrix objects A and B
A <- matrix(c(1, 2, 3, 2, 3, 2, 3, 2, 1), 3, 3, byrow = TRUE)
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    3    2
## [3,]    3    2    1

#perform a row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R2 <- (r1 * -2) + r2
A[2,] <- R2
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0   -1   -4
## [3,]    3    2    1

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R3 <- (r1 * -3) + r3
A[3,] <- R3
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0   -1   -4
## [3,]    0   -4   -8

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R3 <- (r2 * -4) + r3
A[3,] <- R3
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0   -1   -4
## [3,]    0    0    8

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R2 <- (r2 * (-1))
A[2,] <- R2
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0    1    4
## [3,]    0    0    8

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R3 <- (r3 * (1/8))
A[3,] <- R3
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0    1    4
## [3,]    0    0    1

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R2 <- (r3 * (-4)) + r2
A[2,] <- R2
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    0    1    0
## [3,]    0    0    1

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R1 <- (r3 * (-3)) + r1
A[1,] <- R1
A

##      [,1] [,2] [,3]
## [1,]    1    2    0
## [2,]    0    1    0
## [3,]    0    0    1

#perform another row operation
r1 <- A[1,]
r2 <- A[2,]
r3 <- A[3,]

R1 <- (r2 * (-2)) + r1
A[1,] <- R1
A

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

The above series of operations may seem rather tedious - why go through all of these steps to just produce an identity matrix? The main use of this method, is to solve a set of linear equations for 2 or more variables; to determine if the set of equations is consistent; and to assess whether there is a single solution, no solution, or many possible solutions.

There are R packages available that can provide most matrix operation features. For example, in the pracma package, the rref function can be used to to find the row reduced echelon form (RREF) of a matrix using Gaussian elimination. An example is shown below, which was taken from the documentation for the pracma package.

#FROM: https://www.rdocumentation.org/packages/pracma

# generate a 3x3 matrix and solve
A <- matrix(c(1, 2, 3, 1, 3, 2, 3, 2, 1), 3, 3, byrow = TRUE)
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    3    2
## [3,]    3    2    1

rref(A)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

# generate a 3x4 matrix and solve
A <- matrix(data=c(1, 2, 3, 2, 5, 9, 5, 7, 8,20, 100, 200),
            nrow=3, ncol=4, byrow=FALSE)  
A

##      [,1] [,2] [,3] [,4]
## [1,]    1    2    5   20
## [2,]    2    5    7  100
## [3,]    3    9    8  200

rref(A)

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0  120
## [2,]    0    1    0    0
## [3,]    0    0    1  -20

Tip:

The default matrix notation in LaTeX uses the bmatrix:

$\begin{bmatrix}1 & 2 & 3\\4 & 5 & 6\end{bmatrix}$

However the above function does not appear to support the vertical bar. To represent a matrix with a right hand side separated by a vertical bar in LaTeX use the array function surrounded by left and right brackets:

\[ \left[ \begin{array}{ccc|c} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 3 \\ 0 & 0 & 1 & -3 \\ \end{array} \right] \]

Linear Programming with 2 Variables

An application of the plots above to display and solve an optimization problem.

Application: Donut Shop

IN this application exercise, We will use a donut shop as an example. A baker has 2 award winning donut recipes that both use chocolate and rainbow sprinkles. The baker wants to know how much of each donut to bake to maximize profits. An optimization problem can answer this question.

Donuts

#name the object you will be studying
subject <- "Donut"
madeby <- "Bakery"

#name your recipes/mixtures x and y 
recipe_x <- "Chocolate Dream Donut"
recipe_y <- "Birthday Surprise Donut"

#define x and y for plots
x_def <- paste("Number of",recipe_x)
y_def <- paste("Number of",recipe_y)

#two ingredients available to mix in a recipe
ing1_def <- "Chocolate Sprinkles"
ing2_def <- "Rainbow Sprinkles"

#the maximum amount of each ingredient available using same units as below
max_ing1_g <- 3600
max_ing2_g <- 2400

#number of items produced per batch
item_count <- 100

#the amount of ingredients 1 and 2 needed by recipe x using same units as above
requires_ing1_x <- 120
requires_ing2_x <- 40

#the amount of ingredients 1 and 2 needed by recipe y using same units as above
requires_ing1_y <- 60
requires_ing2_y <- 120

#the profit per item on recipe x
recipe_x_profit <- .22

#the profit per item on recipe y
recipe_y_profit <- .18

result <- paste("A ",madeby," has ", max_ing1_g," grams of ", ing1_def," and ", max_ing2_g," grams of ", ing2_def," available per day. Each batch of the ", recipe_x," recipe takes ", requires_ing1_x," grams of ", ing1_def," and ", requires_ing2_x," grams of ", ing2_def,". Each batch of the ", recipe_y," recipe  takes ", requires_ing1_y," grams of ", ing1_def," and ", requires_ing2_y," grams of ", ing2_def,".  The ", recipe_x," generates a profit of ", recipe_x_profit," per ",subject,", while the ", recipe_y," generates ", recipe_y_profit," profit per ", subject,". Each batch produces ", item_count," ",subject,"s. How can the company maximize profit assuming all items will sell out? What will the total profit be?", sep="")

#result

Challenge
A Bakery has 3600 grams of Chocolate Sprinkles and 2400 grams of Rainbow Sprinkles available per day. Each batch of the Chocolate Dream Donut recipe takes 120 grams of Chocolate Sprinkles and 40 grams of Rainbow Sprinkles. Each batch of the Birthday Surprise Donut recipe takes 60 grams of Chocolate Sprinkles and 120 grams of Rainbow Sprinkles. The Chocolate Dream Donut generates a profit of 0.22 per Donut, while the Birthday Surprise Donut generates 0.18 profit per Donut. Each batch produces 100 Donuts. How can the company maximize profit assuming all items will sell out? What will the total profit be?

Challenge

A Bakery has 3600 grams of Chocolate Sprinkles and 2400 grams of Rainbow Sprinkles available per day. Each batch of the Chocolate Dream Donut recipe takes 120 grams of Chocolate Sprinkles and 40 grams of Rainbow Sprinkles. Each batch of the Birthday Surprise Donut recipe takes 60 grams of Chocolate Sprinkles and 120 grams of Rainbow Sprinkles. The Chocolate Dream Donut generates a profit of 0.22 per Donut, while the Birthday Surprise Donut generates 0.18 profit per Donut. Each batch produces 100 Donuts. How can the company maximize profit assuming all items will sell out? What will the total profit be?

In the example above, we have set up the variables in the code block, and then placed them inline into a word problem format. This allows the RMarkdown file to create meaningful interpretations of data in the report body.

The challenge here is to frame the problem and solve it using R and RMarkdown.

First we will define the objective to be maximized: P (Profit)

In this linear programming problem with 2 variables,

$P = 0.22 * X + 0.18 * Y$

Subject to the constraints:

120 * X + 60 * Y <= 3600

40 * X + 120 * Y <= 2400

$X >= 0$

$Y >= 0$

First, let’s see how we can visualize the problem using the charting strategies demonstrated above.

The shaded area plots are useful for this type of problem. We need to see where the set of feasible points is.

Using the above code we can start by drawing the first of the constraints on to the chart.

If $X = 0$ then $Y <= 60$

If $Y = 0$ then $X <= 30$

coords0 <- data.frame(x=0,y=0)

#setup the first set of coordinates
first_x <- max_ing1_g / requires_ing1_y
first_y <- max_ing1_g / requires_ing1_x

#x/y points
x_points <- c(0,first_y)
y_points <- c(first_x,0)
coords1 <- data.frame(x=x_points,y=y_points)

paste("coords1:",coords1)

## [1] "coords1: c(0, 30)" "coords1: c(60, 0)"

#setup the second set of coordinates
second_x <- max_ing2_g / requires_ing2_y
second_y <- max_ing2_g / requires_ing2_x

#x/y points
x_points <- c(0,second_y)
y_points <- c(second_x,0)
coords2 <- data.frame(x=x_points,y=y_points)

paste("coords2:",coords2)

## [1] "coords2: c(0, 60)" "coords2: c(20, 0)"

# generate a 2x3 matrix and solve for the intersection point
A <- matrix(c(requires_ing1_x, requires_ing1_y, max_ing1_g, 
              requires_ing2_x, requires_ing2_y, max_ing2_g), 2, 3, byrow = TRUE)
A

##      [,1] [,2] [,3]
## [1,]  120   60 3600
## [2,]   40  120 2400

#use pracma for the rref
R <- rref(A)   
R

##      [,1] [,2] [,3]
## [1,]    1    0   24
## [2,]    0    1   12

#assign to data frame intersectionpoint
intersectionpoint <- data.frame(x=R[1,3],y=R[2,3])

#set bounds of focus area
upper_bound <- 100
lower_bound <- -50
left_bound <- -50
right_bound <- 100

#list & arrange feasible points
points <- data.frame()
points <- points %>% bind_rows(coords0)
points <- points %>% bind_rows(coords1[2,])
points <- points %>% bind_rows(coords2[1,])
points <- points %>% bind_rows(intersectionpoint)
points <- points %>% arrange(x,y)

#plot
simple_plot <- ggplot() + 
  geom_point(data=coords0, aes(x=x, y=y)) +
  geom_point(data=coords1, aes(x=x, y=y)) +
  geom_point(data=coords2, aes(x=x, y=y)) +
  geom_point(data=intersectionpoint, aes(x=x, y=y),size=2,color="red") +
  geom_text(data=intersectionpoint, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="red") +
  geom_text(data=coords0, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="purple") +
  geom_text(data=coords1, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="darkgreen") +
  geom_text(data=coords2, aes(x=x, y=y, label=paste("{",x,",",y,"}")),color="blue") +
  geom_line(data=coords1, aes(x=x, y=y),color='darkgreen') +
  geom_line(data=coords2, aes(x=x, y=y),color='blue') +
  geom_polygon(data=points, aes(x=x, y=y),fill='black',alpha=.2) +
  geom_hline(yintercept=0) + 
  geom_vline(xintercept=0) + 
  coord_cartesian(xlim = c(left_bound, right_bound),
                  ylim = c(lower_bound, upper_bound), 
                  expand = FALSE) +
  labs(x=x_def,y=y_def)
simple_plot

#select the feasible points. This process may require additional 
#analysis of the plot to determine which points are feasible. 
#0,0 has been omitted, and we are rounding to the nearest whole number.
#in practice it may be necessary to use floor rounding (always round down) to avoid over-utilization

#get first feasible point
location0x <- round(0,0)  
location0y <- round(0,0)
outcome0 <- (location0x * recipe_x_profit + location0y * recipe_y_profit) * item_count
row0 <- data.frame(x=location0x,y=location0y,total=outcome0)

#get first feasible point
location1x <- round(coords1[2,]$x,0)  
location1y <- round(coords1[2,]$y,0)
outcome1 <- (location1x * recipe_x_profit + location1y * recipe_y_profit) * item_count
row1 <- data.frame(x=location1x,y=location1y,total=outcome1)

#get second feasible point
location2x <- round(coords2[1,]$x,0) 
location2y <- round(coords2[1,]$y,0)
outcome2 <- (location2x * recipe_x_profit + location2y * recipe_y_profit) * item_count
row2 <- data.frame(x=location2x,y=location2y,total=outcome2)

#get third feasible point
location3x <- round(intersectionpoint[1,]$x,0)
location3y <- round(intersectionpoint[1,]$y,0)
outcome3 <- (location3x * recipe_x_profit + location3y * recipe_y_profit) * item_count
row3 <- data.frame(x=location3x,y=location3y,total=outcome3)

#build a the table of outcomes
outcomes <- data.frame()
outcomes <- outcomes %>% bind_rows(row0)
outcomes <- outcomes %>% bind_rows(row1)
outcomes <- outcomes %>% bind_rows(row2)
outcomes <- outcomes %>% bind_rows(row3)

#return the table of outcomes ordered by the total
outcomes %>% arrange(-total)

#assign the best outcome to bestoutcome
bestoutcome <- outcomes %>% arrange(-total) %>% slice(1)

#set some variables
bestx <- bestoutcome$x
besty <- bestoutcome$y
besttotal<- bestoutcome$total

#show the result.
result <- paste("To maximize profits, you should make ",bestx," batches of ",recipe_x," and ",besty," batches of ",recipe_y," for a total profit of $",besttotal,sep="")

Result
To maximize profits, you should make 24 batches of Chocolate Dream Donut and 12 batches of Birthday Surprise Donut for a total profit of $744

We have completed the linear programming model and determined the total number of each batch that should be produced to maximize profit given the constraints provided, and plotted the key details on a GGPlot.

Finance

In this section we will set up a few equations in R that will perform some of the standard financial formula operations, similar to the linear programming problem we worked out above.

This section contains several key formulas related to interest, annuities, amortization, and sinking funds.

Simple Interest

The simple interest formula is as it sounds - simple!

Simple Interest Formula

$I = Prt$

Determining the amount owed is a fairly simple addition to this formula: $A = P + I$

#create simple interest function
simpleinterest <- function (P,r,t){
  r <- r/100
  I <- P * r * t
  I
}

#amount of principal
P <- 2000

#% rate of loan
r <- 5

#number of years
t <- 3

#interest (unknown)
I <- 0

#run the function
I <- simpleinterest(P,r,t)

# A = Principal plus Interest
A <- P + I

#write a sentence to summarize the result
result <- paste("The simple interest on a loan of $",P," at an annual interest rate of ",r,"%, after ",t," years, is: ", I, " and the total amount due is: $", A, sep="") 

#result

Result
The simple interest on a loan of $2000 at an annual interest rate of 5%, after 3 years, is: 300 and the total amount due is: $2300

Discounted Loans

Let $r$ be the per annum rate of interest, $t$ the time in years and $L$ the amount of the loan. The proceeds R is is given by:

$R=L-Lrt=L(1-rt)$

#create discounted loan function
discountedloan <- function (L,r,t){
  r <- r/100
  R <- L * (1- (r * t))
  R
}

#amount of loan
L <- 10000

#% rate of loan
r <- 8

#number of years
t <- .25

#run the function
R <- discountedloan(L,r,t)

#write out the result 
result <- paste("The proceeds of a discounted loan of $",L," at an annual interest rate of ",r,"%, after ",t," years is: $", round(R,2),sep="") 

#result

Result
The proceeds of a discounted loan of $10000 at an annual interest rate of 8%, after 0.25 years is: $9800

Amount of a Discounted Loan

The amount needed to repay a discounted loan can be found using:

$(\frac{L}{R})A$

#nested function
discountedloanamount <- function (L,r,t){
  R <-discountedloan(L,r,t)
  A <- L / R
  A <- L * A
  A
}

#run the function
A <- discountedloanamount(L,r,t)

result <- paste("A discounted loan of $",L," at an annual interest rate of ",r,"%, after ",t," years, will require repayment in the amount of: $", round(A,2),sep="") 

#result

Result
A discounted loan of $10000 at an annual interest rate of 8%, after 0.25 years, will require repayment in the amount of: $10204.08

Compound Interest

The formula for compound interest can be defined as:

The amount $A$ after $t$ years due to a principal $P$ invested ant an annual interest rate $r$ compounded $n$ times per year is:

$A=P*(1+\frac{r}{n})^{nt}$

#create discounted loan function
compoundinterest <- function (P,r,t,n){
 r <- r/100
 i <- (1+(r/n))^(n*t)
 A <- P * i
}

#amount of loan
P <- 300000

#% rate of loan
r <- 3.2

#number of years
t <- 30

#number of years
n <- 12

A <- compoundinterest(P,r,t,n)

result <- paste("A loan of $",P," at an annual interest rate of ",r,"%, after ",t," years, compounding ",n," times per year will cost a total of: $", round(A,2),sep="") 

#result

Result
A loan of $300000 at an annual interest rate of 3.2%, after 30 years, compounding 12 times per year will cost a total of: $782508.47

Continuous Compounding Interest

The amount $A$ after $t$ years due to a principal $P$ invested ant an annual interest rate $r$ compounded continuously is:

$A=Pe^{nt}$

Where $e$ is Eulers constant (~ 2.718281827), which can be assigned precisely using the R function exp(1).

#create simple interest function
continuouscompoundinterest <- function (P,r,t){
  r <- r/100
  e <- exp(1)
  I <- P * e ^ r * t
  I
}

#amount of principal
P <- 3000

#% rate of loan
r <- 8

#number of years
t <- 1

#interest (unknown)
I <- 0

#run the function
I <- continuouscompoundinterest(P,r,t)

# A = Principal plus Interest
A <- I

#write a sentence to summarize the result
result <- paste("A loan of $",P," at an annual interest rate of ",r,"%, after ",t," years, compounding continuously will cost a total of: $", round(A,2),sep="") 

#result

Result
A loan of $3000 at an annual interest rate of 8%, after 1 years, compounding continuously will cost a total of: $3249.86

Present Value

The present value $P$ of $A$ dollars to be received after $t$ years assuming a per annum interest rate $r$ compounded $n$ times per year is given by:

$P=A\cdot(1+\frac{r}{n})^{-nt}$

#create simple interest function
presentvalue <- function (A,r,n,t){
  r <- r/100
  i <- r / n
  P <- (1+i)^(-n*t)
  P * A
  
}

#amount we want to get to
A <- 900000

#% rate of return on loan
r <- 6

# number of payments
n <- 12

# time in years
t <- 30

#interest (unknown)
P <- 0

#run the function
P <- presentvalue(A,r,n,t)

#write a sentence to summarize the result
result <- paste("The amount that should be invested now at ",r,"% per annum to reach a savings goal of $",A," after ",t," years with interest compounding ",n," times per year is: $",round(P,2),sep="") 

#result

Result
The amount that should be invested now at 6% per annum to reach a savings goal of $900000 after 30 years with interest compounding 12 times per year is: $149437.74

Present Value with Continuous Compounding

The present value $P$ of $A$ dollars to be received after $t$ years assuming a per annum interest rate $r$ with continuous compounding is given by:

$P=A\cdot(1+\frac{r}{n})^{-nt}$

#create simple interest function
presentvaluecontinuous <- function (A,r,t){
  r <- r/100
  e <- exp(1)
  P <- A * e ^ -(r * t)
  P
  
}

#amount we want to get to
A <- 10000

#% rate of return on loan
r <- 4

# time in years
t <- 2

#interest (unknown)
P <- 0

#run the function
P <- presentvaluecontinuous(A,r,t)

#write a sentence to summarize the result
result <- paste("The amount that should be invested now at ",r,"% per annum to reach a savings goal of $",A," after ",t," years with interest compounding continuously is: $",round(P,2),sep="") 

#result

Result
The amount that should be invested now at 4% per annum to reach a savings goal of $10000 after 2 years with interest compounding continuously is: $9231.16

Amount of an Annuity

Assuming $P$ is the deposit made at the end of each payment period for an annuity paying an interest rate of $i$ per payment period. The amount $A$ of an annuity after $n$ deposits is given by:

$A=P\cdot[\frac{(1+i)^{n} - 1}{i}]$

Assuming that $i$ = $\frac{r}{n}$ where $r$ is the annual rate and $n$ is the number of compounding periods per year.

#create simple interest function
amountofannuity <- function (P,r,y,n){
  r <- r/100
  i <- r / y
  z <- ((1+i)^n)-1
  A <- P * (z/i)
  A
}

#amount of principal
P <- 100

#% rate of loan
r <- 5

#number of compounding period per year
y <- 12

# number of payments
n <- 8

#interest (unknown)
A <- 0

#run the function
A <- amountofannuity(P,r,y,n)

# A = Principal plus Interest
A <- A

#write a sentence to summarize the result
result <- paste("After making ",n," deposits of $",P," each, an annuity paying ",r,"%, compounded ",y," times per year, is worth: $",round(A,2),sep="") 

#result

Result
After making 8 deposits of $100 each, an annuity paying 5%, compounded 12 times per year, is worth: $811.76

Present Value of an Annuity

Suppose an annuity earns at an interest rate of $i$ per payment period. If $n$ withdrawals of $P$ are made at each payment period, the amount $V$ required is:

$V=P\cdot[\frac{1-(1+i)^{-n}}{i}]$

Assuming that $i$ = $\frac{r}{n}$ where $r$ is the annual rate and $n$ is the number of compounding periods per year.

#create simple interest function
presentvalueofannuity <- function (P,r,y,n){
  r <- r/100
  i <- r / y
  z <- (1-(1+i)^-n)
  V <- P * (z/i)
  V
}

#amount of principal
P <- 1000

#% rate of loan
r <- 5

#number of compounding period per year
y <- 1

# number of payments
n <- 4

#interest (unknown)
V <- 0

#run the function
V <- presentvalueofannuity(P,r,y,n)

# A = Principal plus Interest
V <- V

#write a sentence to summarize the result
result <- paste("To make ",n," payments of $",P," each, while receiving the proceeds of an annuity paying ",r,"%, compounding ",y," times per year, the amount needed to invest in the annuity now is: $",round(V,2),sep="") 
options(scipen=999)

#result

Result
To make 4 payments of $1000 each, while receiving the proceeds of an annuity paying 5%, compounding 1 times per year, the amount needed to invest in the annuity now is: $3545.95

Amortization

The payment $P$ required to pay off a loan of $V$ dollars borrowed for $n$ payment periods at a rate of interest $i$ per payment period is given by:

$P=V\cdot[\frac{i}{1-(1+i)^{-n}}]$

Assuming that $i$ = $\frac{r}{n}$ where $r$ is the annual rate and $n$ is the number of compounding periods per year.

#create simple interest function
amortization <- function (V,r,y,n){
  r <- r/100
  i <- r / y
  z <- (1-(1+i)^-n)
  P <- V * (i/z)
  P
}

#amount of principal
V <- 500000

#% rate of loan
r <- 4

#number of compounding period per year
y <- 12

# number of payments
n <- 12 * 30

#interest (unknown)
P <- 0

#run the function
P <- amortization(V,r,y,n)

#write a sentence to summarize the result
result <- paste("The payment required to pay off a loan of ",V," dollars borrowed for ",n," payment periods at a rate of interest ",r," is: $",round(P,2),sep="") 

#result

Result
The payment required to pay off a loan of 500000 dollars borrowed for 360 payment periods at a rate of interest 4 is: $2387.08

Double an Investment

The number of years $t$ required to multiply an investment $m$ times when $r$ is the interest rate compounded $n$ times per year is given by:

$t=\frac{\ln(m)}{n \cdot \ln(1+\frac{r}{n})}$

#create simple interest function
multiplyinvestment <- function (m,r,n){
  r <- r/100
  top <- log(m,base=exp(1))
  bot <- n * log(1+(r/n),base=exp(1))
  t <- (top/bot)
  t
}


# % interest rate
r <- 4

# number of times to multiple initial investment
m <- 2

# number of compound periods per year
n <- 1

# time
t <- 0

#run the function
t <- multiplyinvestment(m,r,n)


#write a sentence to summarize the result
result <- paste("To multiply an investment by ",m," times, when ",r,"% is the interest rate, compounded ",n," times per year, it will take: ", round(t,1) ," years."  ,sep="") 

#result

Result
To multiply an investment by 2 times, when 4% is the interest rate, compounded 1 times per year, it will take: 17.7 years.

Time to Reach a Savings Goal with Continuous Compounding

The number of years $t$ required to multiply an investment $m$ times when $r$ is the interest rate compounded $n$ times per year:

$t=\frac{\ln(A)\cdot\ln(P)}{r}$

#create simple interest function
timetoreachgoal <- function (A,P,r){
  r <- r/100
  ln_a <- log(A,base=exp(1))
  ln_p <- log(P,base=exp(1))
  t <- (ln_a - ln_p) / r
 
}


# goal amount
A <- 100000

# initial investment
P <- 20000

# interest rate
r <- 5

# time
t <- 0

#run the function
t <- timetoreachgoal(A,P,r)

t <- round(t,1)

#write a sentence to summarize the result
result <- paste("It will take ",t," years to reach a savings goal of ",A,", earning ",r,"% interest compounding continuously on a principal investment of $",P,"." ,sep="") 

#result

Result
It will take 32.2 years to reach a savings goal of 100000, earning 5% interest compounding continuously on a principal investment of $20000.

Application: Using sapply to visualize compound interest

One of the key advantages of R is that it is a functional language, so a function we have designed can be applied to a vector and return an object of a specified type. The apply family of functions is much faster than using loops to do the same thing.

Money

We will consider the compound interest formula we defined earlier - lets say we want to build a chart that shows what your return on investment would be after 1, 2, 3…30 years. You could certainly loop through from 1 to 30 and build an output vector to show these values.

But instead we will use the sapply function to provide the same output, using just a couple lines of code. sapply takes as it’s input a vector, and returns a vector as output - this is probably the simplest demonstration of the apply family of functions.

#wrap the compoundinterest in a function that defines all the values other than time
investment <- function(t){
  P <- 10000
  r <- 4
  n <- 12
  compoundinterest(P,r,t,n)
}
#create a vector of years
years <- seq(1:30)

#run sapply and get back a vector with the value of your investment after each year
list_earned <- sapply(years,investment)
list_earned

##  [1] 10407.42 10831.43 11272.72 11731.99 12209.97 12707.42 13225.14 13763.95
##  [9] 14324.72 14908.33 15515.72 16147.85 16805.74 17490.43 18203.02 18944.64
## [17] 19716.47 20519.75 21355.75 22225.82 23131.34 24073.74 25054.54 26075.30
## [25] 27137.65 28243.28 29393.96 30591.51 31837.86 33134.98

#create a data frame with the years and amounts
invest_values <- data.frame(year=years,earned=list_earned)

earned_plot <- ggplot(data=invest_values, aes(x=year, y=earned)) + 
  geom_point(size=1) + 
  geom_line() + 
  labs(
    title="Compound Interest Chart",
    subtitle="Growth of a $10,000 investment, over 30 years at 4%, compunding monthly"
  )
earned_plot

Application: Using mapply to plot a ski run by segment

A more powerful version of the apply family is mapply. This function can take multiple vectors as input and return a list of outputs. For this exercise, we will use our ski slope application above; but this time, we will conisider a more complicated scenario. Instead of plotting just one, straight slope, we want to understand the difficulty of a slope with multiple way points.

To accomplish this, we will pass in to this function a vector of heights and a vector of distances. What we receive back form the mapply function will give us a rating for of each section of the ski slope.

#create a vector of years
height_segments <- c(200,700,400,1200,600,100)
distance_segments <- c(1500,1400,300,1250,750,1000)

#run mapply with out get_rating function defined earlier
list_segments <- mapply(get_rating,height_segments,distance_segments)
list_segments <- data.frame(list_segments)

Now we will use the returned data to plot a line chart which shows the difficulty of each portion of the slope.

#first we will transpose the list with names as rows, to instead assign the names to columns
list_rows <- t(list_segments)
#convert the list to a data frame
df<-as.data.frame(list_rows,stringsAsFactors=FALSE)

#set the data type for each row and fortify with some extra data for our chart, including the distance fromm the top in elevation and length of each section
df$height <- as.numeric(df$height)
df$length <- as.numeric(df$length)
df$rating <- as.character(df$rating)
df$colorcode <- as.character(df$color)
df$angle <- as.numeric(df$angle)
df$distance <- as.numeric(df$distance)

#the cumsum will cumulatively sum the data frame rows for our chart
df$heightfromtop <- cumsum(-df$height)
df$distancefromtop <- cumsum(df$distance)

#the lag will take teh last value from each row to [provide a more convenient way to set up geom_segment objects in our GGplot
df$heightfromtopstart <- lag(df$heightfromtop,1)
df$distancefromtopstart <- lag(df$distancefromtop,1)

#replace the NA values from the first row of teh cumsums
df[is.na(df)] = 0

#peek at teh table
df

We can see that the table now contains all the data we need to create a nice looking plot showing each segment of the ski run.

ski_plot <- ggplot(data=df) +   
  geom_segment(size=1,aes(y=heightfromtopstart, 
                          yend=heightfromtop, 
                          x=distancefromtopstart, 
                          xend=distancefromtop,
                          color=rating)) +
  labs(
    title="Ski Run profile",
    subtitle="Difficulty by segment",
    y="Elevation (ft)",
    x="Distance from Top (ft)"
  )
ski_plot

Economics

The following formulas are drawn from the discipline of finite mathematics, for application to the study of economics. These are just a few of the useful, basic formulas from the vast field of economic theory.

Compute Inflation

When the rate of inflation averages $r$ % over $n$ years, the value of your initial amount of Principal will be less by $(1-r)^n$ over time:

$A=P\cdot(1-r)^n$

#create simple interest function
computeinflation <- function (P,r,n){
  r <- r/100
  A <- P * (1-r)^n
  A
 
}

#principal
P <- 1000

# % avg interest rate
r <- 2

# number of year
n <- 3

#run the function
A <- computeinflation(P,r,n)

#write a sentence to summarize the result
result <- paste("Assuming an average annual inflation rate of ",r,"%, $",P," after " ,n," years will be worth ",round(A,2),"."  ,sep="") 

#result

Result
Assuming an average annual inflation rate of 2%, $1000 after 3 years will be worth 941.19.

Compute the Consumer Price Index (CPI)

If the rate of inflation averages $r$ over $n$ years, the the CPI index after $n$ years is:

$CPI=CPI_0(1+\frac{r}{100})$

#create simple interest function
computecpi <- function (cpi0,r,n){
  cpi <- cpi0 *((1+(r/100))^n)
  cpi
 
}
#CPI start year
cpiyear0 <- 2010

#CPI in 2010
cpi0 <- 217.6

# % avg interest rate
r <- 3

# number of year
n <- 10

# second year
cpiyear1 <- cpiyear0 + n

# time
cpi1 <- 0

#run the function
cpi1 <- computecpi(cpi0,r,n)


#write a sentence to summarize the result
result <- paste("Assuming an average annual inflation rate of ",r,"%, the CPI in ",cpiyear1," (after " ,n," years) would be ",round(cpi1,1)," given a starting CPI of ",cpi0," in ",cpiyear0,"."  ,sep="") 

#result

Result
Assuming an average annual inflation rate of 3%, the CPI in 2020 (after 10 years) would be 292.4 given a starting CPI of 217.6 in 2010.

Leontief Models

The Nobel Prize winning Mathematician / Economist Wassily Leontief developed an important field of economic theory pertaining to how resources must be used in the production of other resources, or require some of their output to be used in the production or delivery of that resource. As an example of how a Leontief model can be used, consider the following scenario…

Application: Building a City on an Island

The director of a development consortium plans to build a city on an island in the middle of the ocean. This island, closed off from the outside world, will need 3 key resources to be constructed: Food, Buildings, and Fuel - which must all be produced by workers on that island. Moreover, all of the productivity of these workers will be 100% used up by this undertaking - this is a “closed” Leontief model- all inputs minus all outputs must equal 0 - no food, fuel or buildings will be exported or imported.

Island

Before agreeing to undertake this venture, a group of of farmers, builders and refiners agrees to pool their efforts to accomplish this goal. Their investors agree to provide a total of 30 billion dollars to this task. They further agree that each party- farmers, builders, and refiners should receive roughly equal pay for the performance of their important duties. So based on the investment available, each party should receive about 10 billion dollars for their efforts in this herculean undertaking… but the director makes an important observation - not all parties will necessarily produce, or use the same amount of resources, and some should receive more or less compensation, according to their true relative contribution to this undertaking.

The director points out that it requires fuel to operate a refinery to make more fuel. A tanker truck also require some fuel to transport crude oil to the refinery, and to bring the finished fuel to its final destination at the gas pump. But the driver and refinery workers also need to eat, and be housed. Thus the farmers must produce food, some of which they consume, and some of which, they provide to the refiners and the builders. The farmers in turn, must use tractors, which will use some of the fuel. But then everyone must also have a place to sleep, thus all of these producers require some quantity of the production output of the builders, and the builders, must also use fuel and food in the production of houses. In this closed system each party is using some of their own productivity, and some of the productivity of the other parties, at an uneven rate.

Determining the amount of fuel, food, and houses needed as an input, to bring a certain amount of fuel, food, and houses to market as an output, requires a Leontief model.

For this demonstration we will consider a simplified version of this economy or closed model that considers just these three resources - food, houses, and fuel, and assumes that all resources are consumed or where input is equal to output. We will consider, in the context of a closed economy, what is the relative value of these three resources based on how much each party produces and uses.

produced_by <- c('used_by_farmer','used_by_builder','used_by_refiner')
food_used_by <- c(.4,.4,.2)
houses_used_by <- c(.4,.3,.3)
fuel_used_by <- c(.3,.2,.5)

io_model <- data.frame(input_output=produced_by,
                       produced_by_farmer=food_used_by,
                       produced_by_builder=houses_used_by,
                       produced_by_refiner=fuel_used_by)

io_kable <- knitr::kable(io_model) %>%
  kable_styling(font_size = 12) %>% 
  row_spec(0, bold = TRUE, italic = TRUE, color = 'black', background = '#CCCCFF')

To represent the table above, we will construct a matrix with the following form, where $A$ is the input output matrix:

$A=\begin{bmatrix}.4 & .4 & .3\\.4 & .3 & .2\\.2 & .3 & .5\end{bmatrix}$

We will now define the variables for each workers pay as:

$x = Farmers Pay$

$y = Builders Pay$

$z = Refiners Pay$

And solve for these variables using the following system of equations:

$\begin{bmatrix}x\\y\\z\end{bmatrix}=\begin{bmatrix}.4 & .4 & .3\\.4 & .3 & .2\\.2 & .3 & .5\end{bmatrix} \begin{bmatrix}x\\y\\z\end{bmatrix}$

The system above is referred to as $X = AX$

Simplification of this system results in: $(I-A)X = 0$

Where $I$ is the 3x3 identity matrix:

$I=\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{bmatrix}$

#diag function returns an identity matrix of dimension specified
identity_matrix <- diag(3)

#concatenate the 3 values in the table created in the previous section
concat_prod <- c(food_used_by,houses_used_by,fuel_used_by)

#convert the vector into a matrix
io_matrix <- matrix(concat_prod,3,3) 

#subtract the input output matrix from the I matrix to get the difference
difference_matrix = identity_matrix - io_matrix

Given that:

$I-\begin{bmatrix}.4 & .4 & .3\\.4 & .3 & .2\\.2 & .3 & .5\end{bmatrix}=\begin{bmatrix}.6 & -.4 & -.3\\-.4 & .7 & -.2\\-.2 & -.3 & .5\end{bmatrix}$

#use the identity matrix as the RHS to determine the relationship of X and Y to Z
R <- pracma::rref(difference_matrix)   
R

##      [,1] [,2]       [,3]
## [1,]    1    0 -1.1153846
## [2,]    0    1 -0.9230769
## [3,]    0    0  0.0000000

#determine the proportion of Y and X to Z
X_proportion_to_Z <- 1/R[1,3]
X_to_Z <- -1/X_proportion_to_Z

Y_proportion_to_Z <- 1/R[2,3]
Y_to_Z <- -1/Y_proportion_to_Z

#set the target amount designated to Z
Z_pay <- 1000

#Determine the amount to be designated to X
X_pay <- Z_pay * X_to_Z
X_pay

## [1] 1115.385

#Determine the amount to be designated to Y
Y_pay <- Z_pay * Y_to_Z
Y_pay

## [1] 923.0769

result <- paste("Assuming that resources are used in the proportions indicated in the table above, farmers should receive $",round(X_pay,2)," and builders should receive $",round(Y_pay,2)," for each $",round(Z_pay,2)," of compensation to refiners."  ,sep="") 

#result
#

io_kable

input_output	produced_by_farmer	produced_by_builder	produced_by_refiner
used_by_farmer	0.4	0.4	0.3
used_by_builder	0.4	0.3	0.2
used_by_refiner	0.2	0.3	0.5

Result
Assuming that resources are used in the proportions indicated in the table above, farmers should receive $1115.38 and builders should receive $923.08 for each $1000 of compensation to refiners.

Sets

In this section we will set up a function to perform set operations and analysis.

setoperations <- function (A,B){
  
#Set teh length variables
A_len <- length(A)
B_len <- length(B)

#get teh union
A_un_B <- unique(c(A,B))

#get the intersection
A_in_B <- A[which((A%in%B))]
B_in_A <- B[which((B%in%A))]

#determine set rules
A_sub_B <- FALSE
B_sub_A <- FALSE
B_super_A <- FALSE
A_super_B <- FALSE
A_psub_B <- FALSE
B_psub_A <- FALSE
B_psuper_A <- FALSE
A_psuper_B <- FALSE
A_disjoint_B <- FALSE

if(length(A_in_B)>0){
  #if n(A) is < n(B)
  if(A_len < B_len){
    #If the elements of A that are in B are all the elements of A
    if(length(A_in_B) == length(A)){
    #A is a proper subset of B
    A_psub_B <- TRUE
    B_psuper_A <- TRUE
    }
  }
  #if n(A) is < n(B)
  if(A_len <= B_len){
    #If the elements of A that are in B are all the elements of A
    if(length(A_in_B) == length(A)){
    #A is a proper subset of B
    A_sub_B <- TRUE
    B_super_A <- TRUE
    }
  }
  
  #if n(A) is < n(B)
  if(B_len < A_len){
    #If the elements of A that are in B are all the elements of A
    if(length(B_in_A) == length(B)){
    #A is a proper subset of B
    B_psub_A <- TRUE
    A_psuper_B <- TRUE
    }
  }
  #if n(A) is < n(B)
  if(B_len <= A_len){
    #If the elements of A that are in B are all the elements of A
    if(length(B_in_A) == length(B)){
    #A is a proper subset of B
    B_sub_A <- TRUE
    A_super_B <- TRUE
    }
  }
  
}
if(length(A_in_B) == 0){
  A_disjoint_B <- TRUE
}

A_intersect_B <- A_in_B
A_union_B <- A_un_B

return_obj <- list(
  A = A,
  B = B,
  A_in_B = A_sub_B,
  B_in_A = B_sub_A,
  A_sub_B = A_sub_B,
  B_sub_A = B_sub_A,
  A_super_B = A_super_B,
  B_super_A = B_super_A,
  A_psub_B = A_psub_B,
  B_psub_A = B_psub_A,
  A_psuper_B = A_psuper_B,
  B_psuper_A = B_psuper_A,
  A_union_B = A_union_B,
  A_intersect_B = A_intersect_B,
  A_disjoint_B = A_disjoint_B

)

return_obj

}

Now we will create a couple of small sets to show how the setoperations function can aid in the analysis of 2 simple sets. IN later demonstrations, will try the same function against some larger sets.

#create two small vectors (sets) A and B, where A is a proper subset of B
A <- c(3,5)
B <- c(1,3,5,7,9)
set_ops <- setoperations(A,B)

From the two vectors, A and B, defined above we can determine the following set properties in LaTeX. The setoperations function we have defined contains all the common set length analysis, equality tests including subsets, supersets, intersections, and unions and tests for determination of disjoint and intersection.

Proper Subset / Superset Tests:

A proper subset is a set where the set being tested must be entirely contained within the set being compared to. However, the set being tested may not be equal to the set being compared to. At least one element must be in the set being compared to, which is not in the set being tested.

In the exaMple given, set A is a proper Subset of Set B, since 3 and 5 are both in 1,3,5,7,9, but 3 and 5 are not equal to 1,3,5,7,9.

A Proper Subset operator can be expressed in LaTeX using the “subset” operator where the set on the left is a proper subset of the set on the right.

For these tests we will examine the returned value of A_psub_B and B_psub_A, as returned by the setoperations function we defined above.

$A \subset B$ TRUE

$B \subset A$ FALSE

A Proper Superset can be expressed using the “supset” operator, where the set on the left is a proper superset of the set on the right.

For these tests we will examine the returned value of A_psuper_B and B_psuper_A, as returned by the setoperations function we defined above.

$A \supset B$ FALSE

$B \supset A$ TRUE

Subset / Superset Tests:

A Subset operator (with possible equality) can be expressed in LaTeX using the “subseteq” operator where the set on the left is a subset of the set on the right.

For these tests we will examine the returned value of A_sub_B and B_sub_A, as returned by the setoperations function we defined above.

$A \subseteq B$ TRUE

$B \subseteq A$ FALSE

A Superset operator (with possible equality) can be expressed using the “supseteq” operator, where the set on the left is a superset of the set on the right.

For these tests we will examine the returned value of A_super_B and B_super_A, as returned by the setoperations function we defined above.

$A \supseteq B$ FALSE

$B \supseteq A$ TRUE

Since A is a proper subset of B it is also a regular subset of B.

Intersection Test

An intersection of two sets can be expressed using the “cap” operator in LaTeX. In this case, the intersection of sets A and B is the same as set A since A is entirely contained in B (A is a proper subset of B).

For this test we will examine the returned value of A_intersect_B as returned by the setoperations function we defined above.

$A \cap B$ 3, 5

Union Test

A union of two sets can be expressed using the “cup” operator in LaTeX. In this case, the union of sets A and B is the same as set B since A is entirely contained in B (A is a proper subset of B).

For this test we will examine the returned value of A_union_B as returned by the setoperations function we defined above.

$A \cup B$ 3, 5, 1, 7, 9

Disjoint Test

The disjoint sets can be indicated in latex by adding a “dot” above the “cup” operator.

To test of the sets are disjoint you can just check the A_disjoint_B property of the returned list.

$\dot\cup$ FALSE

So clearly A and B are not disjoint, since A is a subset of B. If instead we set the two sets to even and odd numbers will find the sets are disjoint - they share no common elements.

#create two small vectors (sets) A and B, where A is a proper subset of B
A <- c(2,4,6,8,10)
B <- c(1,3,5,7,9)
set_ops <- setoperations(A,B)

Testing again if the sets are disjoint using A_disjoint_B property of the returned list we find:

$\dot\cup$ TRUE

Because there are no common elements between the odd and even numbers under ten.

Creating Sets using seq:

In practice its not always very efficient to create sets by listing out every member of your sets. R makes it easy to create sets using methods that can produce specific sets, sets based on rules, and random sets.

OK, so lets imagine that in the example above, instead of listing the odd numbers between 1 and 9, we want to generate all the odd numbers between 0 and 100. THis would be a lot of numbers to type out! For this purpose, R provides the seq function. This function will generate a list of numbers between a lower and upper bound, stepping by a specified amount. IN this test we will instead ask, if I counted from 0 to 100 by 3’s and 0 to 100 by 2s, how many numbers would be shared between both sets?

A <- seq(from = 0, to = 100, by =3)
B <- seq(from = 0, to = 100, by =2)
set_ops <- setoperations(A,B)

In the example above the numbers in A that are also in B are:

$A \cap B$ 0, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84, 90, 96

Note that this is equivalent to answering the question: What are the common multiples of 2 and 3? SO we can see that using pair of sets provides a way to test the arithmetical rules that define basic algebra.

Another way to use seq is to simply specify the starting number and how many items you want to return with the length.out argument. This version will return only 20 elements for each set.

A <- seq(from = 0, by =5, length.out =20)
B <- seq(from = 0, by =10, length.out =20)
set_ops <- setoperations(A,B)

The union of both sets contains:

$A \cup B$ 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190

But the intersection only contains:

$A \cap B$ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90

Because the sets were limited in length only 20 elements - so only the10’s from 0 to 90 are found in both sets.

Creating Sets using rnorm:

In the example below, we will us the rnorm to generate 2 sets with a random normal distribution. This code generates 2 sets of 50 random numbers, centered around the number 100, with a standard deviation of 25.

#Generate two random sets centered around the number 50, with a population mean of 100, and standard deviation of 50
set.seed(1)
A <- unique(round(rnorm(n=50,mean=100,sd=25),0))
B <- unique(round(rnorm(n=50,mean=100,sd=25),0))
set_ops <- setoperations(A,B)

In the example above the numbers in A that are also in B are:

$A \cap B$ 105, 140, 108, 112, 118, 114, 92, 110, 99, 100, 115, 102, 88, 97, 117, 82, 109

So we can see that many of the same numbers were generated, given the relatively low standard deviation of the two sets.

Not that there are also ways to generate sets using other distribution models that are not the normal distribution, such as pnorm, dnorm, and qnorm, which will be explored in later installments.

Tip: Using the set.seed function before a block of code that calls a randomization function like rnorm, is CRITICAL for reproducible code! Set the randomization seed to any number to ensure that your random generator will always produce the same set of numbers when that code is run multiple times or on different machines. Changing the number parameter in the set seed function will change the outcome of the randomization - but be sure you don’t “cherry pick” the random results you want by trying different random seeds. The number passed in as a seed should, in itself, be random!

Application: Random Primes Using the Sieve of Eratosthenes

In this simple application we will use the RcppAlgos package which uses and algorithm based on the Sieve of Eratosthenes to generate a list of all Prime Numbers between 0 and 1000 - and then generate a list of random numbers centered around the number 500, and get the intersection of the two lists.

library(RcppAlgos)
set.seed(99)
A <- primeSieve(1000)
B <- unique(round(rnorm(n=500,mean=500,sd=250),0))
set_ops <- setoperations(A,B)

In this example above the prime numbers in A that are also in a random number of selections sampled from set B are:

$A \cap B$ 19, 43, 89, 149, 157, 197, 211, 229, 239, 269, 271, 281, 283, 293, 307, 313, 317, 337, 347, 383, 389, 397, 401, 409, 421, 433, 439, 443, 457, 461, 463, 467, 491, 499, 509, 523, 569, 587, 607, 613, 641, 677, 709, 787, 853, 967

So we can see that set functions in R can provide many interesting ways to combine vectors of numbers that follow various rules to generate highly unique sets!

Number Systems:

The demonstrations above are great for cases when your contain a “bounded” value - where $U < \infty$ But what if you want to represent a set of unbounded numbers?

In this case, you will need some special notation, and symbols that are used in set theory. These are shown in latex using the blackboard bold modification, which resembles a bold letter that has been “double struck”. You can use the mathbb{} syntax in latex to add blackboard bold to any character.

Here’s a few of the most common examples of these sets, and what they mean:

Real Numbers

$\mathbb{R}$

The set of all real numbers, from negative infinity to positive infinity, including integers, fractions, irrational numbers, and 0.

Now, that’s a pretty big set… it would not fit well into a Venn Diagram!

In addition to $\mathbb{R}$ the mother of all sets - there are many well established conventions for how to represent subsets of $\mathbb{R}$

Natural Numbers

$\mathbb{N}$

The set of all natural numbers from 1 to positive infinity:

$\{1,2,3,...,\infty\}$

Whole Numbers

$\mathbb{W}$

The set of all whole numbers, from 0 to positive infinity.

$\{0,1,2,3,...,\infty\}$

Integers

$\mathbb{X}$

The set of all whole numbers, from negative infinity to positive infinity, including 0.

$\{-\infty,...,-3,-2,-1,0,1,2,3,..., \infty\}$

Rational Numbers

$\mathbb{Q}$

The set of numbers in $\mathbb{R}$, with the exception of irrational numbers (Like the square root of 2)

The above list is by no means comprehensive, but these numbers systems should be useful in describing unbounded sets that contain the more common number sets. Really, anything can be a set- but fear not, we don’t need an infinite number of blackboard bold letters (whew!) - as there is also a convention for mathematical notation syntax related to sets that may not follow such an easily described group.

Descriptive Notation of Sets

Let’s consider cards in a standard deck. In mathematical notation, we don’t really have good way to describe these, as they have symbols like hearts, clubs, diamonds and spades.

We could perhaps write out a very long list of 52 values like:

{2 of Hearts,3 of Hearts,4 of Hearts,…,Ace of Spades}

BUt this is really not very user-friendly, and you are going to waste a lot of graph paper doing so. (Stuff’s expensive, don’t want that, right?!)

So when describing a set that contains all the cards in a standard deck, we just write something that would be easy for someone to understand using notation like this:

{X | x = Any cards in a standard deck of playing cards}

So, in the example above we have noted that “Big X” is the set of “Little x”, where Little x is any one of the cards in a deck of cards.

So little x could be Ace of Spades, 2 of Clubs, King of Diamonds, etc. and Big X is the entire deck.

More examples

{J | j = Any even number}

{P | p = Any prime number less than 100}

Set Related Rules and Formulas

Having explored some of the ways we can construct and represent sets in R and Latex, we will now consider how formulas and rules can be used to confirm the validity and integrity of our sets.

Counting Formula

The counting formula holds that the number of elements in A union B is equal to the number of elements in A, plus the number of elements in B, minus the number of elements in A intersect B:

$n(A \cup B) = n(A) + n(B) - n(A\cap B)$

We will prove that our set operations function is correct by using the counting formula. If the set operation is correct then the same_length should be TRUE between A_union_B and the product of the counting formula.

#check the length of A + B - A intersect B
len_A_union_B_1 <- length(set_ops$A) + length(set_ops$B) - length(set_ops$A_intersect_B) 

#Check the length o A union B
len_A_union_B_2 <- length(set_ops$A_union_B)

#compare the length
same_length <- len_A_union_B_1 == len_A_union_B_2
same_length

## [1] TRUE

The counting formula agrees with the length of A union B.

De Morgan’s Properties

British Mathematician Augustus De Morgan defined a set of properties which provide a mechanism in boolean algebra to describe the relationship of sets and their complements, providing rules to define conjunctions and disjunctions, and ensuring consistent and rigorous rules of inference related to sets.

These are often expressed using the ever popular Venn Diagram, which is explored in the next section.

De Morgan observed that:

The complement of the union of two sets is the same as the intersection of their complements:

$\overline{A \cup B} = \overline{A} \cap \overline{B}$

The complement of the intersection of two sets is the same as the union of their complements:

$\overline{A \cap B} = \overline{A} \cup \overline{B}$

Venn Diagrams

Now we will see how De Morgans principles can be applied by using the Venn Diagram, which is drawn as 2 or more overlapping circles. To facilitate the creation of these diagrams, we will develop a function using our trusty companion in all things graphical, GGPlot!

To accomplish this, we will build upon the setoperations function by drawing 2 circles that are placed in different locations on the GGplot canvas, depending on whether the sets intersect, are disjoint, or are subsets of each other.

venndiagram <- function (A,B){

set <- setoperations(A,B)
  
x_points <- rep(1:10)
y_points <- rep(1:10)
coords <- data.frame(x=x_points,y=y_points)

x_points <- c(3,6)
y_points <- c(4,5)
circles <- data.frame(x=x_points,y=y_points)

if(!set$A_disjoint_B){
  # the groups are not disjoint, so they must overlap
  if(set$A_sub_B | set$B_sub_A){
      #if one set is contained in the other
   
    if(set$A_sub_B){
      circlecolor <- c("A","B")
      circlex <- c(2,2)
      circley <- c(3,3)
      circlesize <- c(1.5,2)
    }else if(set$B_sub_A){
      circlecolor <- c("A","B")
      circlex <- c(2,2)
      circley <- c(3,3)
      circlesize <- c(2,1.5)
    }else{
      circlecolor <- c("A","B")
      circlex <- c(1,3)
      circley <- c(3,3)
      circlesize <- c(2,2)
      
    }
    
      circles <- data.frame(
          x0 = circlex,
          y0 = circley,
          r = circlesize,
          fill = circlecolor
        )
      
      totals <- data.frame(
        x = c(0,4),
        y = c(5,5),
        sides = c('A','B'),
        num = c(length(A),length(B))
      )
      
      counts <- data.frame(
        x = c(0.25,3.75),
        y = c(3,3),
        sides = c('A','B'),
        num = c(length(A)-length(set$A_intersect_B),length(B)-length(set$A_intersect_B))
      )
      
      intersect <- data.frame(
        x = c(2),
        y = c(3),
        sides = c('A∩B'),
        num = c(length(set$A_intersect_B))
      )
      
      U <- length(set$A) + length(set$B) - length(set$A_intersect_B) 
      
      # plot
      ggplot() +
        geom_circle(aes(x0 = x0, y0 = y0, r = r, fill=fill), alpha=.5, data = circles) +
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = totals)+
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = counts)+
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = intersect)  +
        annotate("text", x=-.5,y=5,label = paste("U:",U)) +
       labs(title="Sets with Subsets")
      
  }else{
     #neither set is contained in the other
       
    
    circles <- data.frame(
          c = c("A","B"),
          x0 = c(1,3),
          y0 = c(3,3),
          r = c(2,2)
        )
      
      totals <- data.frame(
        x = c(0,4),
        y = c(5,5),
        sides = c('A','B'),
        num = c(length(A),length(B))
      )
      
      counts <- data.frame(
        x = c(0,4),
        y = c(3,3),
        sides = c('A','B'),
        num = c(length(A)-length(set$A_intersect_B),length(B)-length(set$A_intersect_B))
      )
      
      cap <- TeX('$\\cap$')
      
      intersect <- data.frame(
        x = c(2),
        y = c(3),
        sides = c(paste('A∩B',sep='')),
        num = c(length(set$A_intersect_B))
      )
      
       U <- length(set$A) + length(set$B) - length(set$A_intersect_B) 
      
       
      # plot
      ggplot() +
        geom_circle(aes(x0 = x0, y0 = y0, r = r, fill=c), alpha=.5, data = circles) +
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = totals) +
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = counts) +
        geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = intersect) +
        annotate("text", x=-.5,y=5,label = paste("U:",U)) +
        labs(title="Intersecting Sets")
    
  }


}else{
 # the groups are disjoint, so they do not overlap
 
  circles <- data.frame(
    c = c("A","B"),
    x0 = c(1,5.5),
    y0 = c(3,3),
    r = c(2,2)
  )
  totals <- data.frame(
    x = c(0,4),
    y = c(5,5),
    sides = c('A','B'),
    num = c(length(A),length(B))
  )
  counts <- data.frame(
    x = c(0,4),
    y = c(3,3),
    sides = c('A','B'),
    num = c(length(A)-length(set$A_in_B),length(B)-length(set$A_in_B))
  )
  
 U <- length(set$A) + length(set$B) - length(set$A_intersect_B) 
      
ggplot() +
  geom_circle(aes(x0 = x0, y0 = y0, r = r, fill=c), alpha=.5, data = circles) +
  geom_text(aes(x = x, y = y, label=paste(sides,":",num,sep="")), data = totals) +
  annotate("text", x=-.5,y=5,label = paste("U:",U)) +
  labs(title="Disjoint Sets")

  
}


}

Our venndiagram function above has all the logic we need to draw any possible configuration of 2 sets. IN practice a venn diagram could envisions more than 2 sets, but for this paper we will just focus on 2 sets. Now we will now create some random sets and see if they are intersecting, disjoint, or if one set is a subset of another set.

Further Study: Consider how the function above could be extended to describe more than 2 sets - how could you show 3, 4 or 5 sets by changing the radius of each circle and it’s center point? Where would you need to place the numbers to show the correct intersection zones? Drawing Venn Diagrams is a valuable tool in understanding sets - and this is one area where using a pencil and paper may be a lot easier than writing R code!

Disjoint Sets

In a pair of disjoint sets, none of the elements are shared between the two sets.

#randomly generate a disjoint set
set.seed(12)
A <- unique(round(rnorm(n=5,mean=50,sd=10),0))
B <- unique(round(rnorm(n=5,mean=150,sd=10),0))
venndiagram(A,B)

Looking at the above plot - we can confirm these two sets are disjoint - there are no shared numbers between them. The count of numbers in A is 5, and in B is 4. SInce the intersection of these is 0 - the count of U “the universe” is 9.

## [1] 35 66 40 41 30

## [1] 147 144 149 154

Intersecting Sets

When one or more elements in one set are also in the other set, the sets are intersecting. A venn diagram can help visualize the number of elements that are in set A, set B, or both. Note that the number is the number of elements, not the actual number that occurs in both sets.

#randomly generate a an intersecting set
A <-  c(92,99,97,104,88,112,100,98,96) 
B <-  c( 113,95,101,97,118,92,99,74,109,93)
venndiagram(A,B)

## [1]  92  99  97 104  88 112 100  98  96

##  [1] 113  95 101  97 118  92  99  74 109  93

Looking at the above plot - we can see these two sets overlap - there are 3 shared numbers between them, this is the overlapping region of the Venn diagram. The total count of numbers in A is 9, and in B is 10. The count of U “the universe” is 16 - this follows our counting formulas specified above, because:.

$n(A \cup B) = n(A) + n(B) - n(A\cap B)$

Which agrees with the counts we came up with:

$16 = (9 + 10) - 3$

Set with a Subset

When all of the elements within one set are also in the other set, the sets are intersecting. A Venn diagram can help visualize the number of elements that are contained in the outer superset and the inner subset.

#create a set with a subset
A <- c(1,2,3,4,5,6)
B <- c(3,4,5)
venndiagram(A,B)

## [1] 1 2 3 4 5 6

## [1] 3 4 5

Combinations & Permutations

In this section we will explore the concepts of expected value, using combinations, permutations, and principles of counting.

Multiplication Principle of Counting

When completing a sequence of choices, there are $p$ choices for the first, $q$ for the second, $r$ for the third, and so on, then the multiplication principle of counting may be used to determine how many different ways $D$ that this sequence can be completed.

This formula holds that: $D = p \cdot q \cdot r ...$

This formula is useful for determining things like the number of paths through a maze, the number of item combinations on a menu, or the number of possible members of a committee or focus group, when selecting from multiple pools of candidates.

#set a series of choices that can be combined
choices <- c(4,6,3)

#the total will multiple each choices following the first, byu the sum of all prior products
total <- choices[1]
for (x in choices[-1]){
  total <- total * x
}

# result
result <- paste("A menu with ", choices[1]," appetizer options,  ", choices[2]," entrees, and ", choices[3]," drinks provides a total of ",total," distinct meal choices.",sep="" )
#result

Result
A menu with 4 appetizer options, 6 entrees, and 3 drinks provides a total of 72 distinct meal choices.

Expected Value

To find the expected value, you must define a sample space, determine the probability of each outcome, and assign a payoff (or cost) to each outcome.

$E = (p_1 \cdot m_1) + (p_2 \cdot m_2) + ... + (p_n \cdot m_n)$

Considering the roll of a dice - if you have to pay to roll, and have 1 chance to win, the probability of each outcome is 1 divided by the number of sides on the dice. The list of payouts is 1 times the winning amount, and sides - 1 times cost to play.

If the game results in an expected value of 0 - it is considered a “fair” game. If you played the game enough times, your net earnings would be 0. If the expected value is greater than 0, the game favors the player, but if the value is less than 0, the game favors the house. A game that favors the player will eventually lead to more money won than lost. This is why most casino games are not “fair games” - otherwise the casino would lose money over time.

#expected value function
expected_value <- function (p,m){
 E <- sum(p * m)
 E
}

sides <- 8
side_chance <- 1/sides
payoff <- 3
cost <- -.50

#likelihood of each outcome in teh sample space
p <- rep(side_chance,sides)

#payoff for each outcome
m <- c(payoff,rep(cost,sides-1))
 
E <- expected_value(p,m)

# result
result <- paste("You are playing a dice game. To win $",payoff," you must roll a 1 on a ",sides," sided dice. If it costs $",cost," to play, the expected value is ",round(E,2),sep="" )

Result
You are playing a dice game. To win $3 you must roll a 1 on a 8 sided dice. If it costs $-0.5 to play, the expected value is -0.06

Application: Pot Odds in Texas Hold ’em

In considering casino games - a fair game is one where the expected value is greater than 0. A simple modification to the above function can tell you if you if a game is fair or a bet is good, based on the expected value. In the following code block, we will assess if a final bet in a game of Texas hold ’em is a “good bet” based on pot odds.

Poker

Pot odds in poker aims to assess whether your chance of drawing a winning card is good enough to justify a bet of a certain amount, to win a pot of a certain size. This is just one of the ways that Bayesian probability rules and the Expected Value can be applied to games of chance where the chance to win varies based on the current state of the game.

#expected value function
is_fair <- function (E){
 isfair <- ifelse(E > 0,'YES','NO')
}

#pot_odds function
pot_odds <- function(a,o,w,b){
  
  already_dealt <- a
  cards <- 52 - already_dealt
  outs <- o
  out_chance <- outs/cards
  payoff <- w
  cost <- -b
 
  #create a vector of chances
  p <- rep(out_chance,cards)

  #create a vector of payouts
  m <- c(rep(payoff,outs),rep(cost,cards-outs))
  
  #get teh expected value
  E <- expected_value(p,m)
  
  #call the is_fair evaluation wrapper on the expected value function
  goodbet <- is_fair(E)
}

#Number of dealt cards
already_dealt <- 10
#cards gthat would improve your hand (& Win)
outs <- 9
#size of the pot
payoff <- 200
#cost to call
cost <- 50

#call the function
goodbet <- pot_odds(already_dealt,outs,payoff,cost)

# result
result <- paste("You are playing a game of Texas hold 'em poker, and you are down to one opponent, and one remaining draw (the river). A total of ",already_dealt," cards have already been dealt. It costs $",cost," to call, for a chance to win the $",payoff," pot. You believe that there are ",outs ," cards (outs) left in the deck that will cause you to win. Is this a good bet? According to strict 'pot odds', the answer is: ",goodbet,sep="" )

#result

Result
You are playing a game of Texas hold ‘em poker, and you are down to one opponent, and one remaining draw (the river). A total of 10 cards have already been dealt. It costs $50 to call, for a chance to win the $200 pot. You believe that there are 9 cards (outs) left in the deck that will cause you to win. Is this a good bet? According to strict ’pot odds’, the answer is: YES

#cost to call
cost <- 70

#call the function
goodbet <- pot_odds(already_dealt,outs,payoff,cost)

# result
result <- paste("You are playing a game of Texas hold 'em poker, and you are down to one opponent, and one remaining draw (the river). A total of ",already_dealt," cards have already been dealt. It costs $",cost," to call, for a chance to win the $",payoff," pot. You believe that there are ",outs ," cards (outs) left in the deck that will cause you to win. Is this a good bet? According to strict 'pot odds', the answer is: ",goodbet,sep="" )


#result

Result
You are playing a game of Texas hold ‘em poker, and you are down to one opponent, and one remaining draw (the river). A total of 10 cards have already been dealt. It costs $70 to call, for a chance to win the $200 pot. You believe that there are 9 cards (outs) left in the deck that will cause you to win. Is this a good bet? According to strict ’pot odds’, the answer is: NO

Factorials

To determine the number of possible ordered arrangements of groups of elements, We must use the factorial operator. A factorial is defined as the product of an integer and all the integers below it.

A factorial can be expressed mathematically as:

$n! = n \cdot (n - 1) \cdot (n - 2) \cdot (n - 3) \cdot ... \cdot 3 \cdot 2 \cdot 1$

And the value of $0!$ is 1 per the convention for an Empty Product.

Number of ordered arrangements of r objects chosen from n objects when the n objects are distinct, allowing for repetition

To determine the number of ordered arrangements (Where order matters) out of a group - simply apply the number of objects chosen, as an exponent, to the number of objects in the larger group.

$n^r$

#create combination function
ordered_arrangements <- function (n,r){
  p0 <-  n^r
  p0
}

#number of possible choices
n <- 9
#number of choices used
r <- 7

#call the function
p0 <-  ordered_arrangements(n,r)

# result
result <- paste("The number of possible ",r," digit codes using ",n," possible digits is", p0 ,sep="" )
#result

Result
The number of possible 7 digit codes using 9 possible digits is4782969

Number of permutations of r objects selected from n distinct possibilities without repetition

$P(n,r) = \frac{n!}{(n-r)!}$

#create permutation function
permutations_nonrepeated <- function (n,r){
  p0_top <-  factorial(n)
  p0_bottom <-  factorial(n-r)
  p0 <- p0_top / p0_bottom
}

n <- 26
r <- 4

p0 <- permutations_nonrepeated(n,r)


# result
result <- paste("The number of possible ",r," letter words choosing from ",n," possible letters is ", p0 ," if no letters can be repeated.",sep="" )
#result

Result
The number of possible 4 letter words choosing from 26 possible letters is 358800 if no letters can be repeated.

Tip: The gtools package contains some useful features related to permutations including the permutations function. This function will return a frame of all possible permutations of r objects from n possibilities, with options for repetition.

# all permutations of 2 objects selected from 4 possibilities
n <- 4
r <- 2
perms <- permutations(n = n, 
                      r = r, 
                      repeats.allowed=TRUE)
head(perms)

##      [,1] [,2]
## [1,]    1    1
## [2,]    1    2
## [3,]    1    3
## [4,]    1    4
## [5,]    2    1
## [6,]    2    2

tail(perms)

##       [,1] [,2]
## [11,]    3    3
## [12,]    3    4
## [13,]    4    1
## [14,]    4    2
## [15,]    4    3
## [16,]    4    4

Number of combinations of n distinct objects in groups of r

To find the number of ways that a combination of multiple choices could possibly be arranged from a larger set of possible outcomes, use the following formula:

$C(n,r) = \frac{n!}{(n-r)! \cdot r!}$

#create conbination function
combinations_using_some <- function (n,r){
  p0_top <-  factorial(n)
  p0_bottom <-  factorial(n-r)*factorial(r)
  p0 <- p0_top / p0_bottom
  p0
}

#number of items
n <- 8

#number of items used
r <- 3

p0 <- combinations_using_some(n,r)

# result
result <- paste("Your bike lock has ",r," rotating dials with ",n," letters on each one. The number of ways that you can program this lock is:", p0 ,".",sep="" )

#result

Result
Your bike lock has 3 rotating dials with 8 letters on each one. The number of ways that you can program this lock is:56.

Application: Odds of correctly choosing a Trifecta in a horse race

In the Preakness Stakes, a total of 14 horses race. Imagine that we are betting the Trifecta bet, where we must correctly choose the order of the first 3 horses… How many ways could these 14 horses place in first, second and third place?

You notice that past Trifecta winners have received over 1000:1 payouts - which sounds pretty darn good! But a Trifecta really a good bet? Let’s find out.

This is a fairly simple application of Expected Value - there is only 1 way to win, and every other bet will lose - in other words, if there are $P$ permutations of the first three places then you have $P-1$ chances to lose $1.. and 1 chance to win.

Using the inner product, you could say the expected value is:

$(1*payout) + ((P-1)*cost)$

Where $P$ is the number of permutations of the first 3 places.

Horse Race

num_horses <- 14
num_guessed <- 3
possible_outcomes <- permutations_nonrepeated(num_horses,num_guessed)
odds_correct <- 1/possible_outcomes
chance_correct <- odds_correct * 100

# result
result <- paste("In a horse race with ",num_horses," horses there are ",possible_outcomes," possible ways for the first ",num_guessed," horses to finish - so your chance of guessing exact order of the first ",num_guessed," correctly is : ",round(chance_correct,3),"%. A fair payout for this bet would be: ",possible_outcomes," to 1", sep="" )
#result

Result
In a horse race with 14 horses there are 2184 possible ways for the first 3 horses to finish - so your chance of guessing exact order of the first 3 correctly is : 0.046%. A fair payout for this bet would be: 2184 to 1

Now let’s consider the Kentucky derby where 20 horses run together!

num_horses <- 20
num_guessed <- 3
possible_outcomes <- permutations_nonrepeated(num_horses,num_guessed)
odds_correct <- 1/possible_outcomes
chance_correct <- odds_correct * 100

# result
result <- paste("In a horse race with ",num_horses," horses there are ",possible_outcomes," possible ways for the first ",num_guessed," horses to finish - so your chance of guessing exact order of the first ",num_guessed," correctly is : ",round(chance_correct,3),"%. A fair payout for this bet would be: ",possible_outcomes," to 1", sep="" )
#result

Result
In a horse race with 20 horses there are 6840 possible ways for the first 3 horses to finish - so your chance of guessing exact order of the first 3 correctly is : 0.015%. A fair payout for this bet would be: 6840 to 1

So we can see that even a 2000:1 payout is not especially fair in this case - because there is actually nearly 7000 to 1 odds that you will not choose the outcome correctly.

Application: Choosing a Trifecta (Continued)

Note that the actual payouts of a Trifecta bet are parimutual, like winning a lottery - the amount you win is based on how many $1 bets were placed, divided by how many winners guessed the Trifecta correctly - so you don’t actually know the payout until the race is run.

When the number of horses increases, the odds of choosing the Trifecta correctly goes down - but in this pari-mutual bet - the odds of other betters choosing the same 3 horses as you, also goes down - so presumably, the payout is likely to be higher as well.

This is why choosing less “favored” horses is sometimes popular - because you are betting on the underdogs to make an upset - and presumably it’s less likely that other betters would have selected these horses.

Lets simulate how a real Trifecta bet might work. In this horse race, we will assume just 8 horses are racing - which means we should have a pretty good chance of guessing a trifecta.

num_horses <- 8
num_guessed <- 3
possible_outcomes <- permutations_nonrepeated(num_horses,num_guessed)
odds_correct <- 1/possible_outcomes
chance_correct <- odds_correct * 100

# result
result <- paste("In a horse race with ",num_horses," horses there are ",possible_outcomes," possible ways for the first ",num_guessed," horses to finish - so your chance of guessing exact order of the first ",num_guessed," correctly is : ",round(chance_correct,3),"%. A fair payout for this bet would be: ",possible_outcomes," to 1", sep="" )
#result

Result
In a horse race with 8 horses there are 336 possible ways for the first 3 horses to finish - so your chance of guessing exact order of the first 3 correctly is : 0.298%. A fair payout for this bet would be: 336 to 1

Just to make these horse names memorable - our horses are named after famous music albums of the 60’s, and each horse is handicapped to have different odds of winning.

#assign odds
freak_out <- 3/5
sgt_pepper <- 4/5
electric_lady <- 2/1
pet_sounds <- 4/1
surreal_pillow <- 12/1
cheap_thrills <- 2/5
blues_power <- 1/1
beggars_banquet <-16/1

#save odds and names to vector
horse_handicaps <- c(freak_out,sgt_pepper,electric_lady,pet_sounds,surreal_pillow,cheap_thrills,blues_power,beggars_banquet)
horse_names <- c('freak_out','sgt_pepper','electric_lady','pet_sounds','surreal_pillow','cheap_thrills','blues_power','beggars_banquet')

So we can see that some of the horses like “sgt pepper” are pretty favored and given higher chances to win and pay less when they do. Others like beggars banquet and surreal pillow are considered something of a long shot to win.

Lets imagine that for each person that places a Trifecta bet, they are similarly motivated to win - so they randomly choose three horses with the same likelihood that the are handicapped. This makes it more likely that someone will choose the favored horses.

In this simulation we’re supposing that there are 10,000 bets placed at $1 each. And in an astonishing act of consistency, each better sticks to the same bets for 100 straight races. ALl of this is done so we can determine how our “raw” odds of choosing the right 3 horses compares to our simulated outcomes.

R is particularly useful for simulating random occurrences by using the sample function. This function can take a set of outcomes (like horses crossing the finish line first), and a set of probabilities (like how favored that horse is to win, based on their handicap.)

So we will use the probability as 1 divided by the handicap - exactly the inverse of the handicappers assigned odds.

#set a seed for the following several code blocks to provide consistent results.
set.seed(9)

#loop through 10000 trifecta betters
betters <- c(1:10000)
bets <- data.frame()
for(b in betters){
 thisbet <- sample(x=horse_names, size=3, prob=1/horse_handicaps, replace=FALSE)
 betrow <- data.frame(win=thisbet[1],place=thisbet[2],show=thisbet[3])
 bets <- bets %>% bind_rows(betrow)
}
#peek at some of the bets placed
head(bets,10)

different_bets <- bets %>% summarize(total_distinct = n_distinct(win,place,show)) 
different_bets

We can see that based on this strict interpretation of the handicap - that given 10,000 betters, there are 259 combinations actually selected. A large portion of of these betters are choosing Freak Out, CHeap THrills, and Sgt PEpper to win, show or place - this makes sense as these horses have the best odds… but some people will still choose the underdogs. This is natural human behavior, but there’s also strategy here - if everyone chose the top 3 horses to win - then everyone splits the pot. Betting on the long shots can sometimes pay off better than betting on the favored horses - if you are especially lucky.

Now let’s run 100 simulated races, once again using the handicaps as our vector of possibilities. We are assuming that our handicappers have some kind of supernatural ability and are exceptionally prescient about the odds of each horse to win. Again, this assumption is subject to some variance based on probability, but we will not worry about that wild card here.

In addition - for simplicity, we will also assume that the entire $10,000 race handle that was wagered is paid back out on each race; in reality, the book will keep some portion of the total, similar to the “rake” in poker.

#loop through 10000 trifecta betters
races <- c(1:100)
results <- data.frame()
for(r in races){
 thisresult <- sample(x=horse_names, size=3, prob=1/horse_handicaps, replace=FALSE)
 resultrow <- data.frame(win=thisresult[1],place=thisresult[2],show=thisresult[3])
 results <- results %>% bind_rows(resultrow)
}
#peek at some of the bets placed
head(results,10)

After 100 races - what are our results?

#loop through 100 races
race_num <- 0
outcomes <- data.frame()
races <- c(1:100)

for(x in races){
 race_num <- x
 this_race <- data.frame(results[x,])
 
 this_race_win <- this_race %>% pull(win)
 this_race_place <- this_race %>% pull(place)
 this_race_show <- this_race %>% pull(show)
 
 #count the number of winners in each race
 count_winners <- nrow(bets %>% filter(win==this_race_win & place==this_race_place & show==this_race_show))
 
 this_payout <- 10000 / count_winners
 outcome <- data.frame(race=race_num,
                       winners=count_winners,
                       payout=this_payout,
                       firstplace=this_race_win,
                       secondplace=this_race_place,
                       thirdplace=this_race_show)
 
 outcomes <- outcomes %>% bind_rows(outcome)
}
#peek at some of the bets placed
head(outcomes,10)

Now that we have simulated 100 races - and… we can see why they call it gambling…

outcomes %>% arrange(-payout) %>% slice(1:10)

In these ten races with the fewest winners, it is interesting to note that the riskier bets tended to have long shot horses like surreal pillow, pet sounds, and beggars banquet. Sometimes these horses won, and in some cases, nobody even bet on that combination.

Betting on surreal pillow to win with her 12 to 1 handicap was a long shot - and it paid off with a 10,000 to one payout. But we can also see that in a couple races the winning combination was not chosen by anyone out of 10,000 betters. SO in these cases, nobody wins- except the “book” who keeps all the money that was wagered!

outcomes %>% arrange(payout) %>% slice(1:10)

On the other hand - in these ten races, there were large numbers of winners - since these races all saw the more favored horses like cheap thrills, sgt pepper and freak out winning in varying positions. These were safer bets - and as a result - more people chose them, the bet was split more ways - and thus, they only won payouts between 20 to 1 and 25 to 1.

outcomes %>% group_by(winners,payout,firstplace,secondplace,thirdplace) %>% summarize(occurrences = n_distinct(race)) %>% arrange(-occurrences)

outcomes %>% filter(winners >= 1) %>% summarize(avg_payout = mean(payout), med_payout = median(payout))

Here we can see that the average payout across all the races, (after excluding the ones where there were no winners) - is about $213. So that means that the book is giving (on average) a 213:1 payout vs the 1:336 odds that any one guess will hit the correct Trifecta. Even worse - we can see that the median win is only $54 paid out - so half of the winning bets will return even less than this amount.

Number of permutations of n distinct objects using all n of them:

To find all the permutations of a set of objects when you are using every object in the set, you may simply apply the factorial of the number of objects:

$P(n,n) = n!$

#create permutation function
permutations_using_all <- function (n){
  p0 <-  factorial(n)
  p0
}

#number of items
n <- 6

#call the function
p0 <-  permutations_using_all(n)

# result
result <- paste("The number of arrangements of a ",n," letter word where all the letters are different, using all ",n," letters is ", p0 ,".",sep="" )
#result

Result
The number of arrangements of a 6 letter word where all the letters are different, using all 6 letters is 720.

Application: Rolling the Hard Way in Craps

For a more complex example we will consider the casino dice game of craps. In the game of craps you roll 2 six sided dice, and you may either win an award, win nothing, or lose your bet(s) based on the outcome of the roll. While the game is in the “on” state, you can bet on a specific sum of the 2 dice being rolled. A “hardway” bet means is a bet on a pair of 2’s, 3’s, 4’s or 5’s being rolled. You must roll a hardway before rolling that number the easy way (without a pair) so a bet on the hard 4 (a pair of 2’s) loses if the roller rolls a 3 and a 1 before rolling a 2 and a 2. You also lose your bet if the roller “craps out” by rolling a 7.

Dice

In this sample, the sample space is more complex to determine - we must consider how many ways you can win, and how many ways you can lose. To do so, we will construct a table that represents our sample space with a chance and payout column for each combination of dice roll that is possible.

#what number will you try to roll the hard way?
roll_hard <- 3

#your bet amount
bet <- 2

#how many sided dice?
sides <- 6

#number of dice
numdice <- 2

# build a 36 x 2 frame with all the possible ways to roll these 2 dice
ways_to_roll <- function(n,r){
  ways <- permutations(n = n, r = r, repeats.allowed=TRUE)
  ways
}

#create a function to apply that determines the payout for each outcome
win_lose_draw <- function(v,n,b){
    #the payouts differ depending on the bet
   if(n==3|n==4){
      p <- 9 * b
   }
   if(n==2|n==5){
      p <- 7 * b
   }
   #start from the asumption it's a draw
   payout <- 0
   #if the dice roll is 7, its a crap out
   if(v[1] + v[2] == 7){
     payout <- -b
   }
   #if the dice roll is your numbers the easy way, you lose
   if(v[1] + v[2] == n + n){
     payout <- -b
   }
   #if the dice roll is your numbers the hard way, you win
   if(v[1] == n & v[2] == n){
     payout <- p
   }
   #return payout
   payout
}

#call the ways to roll
dice <- ways_to_roll(sides,numdice)

#now apply the payouts to the ways to roll (returns a vector)
payouts <- apply(dice, 1, win_lose_draw, roll_hard, bet)

#now create a data frame with the roll of d1, d2, and the payout for each roll.
sample_space <- data.frame(d1 = dice[,1], 
                           d2 = dice[,2], 
                           payout = payouts, 
                           chance = 1/nrow(dice))
#view your sample space
sample_space

#calculate expected value
E <- expected_value(sample_space$chance,sample_space$payout)
isfair <- is_fair(E)

#analyze the results
ways_to_lose <- sum(sample_space$payout < 0)
ways_to_win <- sum(sample_space$payout > 0)
ways_to_draw <- sum(sample_space$payout == 0)
max_payout <- max(sample_space$payout)

# result
result <- paste("In the game of craps there are a total of ",nrow(dice)," possible ways that ",numdice," - ",sides,"-sided dice could be thrown. Betting on a hard way, there would be only ",ways_to_win," way to win and ",ways_to_lose," ways to lose, and ",ways_to_draw," possible rolls that would result in a draw. The payout for winning would be $",max_payout," on a bet of only $",bet," but with an expected value of: ",E,", is this a fair bet? - ",isfair,sep="" )

Result
In the game of craps there are a total of 36 possible ways that 2 - 6-sided dice could be thrown. Betting on a hard way, there would be only 1 way to win and 10 ways to lose, and 25 possible rolls that would result in a draw. The payout for winning would be $18 on a bet of only $2 but with an expected value of: -0.0555555555555555, is this a fair bet? - NO

Application (Continued): Simulate from a Sample Set

We can now run simulations in R using the sample function and apply, using our sample space data frame. The sample function enables efficient Bernoulli trials when used with the apply family of functions. All outcomes must have an equal probability for this code to work. There are other ways to use the sample function including the ability to set a vector of probabilities for the sample, which we will explore in a later chapter.

#do a battery of 10 tests
test_list <- as.matrix(seq(1:10))
rolls_per_test <- 1000000

#in each test, run 1,000,000 random dice rolls with the sample space payouts
roll_test <- function(x,ss){
  outcome <- sample(ss$payout,
                    rolls_per_test,
                    replace=TRUE)
  outcome
}

#get the outcomes (LARGE VARIABLE WARNING - DO NOT ATTEMPT TO PRINT/VIEW THIS VAR)
outcomes <- apply(test_list,1,roll_test,sample_space)

#get the sum of each column (10 million total rolls)
list_sums <- colSums(outcomes)
list_actuals <- colSums(outcomes)/rolls_per_test

sum_lost <- sum(list_sums)
low_val <- min(list_actuals)
high_val <- max(list_actuals)

result <- paste("After ",nrow(test_list)," tests of ",format(rolls_per_test,big.mark=',')," rolls, you have lost ",format(sum_lost, nsmall = 0,big.mark=','),". Your tests resulted in losses between ",low_val," and ",high_val," of every $1.",sep="" )

Result
After 10 tests of 1,000,000 rolls, you have lost -566,880. Your tests resulted in losses between -0.06321 and -0.052564 of every $1.

The expected value closely approximates the actual losses in our simulation.

This concludes “A Finite Math Companion in R Markdown” - Thanks for reading, and please stay tuned for more installments:

In “A Probability Companion in R Markdown”, we will explore more fully the subjects of probability, statistics, and additional applications of expected value, including subjects of linear modeling and inference. We will look at how to use R to generate histograms, ogives, box plots, and other statistical charts.

In “A Linear Algebra Companion in R Markdown” we will more thoroughly explore the subject of vector spaces, and how these to visualize vectors geometrically using the principles of Linear Algebra. This guide will explore more about the inner product, number spaces, and the concepts of magnitude and orthogonality.

Appendix: Installing R

R is a free software environment for statistical computing and graphics. It compiles and runs on a wide variety of UNIX platforms, Windows and MacOS. To get the latest version of R, please visit:

https://www.r-project.org/

Packages used in this document are available from the Comprehensive R Archive Network (CRAN):

http://cran.us.r-project.org

For more info on R packages visit:

https://cran.r-project.org/web/packages/

Notes on Publishing RMarkdown

This document was generated by RMarkdown using R code and LaTeX, and rendered to HTML.

Additional formatting has been applied to results variables using the kable and kableExtra packages. These code blocks are not shown in the output, for readability.

A Finite Math Companion in R Markdown

Ryan J. Cooper

3/30/2021

A Finite Math Companion in R Markdown

Foreword

Load packages

Linear Equations

Graphing 2 Points on a Plane

Graphing a Line between 2 Points

Slope of a Line

Graphing a Line with geom_smooth

Graphing a Linear Equation with 3 or more Points with geom_smooth

Graphing a line with arrows with geom_segment

Graphing 2 Lines

Check Intersection of 2 Lines

Application: Rating a Ski Slope

Systems of Linear Equations

Single Linear Expression

Set of Linear Expressions

Matrix Operations

Matrix Equality

Arithmetic Operations on Matrices

Element-Wise Matrix Addition

Element-Wise Matrix Subtraction

Element-Wise Matrix Multiplication

Element-Wise Matrix Division

Scalar Matrix Multiplication

The Dot Product

Matrix Multiplication

Matrices with Unequal Dimensions

Inverse of a Matrix

Application: Cryptography

Gaussian Elimination

Linear Programming with 2 Variables

Application: Donut Shop

Finance

Simple Interest

Simple Interest Formula

Discounted Loans

Amount of a Discounted Loan

Compound Interest

Continuous Compounding Interest

Present Value

Present Value with Continuous Compounding

Amount of an Annuity

Present Value of an Annuity

Amortization

Double an Investment

Time to Reach a Savings Goal with Continuous Compounding

Application: Using sapply to visualize compound interest

Application: Using mapply to plot a ski run by segment

Economics

Compute Inflation

Compute the Consumer Price Index (CPI)

Leontief Models

Application: Building a City on an Island

Sets

Proper Subset / Superset Tests:

Subset / Superset Tests:

Intersection Test

Union Test

Disjoint Test

Creating Sets using seq:

Creating Sets using rnorm:

Application: Random Primes Using the Sieve of Eratosthenes

Number Systems:

Real Numbers

Natural Numbers

Whole Numbers

Integers

Rational Numbers

Descriptive Notation of Sets

Set Related Rules and Formulas

Counting Formula

De Morgan’s Properties

Venn Diagrams

Disjoint Sets

Intersecting Sets

Set with a Subset

Combinations & Permutations