Homework 3

Homework 3#

Coverage: Lectures 11–12
Due: Sunday, September 20, 2026, 11:59 p.m. ET
Total: 100 points

What this homework practices#

Lecture

Assessed ideas

11 — Selecting prior information

Translate support, mean, and variance information into a maximum-entropy Gaussian prior; distinguish prior location from prior concentration; test sensitivity to alternative prior information.

12 — Analytical Bayesian inference

Derive an analytical posterior, quantify posterior uncertainty, design a sample size, form posterior-predictive distributions, simulate replicated data, and make a probability-based decision.

Useful course-book pages are Lecture 11: selecting prior information, continuous maximum-entropy examples, Lecture 12: analytical Bayesian inference, credible intervals, and posterior-predictive distributions.

Normal–Normal and Gamma–Poisson conjugacy are used as guided extensions of Lecture 12. Every required definition and parameterization is supplied below.

Instructions#

  • Complete this notebook in Google Colab or another fresh Python environment.

  • Problem 1 is a short hand-calculation problem. Show the important derivation and interpretation steps in Markdown/LaTeX, or insert one clearly legible image of your handwritten work. If you insert an image, accompany it with a descriptive alternative-text summary in the response cell and type all final mathematical results there. No code is needed.

  • Problem 2 is a scaffolded scientific-computing study. Use the supplied random seeds and do not delete setup, helper, or data-credit cells.

  • The setup cell downloads and verifies the frozen course data automatically. Do not upload or replace the CSV, and do not query the live USGS catalog.

  • Do not mount Google Drive or use Drive/absolute data paths. Saving your Colab notebook copy in Drive is fine. Do not install additional packages.

  • Label plots and include units. Report numerical answers to at least four significant digits unless instructed otherwise.

  • Follow the current course and Brightspace submission policies.

Student details#

  • First name:

  • Last name:

  • Purdue email:

from io import BytesIO
from pathlib import Path
import hashlib
import urllib.error
import urllib.request

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

SEED = 53903
sns.set_theme(style="ticks", context="notebook")
plt.rcParams["figure.dpi"] = 120
plt.rcParams["figure.constrained_layout.use"] = True
np.set_printoptions(precision=6, suppress=True)

DATA_NAME = "southern-california-baja-earthquakes-1900-2025.csv"
EXPECTED_DATA_SHA256 = (
    "a4e7fb0d2a99fc81cf0f309e13a853b6025413e84d61ab82f18cc8c43ec37072"
)
COURSE_DATA_URL = (
    "https://predictivesciencelab.github.io/data-analytics-se/_downloads/"
    "0df04c34eed724ea5eaaf5e918b4496c/"
    "southern-california-baja-earthquakes-1900-2025.csv"
)
LOCAL_DATA_CANDIDATES = (
    Path(DATA_NAME),
    Path("../data/homework") / DATA_NAME,
    Path("lecturebook/data/homework") / DATA_NAME,
)


def sha256_bytes(payload):
    '''Return the hexadecimal SHA-256 digest of a byte string.'''
    return hashlib.sha256(payload).hexdigest()


def load_frozen_course_data():
    '''Load the hash-verified snapshot locally or from the course repository.'''
    rejected = []
    for candidate in LOCAL_DATA_CANDIDATES:
        if candidate.is_file():
            payload = candidate.read_bytes()
            digest = sha256_bytes(payload)
            if digest == EXPECTED_DATA_SHA256:
                return payload, f"verified local course file: {candidate}"
            rejected.append(f"{candidate} (SHA-256 {digest})")

    request = urllib.request.Request(
        COURSE_DATA_URL,
        headers={"User-Agent": "ME539-course-data/2026"},
    )
    try:
        with urllib.request.urlopen(request, timeout=30) as response:
            payload = response.read()
    except (urllib.error.URLError, TimeoutError) as exc:
        rejected_text = (
            f" Rejected local files: {', '.join(rejected)}."
            if rejected
            else ""
        )
        raise RuntimeError(
            "Could not download the frozen Homework 3 course data. "
            "Check the internet connection and rerun this cell."
            + rejected_text
        ) from exc

    digest = sha256_bytes(payload)
    if digest != EXPECTED_DATA_SHA256:
        raise ValueError(
            "Downloaded bytes do not match the frozen course snapshot: "
            f"expected {EXPECTED_DATA_SHA256}, received {digest}."
        )
    return payload, f"verified published course file: {COURSE_DATA_URL}"


DATA_BYTES, DATA_SOURCE = load_frozen_course_data()
print(f"Loaded {len(DATA_BYTES)} verified bytes from {DATA_SOURCE}.")

Problem 1 — Maximum entropy and Bayesian sensor calibration (25 points)#

A sensor’s calibration residual is its reported value minus a known reference value. Let \(b\), measured in millivolts (mV), denote the sensor’s unknown but fixed signed zero offset. Before collecting calibration measurements, the available engineering information is

\[ b\in\mathbb{R}, \qquad \mathbb{E}[b]=b_0, \qquad \mathbb{V}[b]=\tau_0^2. \]

Assume \(\tau_0>0\).

Support and moments do not uniquely determine a prior distribution. Once the variable and reference density are specified, the maximum-entropy principle selects the feasible distribution that introduces the least additional structure beyond the stated information.

1.1 Derive the maximum-entropy prior (7 points)#

Take the reference density to be \(q(b)=1\) relative to Lebesgue measure—that is, ordinary integration with respect to \(db\)—in the stated mV coordinate. Equivalently, write the variance constraint as the second-moment constraint

\[ \mathbb{E}[b^2]=b_0^2+\tau_0^2. \]

Lecture 11’s maximum-entropy result gives a density of the form

\[ \pi(b) = \frac{1}{Z} \exp\left(\eta_1b+\eta_2b^2\right), \qquad b\in\mathbb{R}, \]

where \(Z\) is the normalization constant.

  1. Explain why \(\eta_2<0\) is required. (1 point)

  2. Complete the square and use the two moment constraints to determine \(\eta_1\) and \(\eta_2\). (4 points)

  3. Determine \(Z\) and identify the resulting named prior distribution. (2 points)

You may use

\[ \int_{-\infty}^{\infty} \exp\left[-\frac{a}{2}(x-c)^2\right]\,dx = \sqrt{\frac{2\pi}{a}}, \qquad a>0. \]

Response: Replace this text with your derivation.

1.2 Derive the analytical posterior (10 points)#

Suppose the calibration residuals \(Y_1,\ldots,Y_N\) are conditionally independent given \(b\) and satisfy

\[ Y_i\mid b\sim N(b,\sigma^2), \qquad i=1,\ldots,N, \]

where \(\sigma>0\), \(N\geq1\), and the measurement-noise variance \(\sigma^2\) is known. Let \(y_{1:N}=(y_1,\ldots,y_N)\) denote the observed residuals, and define

\[ \bar y=\frac{1}{N}\sum_{i=1}^N y_i. \]

  1. Write the likelihood as a function of \(b\), retaining all factors that depend on \(b\). (2 points)

  2. Starting from Bayes’ rule,

    \[ p(b\mid y_{1:N}) \propto p(y_{1:N}\mid b)\pi(b), \]

    complete the square and show that

    \[ b\mid y_{1:N}\sim N(b_N,v_N). \]

    Derive explicit expressions for \(b_N\) and \(v_N\). (5 points)

  3. Precision means reciprocal variance. Rewrite \(b_N\) as a weighted average of \(b_0\) and \(\bar y\) whose weights are determined by prior precision and data precision. Verify that the weights sum to one and interpret them. (2 points)

  4. A prior family is conjugate to a likelihood when the posterior belongs to the same distribution family as the prior. Explain why this update is conjugate. (1 point)

You may use

\[ \sum_{i=1}^N(y_i-b)^2 = \sum_{i=1}^N(y_i-\bar y)^2 + N(b-\bar y)^2. \]

Leave the posterior in terms of \(N\) and \(\bar y\); no numerical observations are required.

Response: Replace this text with your posterior derivation and interpretation.

1.3 Design the calibration sample size (5 points)#

Let \(v_*>0\) be a target upper bound for the posterior variance. When \(N=0\), define \(v_0=\tau_0^2\), the prior variance.

  1. Derive the minimum number \(N_{\min}\) of calibration measurements required to ensure

    \[ v_N\leq v_*. \]

    Allow \(N_{\min}=0\) when the prior already satisfies the target. You may use \(\lceil x\rceil\), the smallest integer greater than or equal to \(x\). (3 points)

  2. Evaluate \(N_{\min}\) when

    \[ \tau_0=1.00\ \text{mV}, \qquad \sigma=1.00\ \text{mV}, \qquad v_*=0.0625\ \text{mV}^2. \]

    Verify that one fewer measurement fails to meet the target. (1 point)

  3. Explain why \(N_{\min}\) does not depend on the observed residual values. (1 point)

Response: Replace this text with your sample-size calculation and explanation.

1.4 Predict a future calibration residual (3 points)#

Let \(Y_*\) denote a future residual from the same sensor, with

\[ Y_*\mid b\sim N(b,\sigma^2), \]

independently of the calibration data given \(b\). The posterior predictive distribution describes an unobserved future measurement after averaging over posterior uncertainty in \(b\).

  1. Derive the posterior predictive distribution of \(Y_*\) and its central 95% predictive interval. You may use \(P(-1.96\leq Z\leq1.96)\approx0.95\) for \(Z\sim N(0,1)\). For the numerical design at \(N=15\), also report the predictive variance and the 95% predictive half-width. Because no numerical value of \(\bar y\) is supplied, leave this numerical interval centered at the symbolic posterior mean \(b_{15}\). (2 points)

  2. Interpret the two contributions to its variance. State which contribution can be reduced by collecting more calibration measurements and identify the limiting predictive variance as \(N\to\infty\). (1 point)

Problem 1 credit. The maximum-entropy construction follows Jaynes (1957) and Lecture 11. The Normal–Normal update is standard conjugate Bayesian analysis. The sensor calibration setting, sample-size design question, and numerical values were developed for this course.

Response: Replace this text with your posterior-predictive derivation and interpretation.

Problem 2 — Bayesian earthquake-rate study (75 points)#

The course’s frozen USGS ComCat snapshot contains catalog events with magnitude at least 6.5 from January 1, 1900 through December 31, 2025 in the rectangular study region \(32\le\text{latitude}\le37\) degrees N and \(-122\le\text{longitude}\le-114\) degrees. This rectangle includes parts of California and Baja California; it is not a Southern-California boundary.

Let \(\Delta_t=1\) year denote the known exposure for calendar year \(t\). Model the annual counts as conditionally independent observations

\[ Y_t\mid\lambda,\Delta_t \sim\operatorname{Poisson}(\lambda\Delta_t), \]

where \(\lambda\) is a constant event rate in events/year. This is a compact teaching model, not an official seismic-hazard forecast.

Parameterization note. Throughout the notebook, Gamma distributions use shape \(\alpha\) and rate \(\beta\). SciPy instead expects stats.gamma(a=alpha, scale=1 / beta).

For \(\lambda>0\), the shape-rate density and its moments are

\[ f(\lambda\mid\alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1}e^{-\beta\lambda}, \qquad \mathbb{E}[\lambda]=\frac{\alpha}{\beta}, \qquad \mathbb{V}[\lambda]=\frac{\alpha}{\beta^2}. \]

Gamma–Poisson update card. If \(\lambda\sim\operatorname{Gamma}(\alpha_0,\beta_0)\), total observed count is \(n=\sum_t y_t\), and total known exposure is \(E=\sum_t\Delta_t\), then

\[ \lambda\mid\mathbf y,\boldsymbol\Delta \sim\operatorname{Gamma}(\alpha_N,\beta_N), \qquad \alpha_N=\alpha_0+n, \quad \beta_N=\beta_0+E. \]

For a future exposure of \(T\) years, the corresponding probability of at least one event is

\[ P(Y_*\ge1\mid\mathbf y,\boldsymbol\Delta) =1-\left(\frac{\beta_N}{\beta_N+T}\right)^{\alpha_N}. \]

These formulas are supplied; you are responsible for using and interpreting them, not deriving them.

events = (
    pd.read_csv(BytesIO(DATA_BYTES), parse_dates=["time"])
    .sort_values("time", kind="stable")
    .reset_index(drop=True)
)

YEARS = np.arange(1900, 2026)
annual_counts = (
    events["time"].dt.year.value_counts()
    .reindex(YEARS, fill_value=0)
    .sort_index()
    .to_numpy(dtype=int)
)
EXPOSURE_YEARS = YEARS.size
assert annual_counts.shape == (126,)
events.head()

2.1 Audit and represent the observations (10 points)#

Verify the following without changing the frozen CSV:

  • unique event identifiers;

  • that every event lies inside the stated query window, magnitude threshold, and geographic bounds;

  • missingness in every column; identify which missing values, if any, affect the annual-count likelihood or the frozen-extract audit;

  • the total event count, number of zero-event years, and maximum annual count.

Plot the 126 annual counts. Explain why the zero-event years must remain in the exposure and likelihood. Identify one catalog-specific limitation that should be considered before treating the entire 1900–2025 window as exchangeable. Here, exchangeable means that the joint model is unchanged if the annual observations are permuted; in this setting, it rules out a systematic dependence of the event rate on calendar year.

# YOUR CODE HERE

Response: Replace this text with your audit findings and explanation.

2.2 Prior information and sensitivity (15 points)#

Compare these Gamma shape-rate priors for \(\lambda\):

  • Maximum-entropy mean-only: \(\operatorname{Gamma}(1,8)\); positive support, mean \(0.125\) events/year, and no other quantitative constraint relative to \(q(\lambda)=1\) with respect to \(d\lambda\) in the stated events/year coordinate.

  • Same-mean concentrated: \(\operatorname{Gamma}(4,32)\); hypothetical stronger prior information with the same mean.

  • Lower-rate sensitivity: \(\operatorname{Gamma}(2,40)\); hypothetical alternative with mean \(0.05\) events/year.

Report each prior mean, standard deviation, and central 90% interval in events/year. Plot all three densities on a common, readable rate range. Explain why the first prior is exponential and how the first two priors can have the same mean but express different information.

# YOUR CODE HERE

Response: Replace this text with your prior comparison.

2.3 Posterior inference, credible intervals, and Bayes actions (20 points)#

Implement the supplied Gamma–Poisson update using total count and total exposure. For each prior, report the posterior shape and rate, mean, median, mode, and central 95% credible interval for \(\lambda\). Express rate summaries in both events/year and events/decade. Plot the three posterior densities.

State which point summary is the Bayes action under squared-error loss and which is the Bayes action under absolute-error loss. Report the continuous posterior mode only as a density summary; do not use exact 0–1 loss to justify it for a continuous parameter.

def gamma_poisson_update(alpha, beta, counts, exposure_years):
    '''Return posterior shape and rate for a Gamma–Poisson rate model.'''
    # YOUR CODE HERE
    raise NotImplementedError


# Use the function above to build a posterior summary and density plot.
# YOUR CODE HERE

Response: Replace this text with your posterior comparison and loss-based interpretation.

2.4 Posterior-predictive model check (15 points)#

Use the maximum-entropy prior’s posterior for this check. Generate 50,000 replicated 126-year catalogs by repeating these two steps for each catalog:

  1. draw \(\lambda^{\mathrm{rep}}\) from the posterior;

  2. draw 126 conditionally independent annual counts from \(\operatorname{Poisson}(\lambda^{\mathrm{rep}})\).

Each replication is a synthetic 126-year catalog generated under the posterior-predictive model. The sampled rate \(\lambda^{\mathrm{rep}}\) remains fixed across all 126 years within that catalog and is redrawn for the next catalog. Variation in \(\lambda^{\mathrm{rep}}\) across catalogs propagates posterior uncertainty about the unknown rate; the Poisson draws represent year-to-year count variability at a fixed rate.

Define the discrepancy

\[ T(\mathbf y)=\frac{s_y^2}{\bar y}, \]

where \(s_y^2\) is the sample variance across years. The discrepancy compares the across-year variability of the counts with their average annual level. A Poisson random variable has variance equal to its mean, so \(T=1\) is the natural reference. Values above 1 indicate overdispersion, and values below 1 indicate underdispersion. For a finite, sparse catalog, compare the observed value with its posterior-predictive distribution rather than judging it against 1 alone.

In this snapshot every observed annual count is either 0 or 1. Before simulating, show that if \(N>1\) and \(1\leq n\leq N\) events occur across \(N\) such annual counts, then

\[ T(\mathbf y)=\frac{N-n}{N-1}. \]

Explain one kind of temporal behavior that this observed statistic cannot detect.

The statistic is undefined for a replicated catalog whose total count is zero. Exclude any such catalog and report how many were excluded. Report the observed value and, conditional on a positive replicated total, the central 95% posterior-predictive interval for \(T\) and the Bayesian tail area

\[ P\left( T(\mathbf Y^{\mathrm{rep}})\geq T(\mathbf y) \mid \mathbf y,\, \sum_tY_t^{\mathrm{rep}}>0 \right). \]

Plot the replicated discrepancy distribution with the observed value marked. Interpret what the check can and cannot establish about the homogeneous Poisson model.

def simulate_replicated_catalogs(
    alpha_post,
    beta_post,
    n_years,
    n_replications=50_000,
    seed=SEED,
):
    '''Draw replicated annual-count catalogs from the posterior predictive.'''
    # Return an array with shape (n_replications, n_years).
    # YOUR CODE HERE
    raise NotImplementedError


def dispersion_index(count_matrix):
    '''Return sample-variance/mean along the final axis.'''
    # For a two-dimensional input, return one discrepancy per row.
    # YOUR CODE HERE
    raise NotImplementedError


# Carry out the posterior-predictive calculation. Store the retained
# positive-total replicated discrepancies in finite_discrepancies and define
# observed_discrepancy, excluded_discrepancies, predictive_interval, and
# bayesian_tail_area. Use the same retained sample for the interval, tail
# area, and plot.
# YOUR CODE HERE

# Visualization code is supplied; it uses the quantities computed above.
fig, ax = plt.subplots(figsize=(8.5, 4.5))
ax.hist(
    finite_discrepancies,
    bins=60,
    density=True,
    color="#90cdf4",
    edgecolor="white",
    label="Replicated catalogs",
)
ax.axvline(
    observed_discrepancy,
    color="#c53030",
    linewidth=2.5,
    label=fr"Observed $T={observed_discrepancy:.3f}$",
)
ax.set(
    xlabel=r"Annual-count dispersion index $T$",
    ylabel="Posterior-predictive density",
    title="Posterior-predictive check of annual-count dispersion",
)
ax.legend()
sns.despine(ax=ax)
plt.show()

Response: Replace this text with your posterior-predictive interpretation.

2.5 Ten- and twenty-year probabilities, sensitivity, and communication (15 points)#

For each prior, calculate the exact posterior-predictive probabilities \(P(Y_*\ge1\mid\mathbf y)\) for \(T=10\) and \(T=20\) years. Compare each exact value with the plug-in approximation

\[ 1-\exp[-\mathbb{E}(\lambda\mid\mathbf y)T]. \]

A hypothetical preparedness rule triggers an action when the exact 20-year probability is at least \(0.90\). State the decision under each prior.

Then write 150–200 words for a technically literate decision maker. Your explanation must distinguish expected count from probability of at least one event, report the prior-sensitivity ranges for both the 10- and 20-year horizons, state the action-rule result, and identify at least two concrete limitations—one catalog-specific and one model-specific. Do not present this calculation as an official forecast.

# YOUR CODE HERE

Response (150–200 words): Replace this text with your probability communication.

Problem 2 data source and credit#

The automatically loaded CSV is a frozen extract from the U.S. Geological Survey’s ANSS Comprehensive Earthquake Catalog (ComCat) (catalog DOI 10.5066/F7MS3QZH), queried through the official FDSN event web service. The exact query requested CSV output, start time 1900-01-01, end time 2026-01-01, latitude 32–37 degrees N, longitude −122 to −114 degrees, minimum magnitude 6.5, and ascending event time. The snapshot was retrieved August 17, 2026; only time, latitude, longitude, depth, mag, place, and id were retained and sorted by time.

Frozen-file SHA-256: a4e7fb0d2a99fc81cf0f309e13a853b6025413e84d61ab82f18cc8c43ec37072.

The Bayesian Poisson treatment of uncertain earthquake rates and future occurrence probabilities is informed by Love (2012) and Rotondi (2013). The frozen California–Baja snapshot, prior-information exercise, posterior-predictive check, and 10/20-year comparison were assembled for this assignment.

The setup cell first accepts a matching local course copy and otherwise downloads the frozen Homework 3 course-data file. It verifies the bytes against the checksum before parsing them. The USGS query above documents provenance only; this notebook never executes that mutable live query.

USGS catalogs can be revised. This fixed file is used for reproducibility. The reported quantity is a model-based probability of at least one catalog event under the constant-rate Poisson model, not earthquake risk or an authoritative California forecast. For California’s scientific long-term forecast resources, see the California Geological Survey’s UCERF page.

Submission checklist#

  • Both numbered problems total 100 points.

  • Every code cell runs from top to bottom and the setup reports a verified data source; no manual data upload is required.

  • All figures have labels and units.

  • The 150–200 word communication is included.

  • No response cell still contains placeholder text.

Contents