A big earthquake has just struck. What will the next week bring?¶
EPS 88 · PyEarth. Open your own copy on DataHub: click here.
On the morning of 4 July 2019 a magnitude 6.4 earthquake shook the desert near the town of Ridgecrest, California. For a day and a half it looked like the main event. Then, 34 hours later, a magnitude 7.1 struck close by, the largest earthquake in California in twenty years.
After a large earthquake everybody asks the same thing: is it over, or is another one coming? Nobody can predict the next earthquake. What seismologists can give is a forecast, like a weather forecast: how many earthquakes to expect in the coming days, and the probability that one of them is large.
Today you make that forecast yourself, for the week that began ten days after the magnitude 7.1, and then check it against what happened.
Every place you write something opens with a pencil icon and the words Your turn, and is
followed by a cell that already has the shape of the answer written in, with ... wherever your
code goes. Replace every ..., run every cell so your answers and figures are
saved in the file, then download this notebook itself — the .ipynb file — and upload it
to Gradescope. Not a PDF: the marking reads your notebook, and a PDF cannot be read.
In JupyterLab: File ▸ Download, or right-click the file in the left-hand panel and choose
Download.
Two habits from the first minute. A cell runs when you press Shift+Enter, and the notebook remembers everything it has already run — so when something breaks and you cannot see why, Kernel → Restart Kernel and Run All Cells throws the memory away and rebuilds it from the top. That is never the wrong thing to do.
What you'll be able to do¶
The science. Say what a large earthquake does to the earthquakes around it, how fast its aftershocks die away, and how a forecast for the coming week is made and read.
The statistics. A rate · a probability as a proportion · an expected number · a random experiment and the Monte Carlo method · the distribution it gives, its mean and its standard deviation · the three rules for "not", "and" and "at least one" · a forecast compared with what came.
The code. A proportion in one line, (condition).mean() · random draws you can repeat,
np.random.default_rng and rng.choice · a loop that lives through the same week 2,000 times ·
.mean() and .std() of the results.
Eight places where you write something: five in class, three at home. Each is headed
Your turn, and the cell under it already has the shape of the answer written in, with ...
wherever your code goes. Replace every ...; a comment beside each one says what it is for. A
... left in place is not code, and the cell will not do what you expect until it is replaced.
- What did the mainshock change?
- How many aftershocks will the next week bring?
- How big will they be?
- What is the probability of a magnitude 4 or more next week?
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# house style, set once, so every plot cell below holds only what matters
plt.rcParams.update({"figure.dpi": 110,
"axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True})
CACHE = "https://raw.githubusercontent.com/AI4EPS/EPS88_PyEarth/main/data"
def load(name, start, end, region):
"""Fetch one window of the USGS catalogue; fall back to the copy stored with the course.
`region` is extra query text limiting the search to one area.
"""
# Ask the live archive first. If it is down, or you are offline, read the copy stored with
# the course instead, so the notebook still runs.
try:
return pd.read_csv(f"https://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&orderby=time-asc"
f"&starttime={start}&endtime={end}&minmagnitude=2.5"
+ region)
except Exception as e:
print("live source unreachable, using the cached copy:", type(e).__name__)
return pd.read_csv(CACHE + "/" + f"week05_{name}_{start}_{end}_M2.5.csv")
# Every earthquake of magnitude 2.5 or more around two California earthquakes: Ridgecrest in 2019
# for class, Loma Prieta in 1989 for the homework.
quakes = load("ridgecrest", "2019-01-01", "2020-07-01",
"&minlatitude=35.2&maxlatitude=36.4&minlongitude=-118.2&maxlongitude=-117.0")
loma = load("lomaprieta", "1989-01-01", "1990-10-18",
"&latitude=37.04&longitude=-121.88&maxradiuskm=60")
# One new column in each table, `days`: the time since its largest earthquake, the mainshock.
# Negative means before it. This uses dates, which you are not expected to follow yet.
for table in [quakes, loma]:
times = pd.to_datetime(table["time"])
mainshock = times[table["mag"].idxmax()]
table["days"] = (times - mainshock).dt.total_seconds() / 86400
pd.set_option("display.width", 100) # room for a table of earthquakes to print on one line
print("Ridgecrest:", quakes.shape, " Loma Prieta:", loma.shape)
live source unreachable, using the cached copy: URLError
live source unreachable, using the cached copy: URLError Ridgecrest: (2945, 23) Loma Prieta: (644, 23)
1. What did the mainshock change?¶
An earthquake is a sudden slip on a fault, a break in the rock of the Earth's crust. Its size is the magnitude: each step of 1 is about 32 times the energy, as in the plotting week.
The data. quakes is a piece of the public earthquake catalogue of the US Geological Survey
(USGS): every event of magnitude 2.5 or more in a box about 120 km across around Ridgecrest,
from the start of 2019 to the middle of 2020. A magnitude 2.5 is about the smallest that
people nearby can feel. One row is one event (all but two are earthquakes; the two are quarry
blasts). Today you need four columns:
mag: the magnitude;latitudeandlongitude: the place;days: the time since the magnitude 7.1, in days. Negative means before it; 1.5 means a day and a half after it.
The plot shows every row: time across, magnitude up.
plt.figure(figsize=(7, 3.6))
plt.scatter(quakes["days"], quakes["mag"], s=4)
plt.xlim(-190, 60)
plt.xlabel("days since the magnitude 7.1")
plt.ylabel("magnitude")
plt.title(f"Ridgecrest, six months before to two months after (n = {len(quakes)} in all)")
plt.show()
print(quakes[quakes["mag"] >= 6][["days", "mag", "place"]]) # the two largest
days mag place 25 -1.406991 6.4 Ridgecrest Earthquake Sequence 325 0.000000 7.1 Ridgecrest Earthquake Sequence
For six months there is a dot every week or so: the ordinary activity of the region, called the background. Then come the magnitude 6.4, the magnitude 7.1 34 hours after it, and a cloud of earthquakes that thins out as the days pass. Together they are a sequence, and its pieces have names:
- the mainshock is the largest earthquake of the sequence, here the magnitude 7.1;
- an aftershock is one of the earthquakes that follow it. A large earthquake changes the stress in the rock around it, and parts of the same fault and faults near it that were already close to failing then fail too — most within hours and days, fewer and fewer over the months that follow.
- a foreshock is an earthquake followed by a larger one nearby. The magnitude 6.4 was one, and nobody could tell until the larger one came.
✏️ Your turn 1¶
How many earthquakes a day, before the sequence and after the mainshock?
Why: the two numbers together say how much the mainshock changed the region.
A rate is a count divided by the time it took.
$$\text{rate} = \frac{\text{count}}{\text{time}}$$
30 earthquakes in 10 days is a rate of 30 / 10 = 3 a day. A rate lets you compare periods of different length.
Work out the rate for two groups of rows. Pick each group with a mask on days; len counts
its rows.
before: the quiet period, more than 2 days before the mainshock, which isdays < -2. It is 184 days long. (The cut at 2 days keeps the foreshock out of it.)first_month: the first 30 days after the mainshock, which isdays > 0anddays < 30. The mainshock itself hasdaysequal to 0, so> 0leaves it out.
Use these names, because the self-check looks for them: before, first_month, rate_before,
rate_first_month.
before = quakes[quakes["days"] < -2]
first_month = quakes[(quakes["days"] > 0) & (quakes["days"] < 30)]
rate_before = len(before) / 184
rate_first_month = len(first_month) / 30
print("before:", len(before), "earthquakes in 184 days,", round(rate_before, 2), "a day")
print("first month:", len(first_month), "earthquakes in 30 days,", round(rate_first_month), "a day")
print("the rate went up", round(rate_first_month / rate_before), "times")
before: 24 earthquakes in 184 days, 0.13 a day first month: 2186 earthquakes in 30 days, 73 a day the rate went up 559 times
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (before["days"] < -2).all() and len(before) == (quakes["days"] < -2).sum(), "before keeps every row whose days is below -2, and no other"
assert (first_month["days"] > 0).all() and (first_month["days"] < 30).all() and len(first_month) == ((quakes["days"] > 0) & (quakes["days"] < 30)).sum(), "first_month keeps every row with days above 0 and below 30, and no other: two conditions joined with &"
assert abs(rate_before - len(before) / 184) < 1e-9 and abs(rate_first_month - len(first_month) / 30) < 1e-9, "each rate is the number of rows, len(...), divided by the number of days"
print(f"✓ Your turn 1 — {round(rate_before, 2)} earthquakes a day before, {round(rate_first_month)} a day in the first month")
✓ Your turn 1 — 0.13 earthquakes a day before, 73 a day in the first month
One earthquake a week became 73 a day: the rate went up more than five
hundred times. The map shows where they were: your before in grey, everything after the
mainshock in red.
after = quakes[quakes["days"] > 0] # every row after the mainshock; before holds your rows of Your turn 1
plt.figure(figsize=(6, 5.4))
plt.scatter(after["longitude"], after["latitude"], s=2, color="crimson", label="after the mainshock")
plt.scatter(before["longitude"], before["latitude"], s=14, color="0.35", label="in the six months before")
plt.scatter(-117.6709, 35.6225, s=60, marker="s", color="black")
plt.text(-117.6709 - 0.33, 35.6225 - 0.02, "Ridgecrest")
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title(f"Before (n = {len(before)}) and after (n = {len(after)})")
plt.legend(loc="lower left")
plt.show()
Aftershocks happen on and around the fault that slipped in the mainshock. Most of the red dots draw one band about 50 km long, running from north-west to south-east past the town: that is the fault.
The rate jumped, and since then it has been falling. To forecast the days ahead you need to know how fast it falls.
2. How many aftershocks will the next week bring?¶
Put yourself ten days after the mainshock. You make the forecast with what is known on that day. Then you compare it with what came. From here on, every earthquake in the box counts as an aftershock: at 73 a day against one a week, nearly all of them are.
First count the aftershocks on each of the ten days. Day 1 is the first 24 hours after the
mainshock, which is days from 0 to 1. Day 2 is days from 1 to 2, and so on. The week to
forecast, days 11 to 17, is days from 10 to 17.
In a mask, a period starts with >= and ends with <. For example,
(quakes["days"] >= 1) & (quakes["days"] < 10) reads "from 1 day after the mainshock, up to
but not including 10 days after it". With this rule every aftershock falls in exactly one
period. One exception: a period that starts at the mainshock starts with > 0, to leave out the
mainshock itself, whose days is 0.
plt.hist draws the count of each day as a bar. np.histogram with the same bins gives the
counts as numbers, past_counts, and past_days holds the day numbers that go with them.
first_ten = quakes[(quakes["days"] > 0) & (quakes["days"] < 10)] # the first 10 days after the mainshock, without the mainshock itself
plt.figure(figsize=(7, 3.6))
plt.hist(first_ten["days"], bins=np.arange(0, 11)) # one bin for each day
plt.xlabel("days since the mainshock")
plt.ylabel("aftershocks in the day")
plt.title(f"The first ten days (n = {len(first_ten)} aftershocks)")
plt.show()
past_counts, edges = np.histogram(first_ten["days"], bins=np.arange(0, 11)) # the same counts, as numbers
past_days = np.arange(1, 11) # the day numbers that go with them: 1, 2, ... 10
print("day number: ", past_days)
print("aftershocks that day:", past_counts)
day number: [ 1 2 3 4 5 6 7 8 9 10] aftershocks that day: [1257 170 104 76 86 70 37 29 31 34]
The count falls fast at first, then more and more slowly: 1257 on day 1, 170 on day 2, 76 on day 4, 29 on day 8. From day 2 on, each time the day number doubles, the count roughly halves.
Omori's law. The number of aftershocks in a day falls roughly as one over the time since the mainshock: twice as late, about half as many. Fusakichi Omori found it in Japan in 1894. As an equation:
$$\text{count on day } t = \frac{K}{t}$$
- $t$ is the day number after the mainshock: 1, 2, 3, ...
- $K$ is one number that belongs to the sequence. A larger $K$ means more aftershocks on every day.
Read it with numbers. On day 2 the equation gives $K / 2$ and on day 4 it gives $K / 4$:
twice as late, half as many. In Python, K / past_days works the equation out for all ten days
at once.
An equation that stands in for the data is called a model. Once you know $K$, the model is a forecast: put in a day that has not happened yet.
✏️ Your turn 2¶
Which K fits this sequence, and how many aftershocks should the next seven days bring?
Why: the model can forecast only when its K is known, and K is found from the days already counted.
Step 1: find K. The figure shows the count of each past day as a dot, and the model,
K / past_days, as a red curve. Start with K = 100, run the cell, and change K until the curve
passes through the dots from day 2 on. No value fits every dot; stop when about as many dots
are above the curve as below it. The dot of day 1 stays far above the curve: in the first hours
the count falls too fast for so simple a formula.
Step 2: forecast. Use the same model on the coming days, days 11 to
17. expected_daily is the count the model expects on each of them, $K / t$ for
$t$ = 11, ..., 17, and expected_week is their sum:
$$\text{expected in the week} = \frac{K}{11} + \frac{K}{12} + \dots + \frac{K}{17}$$
Use these names, because the self-check looks for them: K, expected_daily,
expected_week.
# Step 1: find K. From the cell above: past_days holds the day numbers 1 to 10, past_counts the count of each
K = 320
plt.figure(figsize=(7, 3.6))
plt.scatter(past_days, past_counts, label="the count of each day")
plt.plot(past_days, K / past_days, color="firebrick", label="the model, K / day")
plt.yscale("log") # a log axis, so that 1,000 and 30 can both be read
plt.xlabel("day after the mainshock")
plt.ylabel("aftershocks in the day")
plt.title(f"Omori's law with K = {K} (n = {len(past_days)} days)")
plt.legend()
plt.show()
# Step 2: forecast
coming_days = np.arange(11, 18) # the day numbers 11 to 17: the week to forecast
expected_daily = K / coming_days
expected_week = expected_daily.sum()
print("K =", K)
print("expected on each of days 11 to 17:", expected_daily.round())
print("expected in days 11 to 17:", round(expected_week, 1))
K = 320 expected on each of days 11 to 17: [29. 27. 25. 23. 21. 20. 19.] expected in days 11 to 17: 163.4
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert 250 <= K <= 400, "choose K so that the curve passes through the dots from day 2 on: a value between about 250 and 400"
assert len(expected_daily) == 7 and abs(expected_daily[0] - K / coming_days[0]) < 1e-9, "expected_daily is the model on the COMING days, K / coming_days: seven expected counts, for days 11 to 17"
assert abs(expected_week - expected_daily.sum()) < 1e-9, "expected_week is the seven daily counts added up"
print(f"✓ Your turn 2 — with K = {K}, about {round(expected_week)} aftershocks are expected in days 11 to 17")
✓ Your turn 2 — with K = 320, about 163 aftershocks are expected in days 11 to 17
Now compare your forecast with what came. The catalogue holds the real aftershocks of days 11 to 17. The cell counts them day by day and draws them, in orange, beside your forecast.
came = quakes[(quakes["days"] >= 10) & (quakes["days"] < 17)] # the aftershocks that came in days 11 to 17: from 10 to 17 days after the mainshock
came_counts, edges = np.histogram(came["days"], bins=np.arange(10, 18)) # the count of each of days 11 to 17
plt.figure(figsize=(7, 3.6))
plt.scatter(past_days, past_counts, label="the count of each day, days 1 to 10")
plt.plot(past_days, K / past_days, color="firebrick", label="the model, K / day")
plt.plot(coming_days, expected_daily, color="firebrick", ls="--", label="your forecast, days 11 to 17")
plt.scatter(coming_days, came_counts, color="darkorange", label="what came")
plt.yscale("log")
plt.xlabel("day after the mainshock")
plt.ylabel("aftershocks in the day")
plt.title(f"Your forecast and what came (n = {len(came)} aftershocks came)")
plt.legend()
plt.show()
print("your forecast for each of days 11 to 17:", expected_daily.round())
print("what came on each of them: ", came_counts)
print("your forecast for the week:", round(expected_week), " what came:", len(came))
your forecast for each of days 11 to 17: [29. 27. 25. 23. 21. 20. 19.] what came on each of them: [41 31 28 28 19 13 9] your forecast for the week: 163 what came: 169
With K = 320 the forecast for the week is 163, and 169 came. Day by day it is rougher: the first days brought more than the model said (41 against 29 on day 11) and the last days fewer (9 against 19 on day 17).
Think. Is this a good forecast? Which would you give to the people of Ridgecrest on day 10: the count for the week, or the count for each day?
So that everyone's numbers agree, the notebook uses the forecast of 163 from here on, under
the name n_week.
3. How big will they be?¶
About 163 aftershocks are expected next week. Most will be small. What is the probability that an aftershock is large?
What a probability is¶
Probability. The proportion of many tries in which a thing happens.
$$P(A) = \frac{\text{number of tries in which } A \text{ happens}}{\text{number of tries}}$$
$A$ stands for the thing you ask about, such as "the die shows a six", and $P(A)$ is read "the probability of A". It is a number between 0 and 1: 0 means never, 1 means always. It is also written as a percentage: 0.3 is 30%.
- A coin lands heads with probability 0.5: in many tosses, about half are heads.
- A die shows a six with probability 1/6, about 0.17: in many rolls, about one in six is a six.
- A forecast of "30% rain" says: of many days like this one, it rains on about 3 in 10.
A probability does not say what the next try will be. The next roll is a six or it is not. It says what happens over many tries.
Where does the number come from? For a die you can reason it out: six faces, all alike, so 1/6. For an aftershock nobody can. So you measure it: take many aftershocks that have already happened, and find the proportion of them that were large.
Measuring a probability from the data¶
The tries are the aftershocks already recorded. The values you have are called a
sample. They stand for all the aftershocks this sequence can produce, most of which have not
happened yet. sample holds the aftershocks from 1 day to
10 days after the mainshock: days >= 1 and days < 10.
- nothing after 10 days, because the forecast is made on day 10 and can use only what has happened by then;
- not the first 24 hours, because the catalogue of that day is not complete. Right after a large earthquake the waves of many aftershocks overlap, and only the larger ones can be picked out. Counting that day would make large aftershocks look more common than they are.
The cell draws the magnitudes of the sample.
sample = quakes[(quakes["days"] >= 1) & (quakes["days"] < 10)] # from 1 day to 10 days after the mainshock: the first 24 hours are left out
plt.figure(figsize=(9, 3.6))
plt.subplot(1, 2, 1)
plt.hist(sample["mag"], bins=np.arange(2.5, 5.2, 0.1))
plt.axvline(4, color="firebrick")
plt.xlabel("magnitude")
plt.ylabel("number of aftershocks")
plt.title("ordinary axis")
plt.subplot(1, 2, 2)
plt.hist(sample["mag"], bins=np.arange(2.5, 5.2, 0.1))
plt.axvline(4, color="firebrick")
plt.yscale("log")
plt.xlabel("magnitude")
plt.title(f"log axis (n = {len(sample)} aftershocks)")
plt.show()
Small aftershocks far outnumber large ones. The bars fall away quickly to the right. The red line marks magnitude 4, large enough to be felt well beyond the place where it happens. The few bars to its right show only on the log axis.
How to measure a proportion. The cell does it step by step, on the first ten magnitudes of the sample: what proportion of them is 3 or more?
ten = np.array(sample["mag"])[:10] # the first ten magnitudes of the sample
print("ten magnitudes: ", ten)
print("3 or more? ", ten >= 3)
print("how many are True: ", (ten >= 3).sum(), "of 10")
print("the proportion: ", (ten >= 3).mean())
ten magnitudes: [2.83 3.3 2.52 2.65 3.12 3.32 2.9 3.05 2.84 4.53] 3 or more? [False True False False True True False True False True] how many are True: 5 of 10 the proportion: 0.5
A proportion is the mean of a mask. ten >= 3 asks "is this one magnitude 3 or more?" of
every value and answers True or False. 5 of the ten are True: a proportion of
0.5. Python counts True as 1 and False as 0, so .mean() of the mask gives the
proportion in one step.
The proportion in the whole sample is the probability to use for next week. Ten values are too few to trust; the sample has 637. With $M$ for the magnitude of an aftershock:
$$P(M \geq 4) \approx \frac{\text{aftershocks in the sample with } M \geq 4}{\text{aftershocks in the sample}}$$
The sign ≈ means "is about": a proportion measured on a sample is an estimate of the probability, not the probability itself. It also assumes one thing: the sizes of next week's aftershocks will be like the sizes in the sample. You will check that against what came.
Expected number. How many to expect on average: the number of tries times the probability of each.
$$E = n \times p$$
- $n$ is the number of tries.
- $p$ is the probability of each.
- $E$ is the expected number.
With 10 tries and a probability of 0.5, $E$ = 10 × 0.5 = 5. It is an average, so it does not have to be a whole number.
✏️ Your turn 3¶
What is the probability that an aftershock is magnitude 3 or more, and magnitude 4 or more? How many of each should next week bring?
Why: the proportion in the sample is the probability for one aftershock, and with it you can say how many large aftershocks to expect among next week's.
prob_m3andprob_m4: the proportion of the wholesamplethat is magnitude 3 or more, and magnitude 4 or more, as the cell above did on ten values.expected_m3andexpected_m4: the expected number of each. The tries are next week's aftershocks,n_weekof them.
Use these names, because the self-check looks for them: prob_m3, prob_m4, expected_m3,
expected_m4.
n_week = 163 # aftershocks expected next week, days 11 to 17 (Your turn 2)
prob_m3 = (sample["mag"] >= 3).mean()
prob_m4 = (sample["mag"] >= 4).mean()
expected_m3 = n_week * prob_m3
expected_m4 = n_week * prob_m4
print((sample["mag"] >= 3).sum(), "of", len(sample), "are magnitude 3 or more: a probability of", round(prob_m3, 3))
print((sample["mag"] >= 4).sum(), "of", len(sample), "are magnitude 4 or more: a probability of", round(prob_m4, 3))
print("expected next week:", round(expected_m3), "of magnitude 3 or more,", round(expected_m4, 1), "of magnitude 4 or more")
212 of 637 are magnitude 3 or more: a probability of 0.333 13 of 637 are magnitude 4 or more: a probability of 0.02 expected next week: 54 of magnitude 3 or more, 3.3 of magnitude 4 or more
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert abs(prob_m3 - (sample["mag"] >= 3).mean()) < 1e-9 and 0.30 < prob_m3 < 0.37, "prob_m3 is the proportion of the sample with magnitude 3 or more: the mean of a mask on sample, about a third"
assert abs(prob_m4 - (sample["mag"] >= 4).mean()) < 1e-9 and 0.01 < prob_m4 < 0.03, "prob_m4 is the proportion of the sample with magnitude 4 or more, about 0.02"
assert abs(expected_m3 - n_week * prob_m3) < 1e-9 and abs(expected_m4 - n_week * prob_m4) < 1e-9, "an expected number is the number of tries, n_week, times the probability"
print(f"✓ Your turn 3 — a probability of {round(prob_m3, 3)} for magnitude 3 or more and {round(prob_m4, 3)} for 4 or more: expect {round(expected_m3)} and {round(expected_m4, 1)} next week")
✓ Your turn 3 — a probability of 0.333 for magnitude 3 or more and 0.02 for 4 or more: expect 54 and 3.3 next week
A third of the aftershocks are magnitude 3 or more; about one in fifty is magnitude 4 or more (13 of 637). So of next week's 163, expect about 54 of magnitude 3 or more and about 3.3 of magnitude 4 or more.
Now compare with what came. The cell counts the large aftershocks among the real ones of days 11 to 17.
came = quakes[(quakes["days"] >= 10) & (quakes["days"] < 17)] # the aftershocks that came in days 11 to 17: from 10 to 17 days after the mainshock
came_m3 = (came["mag"] >= 3).sum() # how many of them are magnitude 3 or more
came_m4 = (came["mag"] >= 4).sum() # how many are magnitude 4 or more
print("magnitude 3 or more: expected", round(expected_m3), " came", came_m3)
print("magnitude 4 or more: expected", round(expected_m4, 1), " came", came_m4)
print("proportion of magnitude 3 or more: in the sample", round(prob_m3, 3), " in what came", round(came_m3 / len(came), 3))
print("proportion of magnitude 4 or more: in the sample", round(prob_m4, 3), " in what came", round(came_m4 / len(came), 3))
magnitude 3 or more: expected 54 came 62 magnitude 4 or more: expected 3.3 came 3 proportion of magnitude 3 or more: in the sample 0.333 in what came 0.367 proportion of magnitude 4 or more: in the sample 0.02 in what came 0.018
Expected 54, came 62. Expected 3.3, came 3. The proportions held too: of what came, 0.367 was magnitude 3 or more (0.333 in the sample) and 0.018 was magnitude 4 or more (0.02 in the sample). So next week's aftershocks were like the sample in size, as the forecast assumed.
Think. You expected 3.3 aftershocks of magnitude 4 or more, and 3 came. Could none have come? Could 7? The expected number alone cannot say. The next section can.
4. What is the probability of a magnitude 4 or more next week?¶
Expect about 3.3 aftershocks of magnitude 4 or more. But that is an average. A real week holds a whole number of them: 0, 1, 2, 3 or more. What is the probability that it holds at least one?
This section first teaches the tool on simple cases, a die and the number π, and then uses it on the earthquakes.
Random experiments¶
Random experiment. Something you can do again and again whose result chance decides each time: toss a coin, roll a die. One run tells you little. Many runs tell you the probability: the proportion of runs in which a thing happens.
Why run experiments? Some probabilities can be reasoned out: one die shows a six with probability 1/6. Many cannot, or only with hard mathematics: how many of 163 aftershocks will be large? When you cannot work a probability out, run the experiment many times and count. A computer can run a simple experiment thousands of times in a second.
Python runs a random experiment with two tools:
rng = np.random.default_rng()makes a random-number generator.rng.choice(values, size)drawssizeof the values at random. Every value has the same chance at every draw.
The experiment in the cell: roll a die ten times. Run the cell two or three times.
die = [1, 2, 3, 4, 5, 6] # the six faces
rng = np.random.default_rng() # a random-number generator
rolls = rng.choice(die, 10) # roll the die 10 times
print("ten rolls:", rolls)
ten rolls: [6 5 5 1 6 3 2 4 3 1]
Every run gives different rolls, as a real die does. That is what random means. But it also means your numbers differ from your neighbour's, and nobody can check them.
A seed fixes the draws. A computer's random numbers come from a recipe, and the seed is where the recipe starts, so the same seed gives the same draws every time. The seed is the number in the brackets:
np.random.default_rng(88). Run the next cell two or three times.
die = [1, 2, 3, 4, 5, 6] # the six faces
rng = np.random.default_rng(88) # a generator started from the seed 88
rolls = rng.choice(die, 10)
print("ten rolls, seed 88: ", rolls)
rng = np.random.default_rng(88) # start again from the same seed
rolls_again = rng.choice(die, 10)
print("seed 88 once more: ", rolls_again)
more_rolls = rng.choice(die, 10) # go on drawing from the same generator
print("the next ten rolls: ", more_rolls)
ten rolls, seed 88: [4 2 6 5 5 5 6 5 3 4] seed 88 once more: [4 2 6 5 5 5 6 5 3 4] the next ten rolls: [4 1 5 2 5 2 3 2 6 1]
The same seed gives the same rolls, on every run and on every computer. Going on drawing from the generator still gives new values, spread like those of a real die. This notebook fixes the seed wherever it draws, so that your results can be checked.
With replacement, or without? There are two ways to draw.
- With replacement: a value that has been drawn can be drawn again. A die works this way:
a six does not use the six up.
rng.choice(values, size)draws with replacement. - Without replacement: a value that has been drawn cannot come again. Dealing cards works
this way.
rng.choice(values, size, replace=False)draws without replacement.
rng = np.random.default_rng(88) # a generator started from the seed 88
cards = [1, 2, 3, 4, 5, 6]
with_replacement = rng.choice(cards, 6) # every draw can be any of the six values
print("with replacement: ", with_replacement)
without_replacement = rng.choice(cards, 6, replace=False) # a value that is drawn cannot come again
print("without replacement:", without_replacement)
with replacement: [4 2 6 5 5 5] without replacement: [2 5 3 6 4 1]
With replacement some values repeat and others are missing. Without replacement each value comes exactly once, in a random order.
The Monte Carlo method, on a simple case¶
Monte Carlo. Let the computer live through the same situation thousands of times, each time with fresh random draws, and count how often a thing happens. It has three steps:
- Describe one run of the experiment with random draws.
- Repeat it many times.
- Count: the proportion of the runs in which the thing happens.
A simple case: find the number π with random points. The square in the figure below has sides of length 1, so its area is 1. A circle of radius $r$ has area $\pi r^2$, so the quarter circle of radius 1 inside the square has area $\pi / 4$. A point dropped at random in the square falls inside the quarter circle with probability $\pi / 4$, the share of the square that the quarter circle covers. Turned round:
$$\pi \approx 4 \times \frac{\text{points inside the quarter circle}}{\text{all the points}}$$
- One run: drop one point at random.
rng.random(size)draws numbers at random between 0 and 1, one forxand one fory. The point is inside when its distance from the corner, $\sqrt{x^2 + y^2}$, is less than 1, which is the same as $x^2 + y^2 < 1$. - Repeat: 10,000 points.
- Count: the proportion of the points inside. Four times that proportion should be close to π = 3.14159...
rng = np.random.default_rng(88) # a generator started from the seed 88
n_points = 10000
x = rng.random(n_points) # 10,000 random numbers between 0 and 1
y = rng.random(n_points) # and 10,000 more
inside = x ** 2 + y ** 2 < 1 # True where the point is less than 1 from the corner: inside the quarter circle
outside = x ** 2 + y ** 2 >= 1 # True for the other points
plt.figure(figsize=(4.6, 4.6))
plt.scatter(x[inside], y[inside], s=1, color="firebrick")
plt.scatter(x[outside], y[outside], s=1, color="gray")
plt.gca().set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title(f"Random points: inside in red (n = {n_points})")
plt.show()
proportion = inside.mean() # the proportion of the points that are inside
print("inside:", inside.sum(), "of", n_points, " a proportion of", proportion)
print("4 times the proportion:", 4 * proportion)
inside: 7856 of 10000 a proportion of 0.7856 4 times the proportion: 3.1424
Random points found π. No formula for π was used: only random draws and a count.
How many runs are enough? The cell repeats the experiment with 100 points, with 10,000 and with a million.
rng = np.random.default_rng(88) # a generator started from the seed 88
for n_points in [100, 10000, 1000000]:
x = rng.random(n_points) # n_points random numbers between 0 and 1
y = rng.random(n_points)
proportion = (x ** 2 + y ** 2 < 1).mean() # the proportion of the points inside the quarter circle
print("with", n_points, "points: pi is about", round(4 * proportion, 4))
with 100 points: pi is about 3.24 with 10000 points: pi is about 3.1436 with 1000000 points: pi is about 3.141
More runs give a better answer. With 100 points the estimate is rough; with a million it is right to two decimals or more. A proportion from many runs settles on the probability: this is called the law of large numbers.
Two numbers that describe many results¶
Many runs give many numbers. Two numbers describe them.
The mean is the average: add the values up and divide by how many there are.
$$\bar{x} = \frac{x_1 + x_2 + \dots + x_N}{N}$$
The standard deviation. How far a value usually is from the average.
$$s = \sqrt{\frac{(x_1 - \bar{x})^2 + (x_2 - \bar{x})^2 + \dots + (x_N - \bar{x})^2}{N}}$$
- $x_1, x_2, \dots, x_N$ are the values, and $N$ is how many there are.
- $\bar{x}$, read "x bar", is their mean.
- $s$ is the standard deviation. Read the equation from the inside out: take the distance of each value from the mean, square it, average the squares, and take the square root.
A small case: the values 2, 4 and 6. The mean is (2 + 4 + 6) / 3 = 4. The distances from the
mean are −2, 0 and 2, and their squares are 4, 0 and 4. The average of the squares is 8 / 3 =
2.67, and its square root is 1.63. In Python these are .mean() and .std().
values = np.array([2, 4, 6])
print("mean:", values.mean())
print("standard deviation:", round(values.std(), 2))
mean: 4.0 standard deviation: 1.63
One possible week¶
Now the same three steps for the earthquakes.
- One run is one possible week: 163 magnitudes drawn at random from
sample, the aftershocks already recorded. The thing to look for: at least one of magnitude 4 or more. The draw is with replacement. Next week's aftershocks are new earthquakes, not a selection of the old ones: the sample only says how often each size comes up, and a size that came once can come again. - Repeat: 2,000 possible weeks.
- Count: the proportion of the weeks with at least one.
The cell does step 1 twice.
n_week = 163 # aftershocks expected next week, days 11 to 17 (Your turn 2)
rng = np.random.default_rng(1) # a generator started from the seed 1
first_week = rng.choice(sample["mag"], n_week) # one possible week: n_week magnitudes drawn from the sample
print("its first five magnitudes:", first_week[:5])
print("one possible week:", (first_week >= 4).sum(), "of magnitude 4 or more")
second_week = rng.choice(sample["mag"], n_week) # draw again: another possible week
print("another possible week:", (second_week >= 4).sum(), "of magnitude 4 or more")
its first five magnitudes: [3.11 2.69 3.1 3.33 2.58] one possible week: 2 of magnitude 4 or more another possible week: 1 of magnitude 4 or more
Each run gives a different week, with its own number of aftershocks of magnitude 4 or more. One week says little. Steps 2 and 3 are yours.
✏️ Your turn 4¶
In what proportion of 2,000 possible weeks is there at least one aftershock of magnitude 4 or more?
Why: that proportion is the probability the forecast needs. These are steps 2 and 3 of the Monte Carlo method: repeat, then count.
- Draw 2,000 possible weeks with a loop. Each time, draw one week as the cell above did,
count its magnitudes of 4 or more, and keep that count in the list
counts_m4. After the loop,np.arrayturns the list into an array. prob_week_sim: the proportion of the 2,000 counts that are 1 or more. This is the probability.mean_countandstd_count: the mean and the standard deviation of the 2,000 counts. A small standard deviation means the weeks are much alike; a large one means they differ a lot.
Use these names, because the self-check looks for them: counts_m4, prob_week_sim,
mean_count, std_count.
n_week = 163 # aftershocks expected next week (Your turn 2); sample holds the aftershocks already recorded
rng = np.random.default_rng(88) # a generator started from the seed 88
counts_m4 = [] # one count for each possible week
for run in range(2000):
one_week = rng.choice(sample["mag"], n_week)
counts_m4.append((one_week >= 4).sum())
counts_m4 = np.array(counts_m4)
prob_week_sim = (counts_m4 >= 1).mean()
mean_count = counts_m4.mean()
std_count = counts_m4.std()
print("the counts of the first ten weeks:", counts_m4[:10])
print("at least one in", round(prob_week_sim, 3), "of the weeks")
print("mean", round(mean_count, 1), " standard deviation", round(std_count, 1))
the counts of the first ten weeks: [2 3 5 4 8 3 3 1 4 2] at least one in 0.962 of the weeks mean 3.3 standard deviation 1.8
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(counts_m4) == 2000 and counts_m4.max() < 40, "counts_m4 holds one number for each of the 2000 weeks: how many of that week's aftershocks are magnitude 4 or more"
assert 3.08 < counts_m4.mean() < 3.57 and 1.53 < counts_m4.std() < 2.08, "each week is n_week magnitudes drawn from the sample's mag column, and counts_m4 keeps the COUNT of those that are 4 or more: about 3.3 on average"
assert abs(prob_week_sim - (counts_m4 >= 1).mean()) < 1e-9, "prob_week_sim is the proportion of weeks where counts_m4 is 1 or more: the mean of a mask"
assert abs(mean_count - counts_m4.mean()) < 1e-9 and abs(std_count - counts_m4.std()) < 1e-9, "mean_count is counts_m4.mean() and std_count is counts_m4.std()"
print(f"✓ Your turn 4 — at least one in {round(100 * prob_week_sim)}% of 2000 weeks; mean {round(mean_count, 1)}, standard deviation {round(std_count, 1)}")
✓ Your turn 4 — at least one in 96% of 2000 weeks; mean 3.3, standard deviation 1.8
Now look at your 2,000 weeks, and at what came. The cell draws how often each count came up. The red lines mark the mean and one standard deviation on each side of it. The orange line marks the real week: 3 aftershocks of magnitude 4 or more came.
came_m4 = 3 # what came: the aftershocks of magnitude 4 or more in days 11 to 17
plt.figure(figsize=(7, 3.6))
plt.hist(counts_m4, bins=np.arange(0, counts_m4.max() + 2), align="left") # one bar for each whole number
plt.axvline(mean_count, color="firebrick", label="the mean")
plt.axvline(mean_count - std_count, color="firebrick", ls="--", label="one standard deviation each side")
plt.axvline(mean_count + std_count, color="firebrick", ls="--")
plt.axvline(came_m4, color="darkorange", lw=3, label="what came")
plt.xlabel("aftershocks of magnitude 4 or more in the week")
plt.ylabel("number of weeks")
plt.title("2,000 possible weeks (n = 2,000)")
plt.legend()
plt.show()
within_one_std = (counts_m4 > mean_count - std_count) & (counts_m4 < mean_count + std_count)
print("weeks within one standard deviation of the mean:", round(within_one_std.mean(), 2))
for tries in [10, 100, 1000, 2000]:
proportion = (counts_m4[:tries] >= 1).mean() # the proportion among the first weeks only
print("after", tries, "weeks, the proportion with at least one:", round(proportion, 2))
weeks within one standard deviation of the mean: 0.72 after 10 weeks, the proportion with at least one: 1.0 after 100 weeks, the proportion with at least one: 0.96 after 1000 weeks, the proportion with at least one: 0.97 after 2000 weeks, the proportion with at least one: 0.96
What the 2,000 weeks say¶
The histogram is a distribution: it shows how often each outcome came up. Your three numbers describe it.
- The probability of at least one is 0.962. That is every bar but the first, the weeks with none.
- The mean is 3.3, the solid line. It is the expected number of Your turn 3, found a second way.
- The standard deviation is 1.8. 72% of the weeks fall between the dashed lines, at 1.5 and 5.1. So "about 3.3" means "usually between 2 and 5".
What came: 3. The orange line stands among the tallest bars. The real week was an ordinary one of your 2,000 possible weeks.
Think. Had 8 come, would your forecast have been wrong? 8 or more came up in 36 of the 2,000 possible weeks: rare, but possible.
Why so many weeks? The last lines of the printout answer, as the points did for π. After 10 weeks the proportion reads 1; by 100 it is near 0.96, and it stays there. A probability is the value the proportion settles on.
The rules of probability behind it¶
The same probability can be worked out with three rules and no simulation.
Rule 1: not. The probability that a thing does not happen is 1 minus the probability that it does.
$$P(\text{not } A) = 1 - P(A)$$
One aftershock is magnitude 4 or more with probability 0.02, so it is smaller with probability 1 − 0.02 = 0.98.
Rule 2: and. If two things are independent, the probability that both happen is the product of their two probabilities.
$$P(A \text{ and } B) = P(A) \times P(B)$$
A coin lands tails with probability 0.5, so two coins both land tails with probability 0.5 × 0.5 = 0.25: one time in four.
Independent. Knowing how one turned out tells you nothing about the next. Two coins are independent. Two aftershocks nearly are: the size of one says almost nothing about the size of the next.
Two aftershocks are both smaller than magnitude 4 with probability
0.98 × 0.98 = 0.96. Three are, with probability
0.98 × 0.98 × 0.98 = 0.94. For any number of
aftershocks, multiply 0.98 by itself that many times. Python writes this with **:
0.98 ** 3 is 0.98 × 0.98 × 0.98.
Rule 3: at least one. "At least one is magnitude 4 or more" is the opposite of "all of them are smaller". So by rule 1 its probability is 1 minus the probability from rule 2.
$$P(\text{at least one of } n) = 1 - (1 - p)^n$$
- $p$ is the probability that one aftershock is magnitude 4 or more.
- $1 - p$ is the probability that one is smaller (rule 1).
- $(1 - p)^n$ is the probability that all $n$ are smaller (rule 2, used $n$ times).
- In Python:
1 - (1 - p) ** n.
The cell uses the three rules for every number of aftershocks from 1 to 300, and draws both probabilities.
prob_m4 = 13 / 637 # the probability that one aftershock is magnitude 4 or more (Your turn 3)
numbers = np.arange(1, 301) # 1, 2, 3, ... 300 aftershocks
curve_all_below = (1 - prob_m4) ** numbers # rules 1 and 2: every one of them is below magnitude 4
curve_at_least_one = 1 - curve_all_below # rule 3: at least one is magnitude 4 or more
plt.figure(figsize=(7, 3.6))
plt.plot(numbers, curve_all_below, color="gray", label="all are below magnitude 4")
plt.plot(numbers, curve_at_least_one, color="firebrick", label="at least one is magnitude 4 or more")
plt.xlabel("number of aftershocks")
plt.ylabel("probability")
plt.title("The more aftershocks, the more likely a large one")
plt.legend()
plt.show()
for number in [1, 10, 50, 100]:
prob = 1 - (1 - prob_m4) ** number # the three rules in one line
print("among", number, "aftershocks: at least one with probability", round(prob, 2))
among 1 aftershocks: at least one with probability 0.02 among 10 aftershocks: at least one with probability 0.19 among 50 aftershocks: at least one with probability 0.64 among 100 aftershocks: at least one with probability 0.87
Each aftershock is unlikely to be large, but there are many of them. One is magnitude 4 or more with probability 0.02. Among 10 there is at least one with probability 0.19, among 50 with 0.64, among 100 with 0.87. The two curves always add up to 1, by rule 1: either all are smaller, or at least one is not.
✏️ Your turn 5¶
What is the probability of at least one aftershock of magnitude 4 or more in the whole week, and tomorrow alone?
Why: the rules give in one line what the simulation needed 2,000 weeks for, and they work
for any number of aftershocks. Use them first for the week's n_week aftershocks, to compare
with your prob_week_sim, then for tomorrow's n_tomorrow.
Use these names, because the self-check looks for them: prob_below_m4, prob_all_below,
prob_week_rules, prob_tomorrow.
n_week = 163 # aftershocks expected in the week, days 11 to 17
n_tomorrow = 29 # aftershocks expected tomorrow, day 11
prob_m4 = 13 / 637 # the probability that one aftershock is magnitude 4 or more (Your turn 3)
prob_below_m4 = 1 - prob_m4
prob_all_below = prob_below_m4 ** n_week
prob_week_rules = 1 - prob_all_below
prob_tomorrow = 1 - prob_below_m4 ** n_tomorrow
print("the week: all below magnitude 4 with probability", round(prob_all_below, 3), " at least one of 4 or more with probability", round(prob_week_rules, 3))
print("tomorrow: at least one of magnitude 4 or more with probability", round(prob_tomorrow, 3))
the week: all below magnitude 4 with probability 0.035 at least one of 4 or more with probability 0.965 tomorrow: at least one of magnitude 4 or more with probability 0.45
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert abs(prob_below_m4 - (1 - prob_m4)) < 1e-9 and abs(prob_all_below - (1 - prob_m4) ** n_week) < 1e-9, "prob_below_m4 is 1 minus prob_m4, and prob_all_below is prob_below_m4 raised to the power n_week"
assert abs(prob_week_rules - (1 - (1 - prob_m4) ** n_week)) < 1e-9, "prob_week_rules is 1 minus prob_all_below"
assert abs(prob_tomorrow - (1 - (1 - prob_m4) ** n_tomorrow)) < 1e-9, "prob_tomorrow is 1 minus prob_below_m4 raised to the power n_tomorrow"
print(f"✓ Your turn 5 — at least one of magnitude 4 or more: {round(100 * prob_tomorrow)}% tomorrow, {round(100 * prob_week_rules)}% in the week")
✓ Your turn 5 — at least one of magnitude 4 or more: 45% tomorrow, 97% in the week
The rules give 97% for the week; the simulation gave 96%. The two agree because they assume the same things: 163 aftershocks come, each is magnitude 4 or more with probability 0.02, and the sizes are independent. Tomorrow alone is close to an even chance, 45%.
When to use which. The rules are quicker. The simulation also gives the mean and the standard deviation, and it answers questions these rules do not, such as the probability of two or more (homework part 3).
Now compare with what came, day by day. The cell uses the rules for each of the seven days, with the count the model expects on that day, and prints the real number beside it.
came = quakes[(quakes["days"] >= 10) & (quakes["days"] < 17)] # the aftershocks that came in days 11 to 17: from 10 to 17 days after the mainshock
prob_m4 = 13 / 637 # the probability that one aftershock is magnitude 4 or more
for day in [11, 12, 13, 14, 15, 16, 17]:
n_day = 320 / day # the model's count for that day, with K = 320
prob_day = 1 - (1 - prob_m4) ** n_day # the three rules: at least one of magnitude 4 or more on that day
that_day = came[(came["days"] >= day - 1) & (came["days"] < day)] # the real aftershocks of that day
came_day = (that_day["mag"] >= 4).sum() # how many of them are magnitude 4 or more
print("day", day, " probability of at least one:", round(prob_day, 2), " came:", came_day)
day 11 probability of at least one: 0.45 came: 2 day 12 probability of at least one: 0.42 came: 0 day 13 probability of at least one: 0.4 came: 1 day 14 probability of at least one: 0.38 came: 0 day 15 probability of at least one: 0.36 came: 0 day 16 probability of at least one: 0.34 came: 0 day 17 probability of at least one: 0.32 came: 0
In the week: 97%, and 3 came. The nearly certain thing happened.
Tomorrow, day 11: 45%, and 2 came. One day cannot test a probability of 45%: it comes true on about 45 days in 100 and fails on the others. Seven days say a little more. Each day is one try with its own probability, so the expected number of days with at least one is the seven probabilities added up: 2.7, which is 2 or 3 of the 7 days. It came on 2.
Think. How would you check a forecast of "45%" properly? How many days like this one would you need?
Your forecast against the real week¶
| your forecast on day 10 | what came | |
|---|---|---|
| aftershocks in days 11 to 17 | 163 | 169 |
| of magnitude 3 or more | 54 | 62 |
| of magnitude 4 or more | 3.3, usually 2 to 5 | 3 |
| at least one of magnitude 4 or more | a 97% probability | yes |
The official forecast was close to yours. The US Geological Survey published its own on 15 July 2019, for nearly the same week: 23 to 52 aftershocks of magnitude 3 or more, and a probability above 99% of at least one. Your 54 is just above that range (the Survey does not count over exactly this box). Its forecast has the same two ingredients as yours: how fast the count falls, and what proportion of aftershocks are large.
Two limits of your forecast:
- It is rough. The curve runs above some dots and below others. The Survey renews its forecast as the sequence goes on.
- It cannot reach sizes the sample does not hold. The Survey also gave a 30% probability of a magnitude 5 or more. Your sample holds none (its largest is 4.9), so its proportion is 0, and that cannot be right. The transformations week shows how to reach sizes too rare to be in the sample.
One aftershock was larger than all of these: a magnitude 5.5, about 11 months after the mainshock.
The question, answered¶
Ten days after the Ridgecrest mainshock about 160 more aftershocks were due within a week and one in fifty aftershocks was magnitude 4 or more, so at least one more of that size was nearly certain, and 3 came.
Week 5 summary¶
The science¶
A big earthquake has just struck. What will the next week bring?
A large earthquake raises the probability of more earthquakes around it, and that extra rate dies away roughly as one over the time since it. A probability is a proportion of many tries: of what has happened, or of futures you simulate. A forecast is a count and a proportion turned into an expected number and a probability of at least one, and the week that comes can still differ from it.
- Probability. The proportion of many tries in which a thing happens.
- Omori's law. The number of aftershocks in a day falls roughly as one over the time since the mainshock: twice as late, about half as many.
- Expected number. How many to expect on average: the number of tries times the probability of each.
- Random experiment. Something you can do again and again whose result chance decides each time: toss a coin, roll a die.
- Monte Carlo. Let the computer live through the same situation thousands of times, each time with fresh random draws, and count how often a thing happens.
- Standard deviation. How far a value usually is from the average.
- Independent. Knowing how one turned out tells you nothing about the next.
- Aftershock sequence. A large earthquake changes the stress in the rock around it, and parts of the same fault and faults near it that were already close to failing then fail too — most within hours and days, fewer and fewer over the months that follow.
The code¶
| Function | What it does |
|---|---|
| Python | |
len(table) |
how many rows |
for x in things: |
do the same thing once for each item |
range(n) |
the numbers 0 to n-1, to count with |
list.append(x) |
add one item to the end of a list |
round(x, 1) |
trim a long decimal — 9.1 rather than sixteen digits |
| NumPy | |
array.std() |
the standard deviation: how far a value usually is from the average |
np.random.default_rng(seed) |
a random-number generator you can reproduce — the same seed gives the same draws |
rng.choice(values, n) |
n of the values drawn at random; the same one may be drawn more than once |
rng.random(n) |
n numbers drawn at random between 0 and 1 |
np.array(list) |
the container that does arithmetic to every number at once |
| pandas | |
mask.mean() |
the proportion of a mask that is True, in one line |
series.mean() |
the average value |
table[mask] |
keep only the rows where the mask is True |
| Matplotlib | |
plt.plot(x, y) |
join points with a line — here, the model K / day drawn through the daily counts |
Homework¶
Make the same forecast for the Bay Area's own earthquake, and check it against what came. At 5:04 in the afternoon of 17 October 1989 the magnitude 6.9 Loma Prieta earthquake struck the Santa Cruz Mountains, about 100 km south of campus. It is the earthquake that interrupted the 1989 World Series and broke a section of the Bay Bridge.
loma holds every event of magnitude 2.5 or more within 60 km of it, from the start of 1989
to a year after the mainshock: 644 rows, with the same days column, counted from its
own mainshock. The cell below shows the sequence.
The three parts follow class, with less given:
| part | question | as in class |
|---|---|---|
| 1 | How many aftershocks in days 11 to 17? | Your turn 2 |
| 2 | What is the probability of one of magnitude 4 or more? | Your turns 3 and 5 |
| 3 | Was the week that came an ordinary one? | Your turn 4 |
Each part ends by comparing your numbers with the real aftershocks. If you have restarted the kernel, run the setup cell at the top; nothing else from class is needed.
plt.figure(figsize=(7, 3.6))
plt.scatter(loma["days"], loma["mag"], s=5)
plt.xlim(-60, 60)
plt.xlabel("days since the magnitude 6.9")
plt.ylabel("magnitude")
plt.title("Loma Prieta: every event from two months before to two months after")
plt.show()
✏️ Your turn 6¶
Ten days after Loma Prieta, how many aftershocks should days 11 to 17 bring?
Why: it is the first half of the forecast, as in Your turn 2: a K found from the first ten days, then the model added up over the coming week.
- The first ten days.
loma_first_tenholds their aftershocks. The cell then counts them day by day. - Find K. Change
K_lomauntil the red curve passes through the dots from day 2 on. This fit is rougher than Ridgecrest's, because there are fewer aftershocks. - Forecast.
expected_week_lomais the model added up over days 11 to 17.
The last lines of the cell compare your forecast with what came. Choose K from the blue dots alone: on day 10 the orange ones had not happened yet.
Use these names, because the self-check looks for them: loma_first_ten, K_loma,
expected_week_loma.
# Step 1: the first ten days
loma_first_ten = loma[(loma["days"] > 0) & (loma["days"] < 10)]
loma_past_counts, edges = np.histogram(loma_first_ten["days"], bins=np.arange(0, 11)) # the count of each day
past_days = np.arange(1, 11) # the day numbers 1 to 10; loma_past_counts holds the count of each
# Step 2: find K
K_loma = 47
plt.figure(figsize=(7, 3.6))
plt.scatter(past_days, loma_past_counts, label="the count of each day")
plt.plot(past_days, K_loma / past_days, color="firebrick", label="the model, K / day")
plt.yscale("log")
plt.xlabel("day after the mainshock")
plt.ylabel("aftershocks in the day")
plt.title(f"Loma Prieta: Omori's law with K = {K_loma} (n = {len(past_days)} days)")
plt.legend()
plt.show()
# Step 3: forecast
coming_days = np.arange(11, 18) # the day numbers 11 to 17: the week to forecast
expected_week_loma = (K_loma / coming_days).sum()
print("K =", K_loma, " expected in days 11 to 17:", round(expected_week_loma, 1))
# What came (nothing to change below)
came_loma = loma[(loma["days"] >= 10) & (loma["days"] < 17)] # the aftershocks that came in days 11 to 17: from 10 to 17 days after the mainshock
came_counts_loma, edges = np.histogram(came_loma["days"], bins=np.arange(10, 18)) # the count of each of days 11 to 17
plt.figure(figsize=(7, 3.6))
plt.scatter(past_days, loma_past_counts, label="the count of each day, days 1 to 10")
plt.plot(past_days, K_loma / past_days, color="firebrick", label="the model, K / day")
plt.plot(coming_days, K_loma / coming_days, color="firebrick", ls="--", label="your forecast, days 11 to 17")
plt.scatter(coming_days, came_counts_loma, color="darkorange", label="what came")
plt.yscale("log")
plt.xlabel("day after the mainshock")
plt.ylabel("aftershocks in the day")
plt.title(f"Loma Prieta with K = {K_loma}: your forecast and what came (n = {len(came_loma)} came)")
plt.legend()
plt.show()
print("your forecast for each of days 11 to 17:", (K_loma / coming_days).round())
print("what came on each of them: ", came_counts_loma)
print("your forecast for the week:", round(expected_week_loma), " what came:", len(came_loma))
K = 47 expected in days 11 to 17: 24.0
your forecast for each of days 11 to 17: [4. 4. 4. 3. 3. 3. 3.] what came on each of them: [2 3 3 2 3 3 2] your forecast for the week: 24 what came: 18
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (loma_first_ten["days"] > 0).all() and (loma_first_ten["days"] < 10).all() and len(loma_first_ten) > 300, "loma_first_ten keeps the rows of loma with days above 0 and below 10"
assert 25 <= K_loma <= 70, "choose K_loma so that the curve passes through the dots from day 2 on: a value between about 25 and 70"
assert abs(expected_week_loma - (K_loma / coming_days).sum()) < 1e-9, "expected_week_loma is the model on the COMING days, K_loma / coming_days, added up"
print(f"✓ Homework 1 — with K = {K_loma}, about {round(expected_week_loma)} aftershocks are expected in days 11 to 17")
✓ Homework 1 — with K = 47, about 24 aftershocks are expected in days 11 to 17
What the numbers mean. With K = 47 the forecast for the week is 24, and 18 came: the right size, a little high. Day by day the model said 3 or 4 aftershocks and 2 or 3 came. The fit is rougher than Ridgecrest's. From day 2 on, the ten days hold only 107 aftershocks, a few a day, so chance moves each dot a long way: 10 on day 8, then 3 on day 9. Any K from about 25 to 70 passes among the dots, and the forecast changes with it.
✏️ Your turn 7¶
What was the probability of at least one aftershock of magnitude 4 or more in days 11 to 17?
Why: it is the second half of the forecast, as in Your turns 3 and 5: a probability measured as a proportion, an expected number, and the three rules.
sample_loma: the aftershocks from 1 day to 10 days after the mainshock, as in class and for the same two reasons.prob_m4_loma: the proportion of the sample that is magnitude 4 or more.expected_m4_loma: how many of magnitude 4 or more to expect among next week'sn_week_lomaaftershocks.prob_week_loma: the probability of at least one of them, by the three rules.
The count from part 1 is written in the cell, so this part stands on its own. The last lines of the cell compare your numbers with what came.
Use these names, because the self-check looks for them: sample_loma, prob_m4_loma,
expected_m4_loma, prob_week_loma.
n_week_loma = 24 # from part 1: aftershocks expected at Loma Prieta in days 11 to 17
sample_loma = loma[(loma["days"] >= 1) & (loma["days"] < 10)]
prob_m4_loma = (sample_loma["mag"] >= 4).mean()
expected_m4_loma = n_week_loma * prob_m4_loma
prob_week_loma = 1 - (1 - prob_m4_loma) ** n_week_loma
print((sample_loma["mag"] >= 4).sum(), "of", len(sample_loma), "in the sample are magnitude 4 or more: a probability of", round(prob_m4_loma, 2))
print("expected:", round(expected_m4_loma, 1), "of magnitude 4 or more")
print(f"the probability of at least one of magnitude 4 or more: {round(100 * prob_week_loma)}%")
# What came (nothing to change below)
came_loma = loma[(loma["days"] >= 10) & (loma["days"] < 17)] # the aftershocks that came in days 11 to 17: from 10 to 17 days after the mainshock
came_m4_loma = (came_loma["mag"] >= 4).sum() # how many of them are magnitude 4 or more
print("magnitude 4 or more: expected", round(expected_m4_loma, 1), " came", came_m4_loma)
print("proportion of magnitude 4 or more: in the sample", round(prob_m4_loma, 3), " in what came", round(came_m4_loma / len(came_loma), 3))
12 of 107 in the sample are magnitude 4 or more: a probability of 0.11 expected: 2.7 of magnitude 4 or more the probability of at least one of magnitude 4 or more: 94% magnitude 4 or more: expected 2.7 came 2 proportion of magnitude 4 or more: in the sample 0.112 in what came 0.111
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(sample_loma) == ((loma["days"] >= 1) & (loma["days"] < 10)).sum() and (sample_loma["days"] >= 1).all(), "sample_loma keeps every row of loma from 1 day to 10 days after the mainshock"
assert abs(prob_m4_loma - (sample_loma["mag"] >= 4).mean()) < 1e-9, "prob_m4_loma is the proportion of sample_loma with magnitude 4 or more: the mean of a mask"
assert abs(expected_m4_loma - n_week_loma * prob_m4_loma) < 1e-9, "an expected number is the number of tries, n_week_loma, times the probability"
assert abs(prob_week_loma - (1 - (1 - prob_m4_loma) ** n_week_loma)) < 1e-9, "prob_week_loma is 1 minus the probability that all n_week_loma are below magnitude 4"
print(f"✓ Homework 2 — a probability of {round(prob_m4_loma, 2)} for one aftershock, {round(expected_m4_loma, 1)} expected, and a {round(100 * prob_week_loma)}% probability of at least one")
✓ Homework 2 — a probability of 0.11 for one aftershock, 2.7 expected, and a 94% probability of at least one
What the numbers mean. 12 of the 107 aftershocks in the sample were magnitude 4 or more, a probability of 0.11: much higher than Ridgecrest's 0.02. So 2.7 were expected among the 24, with a 94% probability of at least one, and 2 came. The likely thing happened, and the proportion held: 0.111 of what came was magnitude 4 or more. The weak part is how little the proportion rests on: 12 earthquakes, and 6 of them are listed as exactly 4.0, because these magnitudes are given to one decimal. Had the question been "more than 4", the proportion would have halved.
✏️ Your turn 8¶
2 aftershocks of magnitude 4 or more came. Was that an ordinary week under your forecast?
Why: the rules of part 2 give only the probability of at least one. To say whether "2" is ordinary you need the whole distribution, and a simulation gives it, as in Your turn 4.
- Simulate the week 2,000 times. This time the loop is yours to write: the generator,
the
forline, the draw, and whatappendkeeps. One possible week isn_week_lomamagnitudes drawn fromsample_loma. mean_count_lomaandstd_count_loma: the mean and the standard deviation of the 2,000 counts.prob_two_or_more: the proportion of the weeks with 2 or more.
The count from part 1 is written in the cell. Its last lines draw your 2,000 weeks, with what came marked on them.
Use these names, because the self-check looks for them: sample_loma, counts_m4_loma,
mean_count_loma, std_count_loma, prob_two_or_more.
n_week_loma = 24 # from part 1: aftershocks expected at Loma Prieta in days 11 to 17
sample_loma = loma[(loma["days"] >= 1) & (loma["days"] < 10)]
rng = np.random.default_rng(1989)
counts_m4_loma = []
for run in range(2000):
one_week = rng.choice(sample_loma["mag"], n_week_loma)
counts_m4_loma.append((one_week >= 4).sum())
counts_m4_loma = np.array(counts_m4_loma)
mean_count_loma = counts_m4_loma.mean()
std_count_loma = counts_m4_loma.std()
prob_two_or_more = (counts_m4_loma >= 2).mean()
print("mean", round(mean_count_loma, 1), " standard deviation", round(std_count_loma, 1))
print("weeks with 2 or more: a proportion of", round(prob_two_or_more, 2))
# What came (nothing to change below)
came_m4_loma = 2 # what came: the aftershocks of magnitude 4 or more in days 11 to 17
plt.figure(figsize=(7, 3.6))
plt.hist(counts_m4_loma, bins=np.arange(0, counts_m4_loma.max() + 2), align="left") # one bar for each whole number
plt.axvline(mean_count_loma, color="firebrick", label="the mean")
plt.axvline(mean_count_loma - std_count_loma, color="firebrick", ls="--", label="one standard deviation each side")
plt.axvline(mean_count_loma + std_count_loma, color="firebrick", ls="--")
plt.axvline(came_m4_loma, color="darkorange", lw=3, label="what came")
plt.xlabel("aftershocks of magnitude 4 or more in the week")
plt.ylabel("number of weeks")
plt.title("Loma Prieta: 2,000 possible weeks (n = 2,000)")
plt.legend()
plt.show()
mean 2.7 standard deviation 1.5 weeks with 2 or more: a proportion of 0.78
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(sample_loma) == ((loma["days"] >= 1) & (loma["days"] < 10)).sum() and (sample_loma["days"] >= 1).all(), "sample_loma keeps every row of loma from 1 day to 10 days after the mainshock, as in part 2"
assert len(counts_m4_loma) == 2000 and counts_m4_loma.max() < 25, "counts_m4_loma holds one number for each of the 2000 weeks: how many of that week's magnitudes are 4 or more"
assert 2.48 < counts_m4_loma.mean() < 2.9 and 1.31 < counts_m4_loma.std() < 1.78, "each week is n_week_loma magnitudes drawn from sample_loma's mag column, and counts_m4_loma keeps the COUNT of those that are 4 or more: about 2.7 on average"
assert abs(mean_count_loma - counts_m4_loma.mean()) < 1e-9 and abs(std_count_loma - counts_m4_loma.std()) < 1e-9, "mean_count_loma is counts_m4_loma.mean() and std_count_loma is counts_m4_loma.std()"
assert abs(prob_two_or_more - (counts_m4_loma >= 2).mean()) < 1e-9, "prob_two_or_more is the proportion of weeks where counts_m4_loma is 2 or more: the mean of a mask"
print(f"✓ Homework 3 — mean {round(mean_count_loma, 1)}, standard deviation {round(std_count_loma, 1)}; 2 or more in {round(100 * prob_two_or_more)}% of the weeks")
✓ Homework 3 — mean 2.7, standard deviation 1.5; 2 or more in 78% of the weeks
What the numbers mean. The 2,000 possible weeks hold 2.7 aftershocks of magnitude 4 or more on average, with a standard deviation of 1.5: usually between 2 and 4. 78% of the weeks hold 2 or more, and the orange line stands among the tallest bars. So the 2 that came were an ordinary week under the forecast. The simulation also agrees with the rules of part 2: 95% of its weeks hold at least one, against 94%.