Activity: Probability & Inference
Foundations of inference
Worksheet: Probability & Inference
How to use this worksheet
Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.
This worksheet uses real data from Arctic grayling populations in two Alaskan lakes (I3 and I8), focusing on frequency distributions, how sample size affects our view of a population, and how distributions differ among lakes.
Part 1 · Setup and Data
▶ Run this in your Script:
library(patchwork)
library(skimr)
library(tidyverse)
g_df <- read_csv("data/gray_I3_I8.csv")
i3_df <- g_df %>% filter(lake == "I3")
head(g_df)✏️ Your turn: How many lakes are in this dataset, and how many fish were measured in each? ________________________
Part 2 · Summary Statistics
▶ Run this:
stats_df <- g_df %>%
group_by(lake) %>%
summarize(
mean_length = mean(length_mm, na.rm = TRUE),
sd_length = sd(length_mm, na.rm = TRUE),
se_length = sd(length_mm, na.rm = TRUE) / sum(!is.na(length_mm))^.5,
count = sum(!is.na(length_mm)),
.groups = "drop"
) %>%
arrange(desc(count))
stats_df✏️ Your turn: Which lake has more fish measured? Which has the larger mean length? ________________________
Part 3 · Frequency Distributions
A histogram shows how many observations fall into certain ranges (“bins”).
▶ Run this — a basic histogram for lake I3:
i3_df %>%
ggplot(aes(x = length_mm)) +
geom_histogram(binwidth = 10)🚀 If you finish early: Try changing binwidth to 5 and then to 1. How does the appearance of the histogram change?
# Write your code here:▶ Run this — compare the two lakes:
g_df %>%
ggplot(aes(x = length_mm, fill = lake)) +
geom_histogram(binwidth = 10, position = position_dodge2(width = 0.9))
g_df %>%
ggplot(aes(x = length_mm, fill = lake)) +
geom_histogram(binwidth = 10) +
facet_wrap("lake")✏️ Your turn: Do the two lakes’ distributions look like they overlap a lot, or are they fairly distinct? ________________________
Part 4 · From Histograms to Density Plots
A density plot is a smoothed version of the histogram — the proportion of the data under each part of the curve, which sums to 1 over the whole curve.
▶ Run this:
i3_df %>%
ggplot(aes(x = length_mm)) +
geom_density(fill = "blue", alpha = 0.5)▶ Run this — overlay the density on the histogram:
i3_df %>%
ggplot(aes(x = length_mm)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 2,
fill = "lightblue", alpha = 0.7) +
geom_density(color = "blue", linewidth = 1)Part 5 · Area Under the Density Curve
The whole point of a probability density is that area under it = probability. Let’s confirm the total area is ~1, then find the area for a specific range.
▶ Run this — confirm the total area is ~1:
calculate_density_area <- function(data_vector) {
data_vector <- data_vector[!is.na(data_vector)]
dens <- density(data_vector)
# Trapezoidal rule
dx <- diff(dens$x)
y_avg <- (dens$y[-1] + dens$y[-length(dens$y)]) / 2
area <- sum(dx * y_avg)
return(area)
}
i3_data <- i3_df %>% pull(length_mm)
area_value <- calculate_density_area(i3_data)
i3_df %>%
ggplot(aes(x = length_mm)) +
geom_density(fill = "blue", alpha = 0.4) +
geom_area(stat = "density", fill = "red", alpha = 0.3) +
labs(title = "Area Under Probability Density Function = 1",
subtitle = paste("Calculated area =", round(area_value, 4)),
x = "Length (mm)", y = "Density")▶ Run this — find the probability that a fish falls in a specific length range (you don’t need to understand every line of this function — it’s here to play with):
lower_bound <- 320 # change these
upper_bound <- 350
i3_fish <- i3_df %>%
filter(!is.na(length_mm))
calculate_probability <- function(data_vector, lower_bound, upper_bound) {
dens <- density(data_vector)
indices <- which(dens$x >= lower_bound & dens$x <= upper_bound)
if (length(indices) <= 1) return(0)
x_values <- dens$x[indices]
y_values <- dens$y[indices]
widths <- diff(x_values)
avg_heights <- (y_values[-1] + y_values[-length(y_values)]) / 2
sum(widths * avg_heights)
}
probability <- calculate_probability(i3_fish$length_mm, lower_bound, upper_bound)
total_area <- calculate_probability(i3_fish$length_mm,
min(i3_fish$length_mm),
max(i3_fish$length_mm))
density_data <- density(i3_fish$length_mm)
density_df <- data.frame(x = density_data$x, y = density_data$y)
highlight_df <- density_df %>% filter(x >= lower_bound & x <= upper_bound)
ggplot(i3_fish, aes(x = length_mm)) +
geom_density(fill = "lightblue", alpha = 0.5) +
geom_area(data = highlight_df, aes(x = x, y = y), fill = "darkred", alpha = 0.7) +
geom_vline(xintercept = lower_bound, linetype = "dashed", color = "red") +
geom_vline(xintercept = upper_bound, linetype = "dashed", color = "red") +
labs(
title = "Probability Distribution of Fish Lengths",
subtitle = paste0("Probability of fish between ", lower_bound,
" and ", upper_bound, " mm = ",
round(probability * 100, 1), "%"),
caption = paste("Total area under the curve =", round(total_area, 3)),
x = "Fish Length (mm)", y = "Density"
) +
annotate("text", x = (lower_bound + upper_bound)/2,
y = max(density(i3_fish$length_mm)$y) * 0.7,
label = paste0("Area = ", round(probability, 3)),
color = "white", size = 4) +
theme_minimal() +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "darkred"))✏️ Your turn: Change lower_bound and upper_bound to a range of your choosing. What percentage of fish fall in your new range? ________________________
Part 6 · Z-Scores — a Shortcut for Area
Integrating the area under the curve every time is a pain. Converting to Z-scores gives us a shortcut.
▶ Run this:
mean_length <- mean(i3_df$length_mm, na.rm = TRUE)
sd_length <- sd(i3_df$length_mm, na.rm = TRUE)
i3_df <- i3_df %>%
mutate(z_score = (length_mm - mean_length) / sd_length)
head(i3_df)
z_fish_plot <- i3_df %>%
ggplot(aes(x = z_score)) +
geom_histogram()
z_fish_plot▶ Run this — the proportion within 1 SD, the easy way:
# Proportion within 1 standard deviation = sum of |Z| <= 1, divided by n
within_1sd <- sum(abs(i3_df$z_score) <= 1, na.rm = TRUE) / sum(!is.na(i3_df$z_score))
cat("Proportion within 1 SD:", round(within_1sd * 100, 1), "%\n")💡 In a true normal distribution: ~68% of data within ±1σ, ~95% within ±2σ (really 1.96σ), ~99.7% within ±3σ. This will vary somewhat if the distribution isn’t perfectly normal.
▶ Run this — using R’s built-in normal-distribution functions instead of integrating by hand:
# For the standard normal distribution (mean = 0, sd = 1):
z_value <- 1.22
prob_left <- pnorm(z_value) # area to the left
prob_right <- 1 - pnorm(z_value) # area to the right
prob_between <- pnorm(2) - pnorm(-2) # area between two z-values
z_for_95_percent <- qnorm(0.888) # z-value for a given probability
print(prob_left)
print(prob_right)
print(prob_between)
print(z_for_95_percent)▶ Run this — apply it to a real question: how long does a fish need to be before it’s in the top 5% (unlikely to catch)?
top_5_percent_z <- qnorm(0.95) # z-score for the 95th percentile
unlikely_length <- mean_length + (top_5_percent_z * sd_length)
cat("Only 5% of fish are longer than:", round(unlikely_length, 1), "mm\n")
cat("This corresponds to z-score:", round(top_5_percent_z, 3), "\n")✏️ Your turn: Using pnorm(), what proportion of I3 fish (assuming normality) are longer than 300mm? (Hint: you’ll need to convert 300mm to a z-score first.)
# Write your code here:🚀 If you finish early: Compare within_1sd (computed directly from the data) to pnorm(1) - pnorm(-1) (the theoretical value for a perfect normal distribution). How close are they?
# Write your code here:Part 7 · Comparing a Sample Mean to an Expected Mean
Did this sample come from a population with a specific mean? Let’s practice a one-sample t-test.
▶ Run this — test whether the mean fish length in Lake I3 differs from 260mm:
i3_df <- g_df %>% filter(lake == "I3")
i3_mean <- mean(i3_df$length_mm, na.rm = TRUE)
cat("Mean:", round(i3_mean, 1), "mm\n")
t_test_result <- t.test(i3_df$length_mm, mu = 260)
t_test_result✏️ Your turn: What is the null hypothesis here? Do you reject or fail to reject it at α = 0.05? ________________________
Part 8 · Comparing Two Means
Formulating hypotheses: for the research question “Are fish in Lake I8 longer than fish in Lake I3?”, write the null and alternative hypotheses.
H0 = ________________________________________________________
Ha = ________________________________________________________▶ Run this:
# H0: mu_I3 >= mu_I8, Ha: mu_I3 < mu_I8
t_test_result <- t.test(length_mm ~ lake, data = g_df,
alternative = "less")
t_test_result✏️ Your turn: Based on this t-test, what can you conclude about the difference in fish length between the two lakes? ________________________
🚀 If you finish early: Rerun the test with alternative = "two.sided" instead of "less". Does the p-value change? Does your conclusion change?
# Write your code here:Review and checkpoint
At this point you can:
📤 What to turn in before next class
Upload both of these to the course management system:
- Your code — the
scripts/folder (or just05_probability_inference.R) - This worksheet, with your written answers
Getting unstuck
When code breaks — and it will, that is normal:
- Read the error message out loud. R usually names the line and the problem.
- Check the usual suspects: did you run
library(tidyverse),library(patchwork), andlibrary(skimr)? Spelling? A missing)or%>%at the start of a line? ?function_nameopens the built-in help page.- Bring the exact error (copy-paste it) to class or office hours.
💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.