Week 2 Discussion: Writing Functions in R and Bayes’ Theorem
I. Classes and data structures in R
A class describes what kind of object R is working with and helps determine how functions treat that object. For example, TRUE has class logical, 3L has class integer, 3.5 has class numeric, and flower has class character. Classes also describe more specialized objects, such as a factor for categories or a data.frame for a dataset. I think of class as the object’s meaning to R while typeof() identifies its underlying storage type. These descriptions can differ: a number such as 3.5 has class numeric but type double.
A data structure describes how values are organized. An atomic vector, such as c(2, 4, 6) holds values of one type. A matrix arranges values of one type into rows and columns. A list can contain different kinds of objects, such as a number, a character vector, and a matrix. A data frame organizes observations into rows and variables into columns, allowing different columns to have different types. A factor represents categories using integer codes and level labels, while a table summarizes counts, such as the number of flowers in each species. Classes and structures overlap: a data frame is a structure that also has class data.frame.
I chose the iris dataset which is included with R and contains 150 flowers, four measurements in centimeters, and a species column.
# Load the built-in dataset; no external download is needed.
data("iris", package = "datasets")
# Compare the high-level class with the underlying storage type.
class(iris)## [1] "data.frame"typeof(iris)## [1] "list"class(iris$Sepal.Length)## [1] "numeric"typeof(iris$Sepal.Length)## [1] "double"class(iris$Species)## [1] "factor"typeof(iris$Species)## [1] "integer"# Inspect the dimensions, column classes, and example values.
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...The dataset has class data.frame and type list. Sepal.Length has class numeric and type double, so R stores its measurements as double-precision numbers. Species has class factor and type integer because its category labels are represented internally by integer codes. The str() output confirms 150 observations, four numeric columns, and one factor with three levels. class() and typeof() give complementary information rather than contradictory answers.
II. Reading and writing functions
To examine existing functions I uaws sd() and IQR(). R is case-sensitive, so the second function is spelled IQR not iqr. Referring to a function without calling it exposes its definition; the chunk below prints those definitions as readable text.
# Display the function definitions rather than evaluate a statistic.
cat(deparse(stats::sd), sep = "\n")## function (x, na.rm = FALSE)
## sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x),
## na.rm = na.rm))cat("\n\n")cat(deparse(stats::IQR), sep = "\n")## function (x, na.rm = FALSE, type = 7)
## diff(quantile(as.numeric(x), c(0.25, 0.75), na.rm = na.rm, names = FALSE,
## type = type))In sd(), the expression sqrt(var(…)) suggests that R first calculates variance and then takes its square root. The result is the sample standard deviation, which measures spread around the mean in the original units. The variance calculation uses the denominator \(n-1\), giving \(s=\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2/(n-1)}\).
In IQR() quantile(…, c(0.25, 0.75)) finds the first and third quartiles and diff() subtracts the first from the third. The calculation is \(IQR=Q_3-Q_1\) the width of the middle 50% of the data. Its default type = 7 specifies how R calculates the quantiles. Both functions are appropriate for a numeric vector such as iris$Sepal.Length species labels would not be a meaningful input for these calculations.
The next chunk applies both functions to sepal length and also displays the quartiles so that the IQR calculation can be checked directly.
# Select the numeric vector of sepal lengths, measured in centimeters.
sepal_length <- iris$Sepal.Length
# Calculate the sample std and iqr.
sd(sepal_length)## [1] 0.8280661quantile(sepal_length, probs = c(0.25, 0.75))## 25% 75%
## 5.1 6.4IQR(sepal_length)## [1] 1.3The sample standard deviation is approximately 0.8281 cm. The quartiles are 5.1 cm and 6.4 cm, so the IQR is 1.3 cm.
I calculated the arithmetic mean by adding the values and dividing by their count. The function accepts a numeric vector and optionally removes missing values. The default na.rm = FALSE a missing value makes the result missing. na.rm = TRUE the denominator counts only the values retained. If no values remain the calculation returns NaN because an empty vector has no arithmetic mean.
# Calculate an arithmetic mean without using mean()
my_mean <- function(x, na.rm = FALSE) {
# Require a numeric vector rather than a data frame or matrix.
if (!is.numeric(x) || !is.null(dim(x))) {
stop("x must be a numeric vector.")
}
# Remove missing values only when requested by the caller.
if (na.rm) {
x <- x[!is.na(x)]
}
# Divide the sum by the number of values actually used.
return(sum(x) / length(x))
}
# Compare the custom result built-in mean on the same data.
my_mean(iris$Sepal.Length)## [1] 5.843333mean(iris$Sepal.Length)## [1] 5.843333isTRUE(all.equal(my_mean(iris$Sepal.Length), mean(iris$Sepal.Length)))## [1] TRUE# Demonstrate how the missing-value option changes the calculation.
my_mean(c(2, NA_real_, 4))## [1] NAmy_mean(c(2, NA_real_, 4), na.rm = TRUE)## [1] 3Both mean functions return approximately 5.8433 cm and the comparison returns TRUE. For the small missing-value example the default result is NA, while removing the missing value gives \((2+4)/2=3\).
I also tried the supplied myglimpse() example on iris$Species. This function counts missing values separately, builds a frequency table, adds its total, and combines the results into one named vector.
# Reproduce the supplied function with comments explaining each step.
myglimpse <- function(x) {
a <- sum(is.na(x)) # Count missing observations.
b <- addmargins(table(x)) # Count categories and add their sum.
temp <- c(a, b) # Combine missing and category counts.
names(temp) <- c("Missing", unique(names(b)))
return(temp) # Return the named summary vector.
}
# Summarize the categorical species column, not the entire data frame.
myglimpse(iris$Species)## Missing setosa versicolor virginica Sum
## 0 50 50 50 150There are 0 missing values, 50 flowers in each species, and 150 nonmissing observations in total. table() excludes missing values by default, the Sum entry counts nonmissing observations; Missing reports any excluded observations separately.
III. Bayes theorem
Bayes theorem updates the probability of an event after new evidence is observed. \(P(A)\) is the prior probability. \(P(B\mid A)\) describes how likely the evidence is if \(A\) occurs and \(P(A\mid B)\) is the updated probability with \(P(B)>0\).
Suppose Box A and Box B are equally likely to be selected and their red-ball proportions are 80% and 20% respectively. If a ball drawn at random from the selected box is red, the probability that Box A was selected is \((0.80\times0.50)/(0.80\times0.50+0.20\times0.50)=0.80\), or 80%. Observing red therefore raises the probability of Box A from 50% to 80%.
\[ P(A\mid B)=\frac{P(B\mid A)P(A)}{P(B)} =\frac{P(B\mid A)P(A)}{P(B\mid A)P(A)+P(B\mid A^c)P(A^c)}. \]
$$
P(A\mid B)=\frac{P(B\mid A)P(A)}{P(B)}
$$IV. The full parking garage
Let \(F\) mean that the garage is full. First, add the probabilities of all three ways the garage can be full:
\[ P(F)=P(A)P(F\mid A)+P(S)P(F\mid S)+P(N)P(F\mid N) =0.35(0.25)+0.20(0.70)+0.45(0.05)=0.25. \]
Apply Bayes theorem:
\[ P(S\mid F)=\frac{P(F\mid S)P(S)}{P(F)} =\frac{0.70(0.20)}{0.25}=0.56=56\%. \]
Knowing the garage is full raises the probability of a sporting event from 20% to 56%. Among full-garage evenings, 56% involve a sporting event. This differs from the given 70% which describes the probability of a full garage given a sporting event which reverses the condition changes the probability.
The next chunk calculates each joint probability by multiplying along its tree path, and divides the sporting-and-full probability by the total probability of a full garage. It also checks that the six possible outcomes account for all evenings.
# Enter the probabilities from Guided Practice 3.43.
event_names <- c("Academic event", "Sporting event", "No event")
prior <- c(0.35, 0.20, 0.45)
full_given_event <- c(0.25, 0.70, 0.05)
# The second branch from each event is the complementary outcome.
spaces_given_event <- 1 - full_given_event
# Multiply each prior by its conditional probability to obtain a leaf.
joint_full <- prior * full_given_event
joint_spaces <- prior * spaces_given_event
# Sum all full-garage leaves, then apply Bayes' theorem.
prob_full <- sum(joint_full)
sporting_index <- match("Sporting event", event_names)
prob_sporting_given_full <- joint_full[sporting_index] / prob_full
# Show the branch probabilities and the two leaves for each event.
tree_probabilities <- data.frame(
Event = event_names,
Prior = prior,
Full_given_event = full_given_event,
Spaces_given_event = spaces_given_event,
Event_and_full = joint_full,
Event_and_spaces = joint_spaces
)
print(tree_probabilities, row.names = FALSE)## Event Prior Full_given_event Spaces_given_event Event_and_full
## Academic event 0.35 0.25 0.75 0.0875
## Sporting event 0.20 0.70 0.30 0.1400
## No event 0.45 0.05 0.95 0.0225
## Event_and_spaces
## 0.2625
## 0.0600
## 0.4275print(c(P_full = prob_full, P_sporting_given_full = prob_sporting_given_full))## P_full P_sporting_given_full
## 0.25 0.56# Check that the priors and the six joint probabilities each sum to 1.
stopifnot(
isTRUE(all.equal(sum(prior), 1)),
isTRUE(all.equal(sum(joint_full) + sum(joint_spaces), 1))
)The tree below starts with the parking garage, branches into the
three event types, and then shows whether spaces are full or available.
Each leaf includes its joint probability. I use the diagram
package to draw a simple black-and-white tree.
# Install diagram once if needed: install.packages("diagram")
library(diagram)
# Label the root, event types, and six outcomes.
events <- c("Academic", "Sporting", "No Event")
leaf_names <- paste(rep(events, each = 2), c("Full", "Available"))
leaf_prob <- as.vector(rbind(joint_full, joint_spaces))
node_names <- c("Parking Garage", events,
paste(leaf_names, sprintf("%.2f%%", 100 * leaf_prob), sep = "\n"))
# A matrix entry connects a parent column to a child row.
branches <- matrix("", nrow = 10, ncol = 10)
branches[2:4, 1] <- paste0(100 * prior, "%")
for (i in 1:3) {
branches[(2 * i + 3):(2 * i + 4), i + 1] <-
paste0(100 * c(full_given_event[i], spaces_given_event[i]), "%")
}
# Draw one root, three event nodes, and six outcome nodes.
draw_probability_tree <- function() {
par(mar = c(0, 0, 0, 0))
plotmat(branches, pos = c(1, 3, 6), name = node_names,
absent = "", curve = 0, latex = TRUE, box.type = "ellipse",
box.size = 0.075, box.prop = 0.30, box.cex = 1.0,
box.col = "white", shadow.size = 0, lwd = 1,
arr.type = "triangle", arr.pos = 0.62,
arr.length = 0.12, cex.txt = 1.0)
}The following chunk displays the tree and saves it as
bayes_tree.png for attachment.
# Show the tree in the report.
invisible(draw_probability_tree())
Probability tree for the parking garage. Leaf percentages are joint probabilities.
# Save the same tree as a high-resolution image.
png("bayes_tree.png", width = 2400, height = 800, res = 200)
invisible(draw_probability_tree())
invisible(dev.off())