Introduction to Probability: Completed Tutorial
Sonia Vasa
July 20, 2026
Overview
This tutorial uses empirical probability, the binomial distribution, the geometric distribution, and the normal distribution. The three datasets cover daily PM2.5 estimates across North Carolina, daily rainfall in Chapel Hill, and January maximum temperatures at RDU.
Empirical Probability of Poor Air Quality
Q1. What does each row represent?
Each row represents the PM2.5 value for one date at one geographic
point. I can tell because the same coordinates repeat across dates,
while pm_val changes. One confusing detail is that the
columns are mislabeled: lat contains negative longitude
values and lon contains positive latitude values.
Q2. What needs to happen before calculating probability?
First, each daily observation needs to be classified as either
moderate-or-worse (pm_val > 12) or not. Then the data
need to be grouped by location so the number of flagged days can be
divided by the total number of observed days at that point.
glimpse(pm_nc)## Rows: 548,100
## Columns: 5
## $ year <dbl> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 201…
## $ date <chr> "01JAN2019", "01JAN2019", "01JAN2019", "01JAN2019", "01JAN2019"…
## $ pm_val <dbl> 5.587, 5.731, 5.738, 5.737, 5.770, 5.388, 5.629, 5.863, 4.550, …
## $ lat <dbl> -79.41953, -79.45008, -79.41730, -79.37631, -79.31452, -79.5129…
## $ lon <dbl> 36.09192, 36.04904, 36.05481, 36.06984, 36.04647, 36.13794, 35.…pm_nc |> select(year, date, pm_val, lat, lon) |> slice_head(n = 6)## # A tibble: 6 × 5
## year date pm_val lat lon
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 2019 01JAN2019 5.59 -79.4 36.1
## 2 2019 01JAN2019 5.73 -79.5 36.0
## 3 2019 01JAN2019 5.74 -79.4 36.1
## 4 2019 01JAN2019 5.74 -79.4 36.1
## 5 2019 01JAN2019 5.77 -79.3 36.0
## 6 2019 01JAN2019 5.39 -79.5 36.1Q3. Calculate empirical probability by location
flagged_pm <- pm_nc |>
mutate(aq_flag = if_else(pm_val > 12, 1L, 0L))
grouped_pm <- flagged_pm |>
group_by(lat, lon) |>
summarise(
tot_flag = sum(aq_flag, na.rm = TRUE),
tot_ob = sum(!is.na(pm_val)),
.groups = "drop"
)
pm_prob <- grouped_pm |>
mutate(probability = tot_flag / tot_ob)
pm_prob |>
summarise(
locations = n(),
minimum = min(probability),
median = median(probability),
maximum = max(probability)
) |>
mutate(across(minimum:maximum, scales::percent, accuracy = 0.1)) |>
knitr::kable()| locations | minimum | median | maximum |
|---|---|---|---|
| 300 | 1.9% | 14.8% | 20.7% |
There are 300 locations, each with 1,827 daily observations. Across locations, the estimated probability ranges from about 1.9% to 20.7%, with a median of about 14.8%.
Q4. Convert the results to a spatial object
The tutorial’s coordinate order is correct for this file because the
mislabeled lat column actually contains longitude (x) and
lon contains latitude (y).
grouped_pm_sf <- pm_prob |>
st_as_sf(coords = c("lat", "lon"), crs = 4326, remove = FALSE)
grouped_pm_sf## Simple feature collection with 300 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -84.12346 ymin: 33.91284 xmax: -75.63916 ymax: 36.52073
## Geodetic CRS: WGS 84
## # A tibble: 300 × 6
## lat lon tot_flag tot_ob probability geometry
## * <dbl> <dbl> <int> <int> <dbl> <POINT [°]>
## 1 -84.1 35.2 132 1827 0.0722 (-84.12346 35.19297)
## 2 -83.6 35.2 116 1827 0.0635 (-83.57001 35.20287)
## 3 -83.5 35.5 118 1827 0.0646 (-83.51175 35.53771)
## 4 -83.4 35.3 126 1827 0.0690 (-83.38509 35.28114)
## 5 -83.2 35.1 114 1827 0.0624 (-83.22211 35.09585)
## 6 -83.0 35.3 96 1827 0.0525 (-83.03616 35.27793)
## 7 -82.7 35.5 120 1827 0.0657 (-82.71969 35.54354)
## 8 -82.7 35.2 124 1827 0.0679 (-82.65454 35.23311)
## 9 -82.6 35.6 135 1827 0.0739 (-82.60949 35.5707)
## 10 -82.6 35.7 128 1827 0.0701 (-82.57899 35.7072)
## # ℹ 290 more rowsQ5. Histogram and map
ggplot(pm_prob, aes(x = probability)) +
geom_histogram(bins = 15, fill = "#2563EB", color = "white") +
scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(
title = "Distribution of Moderate-or-Worse PM2.5 Probability",
x = "Empirical probability at monitoring location",
y = "Number of locations"
)
nc_state <- tryCatch(
tigris::states(cb = TRUE, year = 2024, class = "sf") |>
filter(STUSPS == "NC") |>
st_transform(4326),
error = function(e) NULL
)## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%pm_map <- ggplot()
if (!is.null(nc_state)) {
pm_map <- pm_map +
geom_sf(data = nc_state, fill = "grey98", color = "grey35", linewidth = 0.5)
}
pm_map +
geom_sf(
data = grouped_pm_sf,
aes(color = probability),
size = 1.8,
alpha = 0.9
) +
scale_color_viridis_c(
labels = scales::percent_format(accuracy = 1),
limits = c(0, max(pm_prob$probability)),
name = "Probability\nPM2.5 > 12"
) +
coord_sf(xlim = c(-84.6, -75.2), ylim = c(33.7, 36.8), expand = FALSE) +
labs(
title = "Empirical Probability of Moderate-or-Worse PM2.5 Days",
x = "Longitude",
y = "Latitude"
)
The map shows the highest estimated probabilities in the central and western-central part of the state. The eastern and far-western locations generally have lower probabilities.
Binomial Distribution: Chapel Hill Rainfall
ch_rainfall_binary <- ch_rainfall |>
mutate(rain = if_else(prec_in > 0, 1L, 0L))
total_days <- sum(!is.na(ch_rainfall_binary$prec_in))
total_rainfall_days <- sum(ch_rainfall_binary$rain, na.rm = TRUE)
p_rain <- total_rainfall_days / total_days
probability_days <- tibble(
num_days = 0:10,
probability = dbinom(0:10, size = 10, prob = p_rain)
)
probability_days |>
ggplot(aes(x = num_days, y = probability)) +
geom_line(color = "#0F766E", linewidth = 1) +
geom_point(color = "#0F766E", size = 2.5) +
scale_x_continuous(breaks = 0:10) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(
title = "Binomial Model for Rainy Days in a 10-Day Period",
x = "Number of rainy days",
y = "Probability"
)
Q6. Interpret the binomial distribution
most_likely_rain_days <- probability_days$num_days[which.max(probability_days$probability)]
probability_zero_rain_days <- dbinom(0, size = 10, prob = p_rain)
tibble(
result = c(
"Empirical probability of rain on a day",
"Most likely number of rainy days out of 10",
"Probability of zero rainy days out of 10"
),
value = c(
scales::percent(p_rain, accuracy = 0.1),
as.character(most_likely_rain_days),
scales::percent(probability_zero_rain_days, accuracy = 0.1)
)
) |>
knitr::kable()| result | value |
|---|---|
| Empirical probability of rain on a day | 33.7% |
| Most likely number of rainy days out of 10 | 3 |
| Probability of zero rainy days out of 10 | 1.6% |
The most likely outcome is 3 rainy days in a 10-day period. The chance of getting no rain at all during the 10 days is only about 1.6%.
The cumulative probability of three rainy days or fewer and its complement are:
less_or_equal_3_days <- pbinom(3, size = 10, prob = p_rain)
more_than_3_days <- 1 - less_or_equal_3_days
tibble(
outcome = c("3 rainy days or fewer", "More than 3 rainy days"),
probability = scales::percent(
c(less_or_equal_3_days, more_than_3_days),
accuracy = 0.1
)
) |>
knitr::kable()| outcome | probability |
|---|---|
| 3 rainy days or fewer | 54.9% |
| More than 3 rainy days | 45.1% |
Mini Challenge: Rain Over 0.5 Inches
ch_heavy_rain <- ch_rainfall |>
mutate(heavy_rain = if_else(prec_in > 0.5, 1L, 0L))
p_heavy_rain <- mean(ch_heavy_rain$heavy_rain, na.rm = TRUE)
heavy_rain_days <- tibble(
num_days = 0:10,
probability = dbinom(0:10, size = 10, prob = p_heavy_rain)
)
most_likely_heavy_days <- heavy_rain_days$num_days[
which.max(heavy_rain_days$probability)
]
prob_at_least_3_heavy <- 1 - pbinom(
2,
size = 10,
prob = p_heavy_rain
)
heavy_rain_days |>
ggplot(aes(x = num_days, y = probability)) +
geom_col(fill = "#EA580C", color = "white") +
scale_x_continuous(breaks = 0:10) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(
title = "Days With More Than 0.5 Inches of Rain",
x = "Number of heavy-rain days in 10 days",
y = "Probability"
)
tibble(
result = c(
"Empirical probability of more than 0.5 inches on a day",
"Most likely number of heavy-rain days out of 10",
"Probability of at least 3 heavy-rain days out of 10"
),
value = c(
scales::percent(p_heavy_rain, accuracy = 0.1),
as.character(most_likely_heavy_days),
scales::percent(prob_at_least_3_heavy, accuracy = 0.1)
)
) |>
knitr::kable()| result | value |
|---|---|
| Empirical probability of more than 0.5 inches on a day | 9.1% |
| Most likely number of heavy-rain days out of 10 | 1 |
| Probability of at least 3 heavy-rain days out of 10 | 5.6% |
The most likely number is 1 day with more than 0.5 inches of rain. The chance of at least 3 such days in a 10-day period is about 5.6%.
Both assigned questions count events in a fixed 10-day period, so the binomial distribution is the correct model. A geometric distribution instead describes the wait until the first heavy-rain day. As an extra check, it estimates an average wait of about 11 days and a 61.6% chance that the first heavy-rain day occurs within the next 10 days.
expected_wait_days <- 1 / p_heavy_rain
prob_first_heavy_within_10 <- pgeom(9, prob = p_heavy_rain)
tibble(
geometric_result = c(
"Expected days through first heavy-rain day",
"Probability first heavy-rain day occurs within 10 days"
),
value = c(
round(expected_wait_days, 1),
scales::percent(prob_first_heavy_within_10, accuracy = 0.1)
)
) |>
knitr::kable()| geometric_result | value |
|---|---|
| Expected days through first heavy-rain day | 11 |
| Probability first heavy-rain day occurs within 10 days | 61.6% |
Normal Distribution: January Temperatures at RDU
rdu_jan_complete <- rdu_jan |>
filter(!is.na(max_temp))
jan_mean <- mean(rdu_jan_complete$max_temp)
jan_sd <- sd(rdu_jan_complete$max_temp)
jan_norm <- tibble(
temp = seq(
min(rdu_jan_complete$max_temp) - 3,
max(rdu_jan_complete$max_temp) + 3,
length.out = 400
),
density = dnorm(temp, mean = jan_mean, sd = jan_sd)
)
ggplot(rdu_jan_complete, aes(x = max_temp)) +
geom_histogram(
aes(y = after_stat(density)),
bins = 16,
fill = "#BFDBFE",
color = "white"
) +
geom_line(
data = jan_norm,
aes(x = temp, y = density),
inherit.aes = FALSE,
color = "#2563EB",
linewidth = 1
) +
labs(
title = "January Maximum Temperatures at RDU",
x = "Maximum temperature (F)",
y = "Density"
)
The center of the observed distribution is fairly close to a bell shape, but the fitted curve does not capture every part of the distribution equally well.
qqnorm(
rdu_jan_complete$max_temp,
pch = 19,
col = "#2563EB",
main = "Normal Q-Q Plot for January Maximum Temperature"
)
qqline(rdu_jan_complete$max_temp, col = "#EA580C", lwd = 2)
The middle of the Q-Q plot tracks the reference line reasonably well, while both tails bend away from it. That matters most when the normal model is used for rare temperatures.
As a check on the tutorial example, the normal model estimates a 22.4% probability of a January maximum at or above 60 F.
prob_60_normal <- 1 - pnorm(60, mean = jan_mean, sd = jan_sd)
scales::percent(prob_60_normal, accuracy = 0.1)## [1] "22.4%"Mini Challenge: Maximum Temperature of 83 F or Above
prob_83_normal <- 1 - pnorm(83, mean = jan_mean, sd = jan_sd)
prob_83_empirical <- mean(rdu_jan_complete$max_temp >= 83)
temperature_summary <- tibble(
measure = c(
"Number of January observations",
"Mean January maximum",
"Standard deviation",
"Highest observed January maximum",
"Normal-model probability of 83 F or above",
"Empirical probability of 83 F or above"
),
value = c(
nrow(rdu_jan_complete),
paste0(round(jan_mean, 1), " F"),
paste0(round(jan_sd, 1), " F"),
paste0(max(rdu_jan_complete$max_temp), " F"),
scales::percent(prob_83_normal, accuracy = 0.001),
scales::percent(prob_83_empirical, accuracy = 0.001)
)
)
temperature_summary |>
knitr::kable()| measure | value |
|---|---|
| Number of January observations | 619 |
| Mean January maximum | 51.2 F |
| Standard deviation | 11.6 F |
| Highest observed January maximum | 78 F |
| Normal-model probability of 83 F or above | 0.300% |
| Empirical probability of 83 F or above | 0.000% |
tibble(
method = c("Normal model", "Observed data"),
probability = c(prob_83_normal, prob_83_empirical)
) |>
ggplot(aes(x = method, y = probability, fill = method)) +
geom_col(width = 0.6, show.legend = FALSE) +
geom_text(
aes(label = scales::percent(probability, accuracy = 0.001)),
vjust = -0.5,
fontface = "bold"
) +
scale_fill_manual(values = c("#2563EB", "#EA580C")) +
scale_y_continuous(
labels = scales::percent_format(accuracy = 0.1),
limits = c(0, 0.004)
) +
labs(
title = "Probability of a January Maximum at or Above 83 F",
x = NULL,
y = "Probability"
)
Discussion
The normal model gives a probability of about 0.300%, while the empirical probability is 0% because none of the 619 observed January days reached 83 F; the highest observed value was 78 F.
The estimates differ because the normal curve is continuous and extends beyond the observed range, so it always assigns some probability to extremely high values. The empirical estimate only counts what actually occurred in this sample. With a rare event, even hundreds of observations may not include a single occurrence.
The main problem with the normal model here is that rare-event estimates depend on the tails, and the Q-Q plot shows that the tails are not very normal. The model can therefore overestimate or underestimate extreme-temperature risk. Temperature data can also change over time and may not be identically distributed across years, so the fitted historical average may not represent current or future January conditions.
Conclusion
The empirical calculations summarize what happened in the observed data, while the theoretical distributions extend those patterns to events that may not have appeared in the sample. The binomial model gives useful short-term rainfall probabilities, but its independence and constant-probability assumptions are simplifications. The normal model works reasonably near the middle of the January temperature distribution, but its estimate for 83 F should be treated cautiously because it relies entirely on the fitted tail.