Bayes Rule Notebook
QSCI 497, April 6, 2026
1. Review / Helper Material
Bayes Rule is defined as \[ P(A|B) = \frac{P(B|A) \times P(A)}{P(B)} \] Where the notation \(P(A)\) refers to the probability of A, and \(P(A|B)\) refers to the probability of A given B – the “|” is pronounced as “given” and denotes a conditional probability.
Each term in Bayes Rule is defined:
\(P(A)\) - prior probability, representing our current knowledge
\(P(B|A)\) - likelihood, representing some observation or prediction
\(P(A|B)\) - posterior probability, representing our knowledge as improved by some observation
\(P(B)\) is the probability of the
evidence. For today, it is assumed to be a constant and will be ignored.
Therefore Bayes rule is simplified to: \[
P(A|B) \propto P(B|A) \times P(A)
\]
Normalization
It is necessary to normalize all probabilities such that \[ \int P(x) dx = 1 \] In other words, the sum of probabilities of all possible outcomes must be equal to 1.
2. Prior PDF
In class I used the historic record of observed hurricane tracks to propose a prior PDF:
I used a Gaussian probability distribution, defined as \[ P(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \] Where \(\mu\) = 1200 km and \(\sigma\) = 1500 km. With some location references, that is centered near New Orleans. Let’s plot that now:
# Define the Gaussian function
gaussian <- function(x, mean, stdev) {
(1 / (stdev * sqrt(2 * pi))) * exp(-((x - mean)^2) / (2 * stdev^2))
}
# Sequence of x values
x <- seq(-2, 5.5, by = 0.1)
# Calculate Gaussian (prior)
prior <- gaussian(x, mean = 1.2, stdev = 1.5)
# Alternatively: prior <- dnorm(x, mean = 1.2, sd = 1.5)
# Plot
plot(x, prior, type = "l", col = "blue",
xlab = "Distance Along-Coast [thousand km]",
ylab = "Probability Density [1/thousand km]")3. Forecast 1: September 28, 2022
Given this pressure map, what is the most likely path for a hurricane to take?
(Plot is the height of the 500 mb contour (which is on average around 5 km high, a representative area for hurricane steering winds) – high pressure areas means the contour is elevated, and vice versa.)
Interactive forecast diagram: https://scied.ucar.edu/interactive/forecast-hurricane
Stop: pick a landfall spot first…
# Define Gaussian function
gaussian <- function(x, mean, stdev) {
(1 / (stdev * sqrt(2 * pi))) * exp(-((x - mean)^2) / (2 * stdev^2))
}
# x values
x <- seq(-2, 5.5, by = 0.1)
# Prior distribution
prior <- gaussian(x, mean = 1.2, stdev = 1.5)
# Likelihood distribution
landfall_spot <- 0 # Replace 0 with your actual value
landfall_deviation <- 0.5
likelihood <- gaussian(x, mean = landfall_spot, stdev = landfall_deviation)
# Plot both distributions
plot(x, prior, type = "l", col = "blue", lwd = 2,
xlab = "Distance Along-Coast [thousand km]",
ylab = "Probability Density [1/thousand km]",
ylim = range(prior, likelihood))
lines(x, likelihood, col = "red", lwd = 2)
legend("topright", legend = c("Prior", "Likelihood"),
col = c("blue", "red"), lwd = 2)Computing the Posterior
The posterior is the product of the two PDFs, normalized to 1. Let’s start with the product:
posterior <- prior * likelihoodTo normalized it, we need to make sure the integral adds up to 1.
# If not already installed:
# install.packages("pracma")
library(pracma)
# Trapezoidal integration
posterior_integral <- trapz(x, posterior)
# Normalized posterior
posterior_normed <- posterior / posterior_integralNow let’s plot it!
# Plot the prior first
plot(x, prior, type = "l", col = "blue", lwd = 2,
xlab = "Distance Along-Coast [thousand km]",
ylab = "Probability Density [1/thousand km]",
ylim = range(prior, likelihood, posterior_normed))
# Add likelihood and posterior
lines(x, likelihood, col = "red", lwd = 2)
lines(x, posterior_normed, col = "green", lwd = 2)
# Add a legend
legend("topright", legend = c("Prior", "Likelihood", "Posterior"),
col = c("blue", "red", "green"), lwd = 2)What is the maximum of the posterior? This is the most likely landfall location…
# Get the index of the maximum value in the posterior
max_index <- which.max(posterior)
# Find the corresponding distance along the coast
x[max_index]What city does this correspond to?
4. Forecast 2: October 27, 2012
Now on your own! With the following pressure map:
# Define Gaussian function
gaussian <- function(x, mean, stdev) {
(1 / (stdev * sqrt(2 * pi))) * exp(-((x - mean)^2) / (2 * stdev^2))
}
# x values
x <- seq(-2, 5.5, by = 0.1)
# Prior distribution
prior <- gaussian(x, mean = 1.2, stdev = 1.5)
# Likelihood distribution
landfall_spot <- 0 # Replace 0 with your actual value
landfall_deviation <- 0.5
likelihood <- gaussian(x, mean = landfall_spot, stdev = landfall_deviation)
# Plot prior and likelihood
plot(x, prior, type = "l", col = "blue", lwd = 2,
xlab = "Distance Along-Coast [thousand km]",
ylab = "Probability Density [1/thousand km]",
ylim = range(prior, likelihood))
lines(x, likelihood, col = "red", lwd = 2)
# Add legend
legend("topright", legend = c("Prior", "Likelihood"),
col = c("blue", "red"), lwd = 2)# Compute unnormalized posterior
posterior <- prior * likelihood
# If not already installed:
# install.packages("pracma")
library(pracma)
# Trapezoidal integration
posterior_integral <- trapz(x, posterior)
# Normalized posterior
posterior_normed <- posterior / posterior_integral# Plot the prior first
plot(x, prior, type = "l", col = "blue", lwd = 2,
xlab = "Distance Along-Coast [thousand km]",
ylab = "Probability Density [1/thousand km]",
ylim = range(prior, likelihood, posterior_normed))
# Add likelihood and posterior
lines(x, likelihood, col = "red", lwd = 2)
lines(x, posterior_normed, col = "green", lwd = 2)
# Add a legend
legend("topright", legend = c("Prior", "Likelihood", "Posterior"),
col = c("blue", "red", "green"), lwd = 2)# Get the index of the maximum value in the posterior
max_index <- which.max(posterior)
# Find the corresponding distance along the coast
x[max_index]5. Bonus Section
The technique you just used is called Maximum a Posteriori (or Bayesian MAP) – go on google scholar and search “maximum a posteriori + your subject of interest” and see what comes up.