Clustering · Part A: the methods, on a planted truth

An earthquake catalogue lists where and when each event happened, and nothing else. Which events belong together — on one fault, in one aftershock sequence, in one swarm — and how would you know that a group you found is real?

Session 1, Regression & uncertainty, fitted models to data that came with answers: every training example had a label, and the error could be measured against it. A catalogue has no such labels. No one has marked which fault each earthquake broke, or which events are aftershocks of which. Grouping events without labels is clustering, the main kind of unsupervised learning: the method must decide from the features alone which events are similar enough to share a group.

Without labels, no answer is correct in itself. Every clustering method makes three choices for you, and the answer is only as good as those choices:

1. What counts as close? Which features describe an event (position, time, mechanism) and how the distance between two events is measured. When the features have different units, this includes how many kilometres a day or a degree is worth (A.9).
2. What counts as a group? A cluster can be defined by a centre, by a probability model, by a chain of close neighbours, or by a region of high density. The four definitions give four different answers on the same data.
3. At what scale? Every method has one parameter that sets how large or how many the clusters are, and the data do not settle it.

And one test that every result must pass: is it real? Every method returns clusters, even from random data.

The questions this notebook answers

1. What does a clustering method decide for you when the data have no labels? (A.1)
2. How do the four definitions of a cluster — centre, probability, connection and density — work, and where does each fail? (A.2–A.6)
3. Every method has a parameter the data do not settle. How do you choose it? (A.7)
4. Every method finds clusters in random data. How do you know a cluster is real? (A.8)
5. When origin time joins position, what does "close" mean? (A.9)
6. How do seismologists group earthquakes by space, time and magnitude at once, and tell an aftershock sequence from a swarm? (A.10)

All examples here are synthetic: we plant the groups, so every method's answer can be scored against the truth. On a real catalogue there is no truth to score against, which is why the methods are learned here first. Part B applies the same methods to the Geysers focal-mechanism catalogue, where the questions of scale and reality become the whole problem.

A.1 · A map with planted groups

We build one map of epicentres, in kilometres, and use it for every method in A.2–A.8. It holds five kinds of structure, each chosen because some method handles it badly:

• two long, thin faults, parallel and about two kilometres apart, the way seismicity lines up on a fault seen from above;
• a dense, compact patch, like a cluster of induced seismicity around one injection well;
• a broad, sparse cluster, a group of the same size in events but spread over a much larger area, so its density is far lower;
• a thin line of events joining the patch to the sparse cluster, such as a minor structure between two active volumes;
• background events, scattered at random over the whole map, which belong to no group.

The four groups are the truth a method should find. The thin line and the background belong to none of them; we give them the label $-1$, "no cluster", which is the convention every method below uses for events it leaves out.

Each event also carries a strike, the compass direction (in degrees clockwise from north) of the fault plane it broke. Events on a fault share that fault's strike, give or take some scatter; events in the patch share another; background events have random strikes. The methods never see it. We keep it back to test, in A.8, whether clusters found from position alone are coherent in something they were not built from.

1045 events: group 0: 300, group 1: 150, group 2: 200, group 3: 120, no group: 275
the two faults are 2.1 km apart

The helper below draws a map the same way for every method: events in no cluster in light grey, every cluster in its own colour, east and north at the same scale so that a round cluster looks round.

A method succeeds on this map if it finds the four coloured groups and leaves the grey events out. We score every method with two numbers that say exactly that.

• Groups found. A planted group counts as found when one cluster holds most of its events, and most of that cluster's planted events belong to the group, so that two groups merged into one cluster count as neither being found. Background events inside the cluster are ignored here; they are counted by the second number.
• Background in a cluster. The share of the grey events that the method put into some cluster. A method that must assign every event scores 100%.

The two thresholds that define "most" are printed below.

a group is found if one cluster holds at least 70% of it, and less than 30% of that cluster's planted events come from other groups
the truth itself: 4 of 4 groups found, 0% of background in a cluster

A.2 · A cluster as a centre: k-means

The simplest definition: a cluster is the set of events closest to one centre. Choose a number of clusters $k$. Each cluster $j$ has a centre $\mu_j$, a point on the map, and each event $x_i$ belongs to the cluster whose centre is nearest. Write $z_i$ for the cluster event $i$ is assigned to. A good set of centres keeps every event close to its own centre, so k-means minimises

$$J = \sum_{i=1}^{N} \lVert x_i - \mu_{z_i} \rVert^2 ,$$

the sum, over all $N$ events, of the squared distance from each event to its centre.

Finding the best centres and the best assignment together is hard, but each half is easy when the other is fixed. Lloyd's algorithm (Lloyd 1982) alternates the two halves; MacQueen (1967) gave the method its name, k-means:

1. Assign: with the centres fixed, put each event in the cluster of its nearest centre. No other assignment gives a smaller $J$.
2. Update: with the assignment fixed, move each centre to the mean of its events. The mean is the point that minimises the sum of squared distances to a set of points, so again $J$ cannot grow.

Repeat until no event changes cluster. Because neither step can increase $J$ and there are only finitely many assignments, the loop always stops. It is the whole algorithm:

k = 4, seed 0: J = 3178 km²; 2 of 4 groups found, 100% of background in a cluster

The figure follows one run: the starting centres (four events chosen at random), the state after a single assign-and-update step, and the converged answer. Black crosses are the centres.

Two properties of k-means are visible in the converged map, and both follow from its definition:

• Every event is put in a cluster. The background events are shared out among the four groups, because "no cluster" is not an option.
• Clusters are convex tiles of similar width. The boundary between two clusters is the line halfway between their centres, so the map is cut into straight-edged tiles, each convex: the straight line between any two points of a tile stays inside it, so a tile cannot bend. A long fault is cut across into pieces, and the two parallel faults can share a tile.

Where it starts can matter

Lloyd's algorithm only ever moves downhill in $J$, so it stops at the first point from which no single step helps: a local minimum, which need not be the lowest $J$ there is. Different starting centres can stop in different places. The standard remedy is to run from many starting points and keep the run with the lowest $J$; scikit-learn's KMeans does this (its n_init argument), and also spreads the starting centres apart (k-means++), which makes a bad start less likely. On this map, with $k = 4$, almost every start reaches the same answer:

our code: lowest J = 3178 km², 18 of 20 seeds within 1 km² of it; highest J = 4395 km²
scikit-learn, n_init = 20:  J = 3178 km²

So the poor score against the truth is not bad luck in the start. It is what k-means is: the best set of four round tiles still cuts the faults across and shares out the background. A better start finds a lower $J$; it cannot change what k-means can represent. The next method can represent more, and there the start will matter a great deal.

A.3 · A cluster as a probability distribution: the Gaussian mixture

k-means has hard edges and treats every cluster as round and equally wide. A Gaussian mixture model (GMM) relaxes all three. It is a model of how the data were generated:

1. Each event first picks a cluster $j$, with probability $\pi_j$ (the weights, which sum to 1).
2. Its position is then drawn from that cluster's 2-D Gaussian (normal) distribution, with a centre $\mu_j$ and a $2 \times 2$ covariance matrix $\Sigma_j$. The covariance sets the shape of the cluster: how wide it is in each direction, and how it is oriented. Its contours of equal probability are ellipses.

The cluster each event picked, $z_i$, is never observed. A quantity the model needs but the data do not show is a latent variable. Summing over the cluster an event might have come from gives the probability density of one position,

$$p(x) = \sum_{j=1}^{k} \pi_j \, \mathcal{N}(x \mid \mu_j, \Sigma_j),$$

where $\mathcal{N}(x \mid \mu, \Sigma)$ is the Gaussian density. We fit the weights, centres and covariances by maximum likelihood, as in Session 1: choose them to make the observed events as probable as possible, that is, to maximise the log-likelihood $\sum_i \log p(x_i)$.

Expectation–maximisation

If we knew each event's cluster, the fit would be easy: each Gaussian would be fitted to its own events (their mean and their covariance), and each weight would be the fraction of events in that cluster. We do not know the clusters, but given a current guess of the parameters we can compute how probable each cluster is for each event. Expectation–maximisation (EM; Dempster, Laird & Rubin 1977) alternates these two steps, exactly as Lloyd's algorithm did:

E-step. With the parameters fixed, compute the responsibility of cluster $j$ for event $i$: the probability that the event came from that cluster, by Bayes' rule,

$$\gamma_{ij} = \frac{\pi_j \, \mathcal{N}(x_i \mid \mu_j, \Sigma_j)} {\sum_{l=1}^{k} \pi_l \, \mathcal{N}(x_i \mid \mu_l, \Sigma_l)} .$$

Each event's responsibilities sum to 1 over the clusters. They are a soft assignment: an event between two clusters can belong half to each.

M-step. With the responsibilities fixed, refit each cluster to all events, each weighted by its responsibility. Write $N_j = \sum_i \gamma_{ij}$ for the effective number of events in cluster $j$. Then

$$\pi_j = \frac{N_j}{N}, \qquad \mu_j = \frac{1}{N_j}\sum_i \gamma_{ij}\, x_i, \qquad \Sigma_j = \frac{1}{N_j}\sum_i \gamma_{ij}\,(x_i - \mu_j)(x_i - \mu_j)^\top .$$

These are the fraction of events, the mean and the covariance — the same formulas as with known clusters, with each event counted in proportion to its responsibility.

Each EM step can only raise the log-likelihood (the reason is the Optional subsection below), so EM, like Lloyd's algorithm, climbs to a local maximum that depends on where it starts.

The code computes the Gaussian densities in logarithms, because the density of a distant event under a narrow Gaussian is too small to store as a number; logsumexp adds probabilities that are stored as logarithms.

seed 0:  76 EM iterations, log-likelihood rose every step: True, final log-likelihood = -4182
seed 6:  60 EM iterations, log-likelihood rose every step: True, final log-likelihood = -3871
seed 2:  35 EM iterations, log-likelihood rose every step: True, final log-likelihood = -3625

How many clusters? The map has four planted groups and a background that fills the whole map. A mixture has no "no cluster" label, but it can spend one broad Gaussian on the background. So we fit $k = 5$, and A.7 asks how $k$ could be chosen without knowing the truth.

Where EM starts matters

EM climbs the likelihood the way Lloyd's algorithm descended $J$, and stops at a local maximum. With more freedom per cluster there are more places to stop. The three runs above started from three seeds; each ellipse below is one cluster's Gaussian, drawn at two standard deviations.

seed 0: 2 of 4 groups found, 100% of background in a cluster
seed 6: 2 of 4 groups found, 100% of background in a cluster
seed 2: 4 of 4 groups found, 100% of background in a cluster

The starts end in different answers, and the run with the highest likelihood is also the one that finds the most groups. The likelihood needs no labels, so on real data too the rule is: run EM from many starts and keep the highest likelihood (scikit-learn's GaussianMixture does this with n_init). Nothing guarantees that the highest of a few starts is the highest there is.

We keep the best run, seed 2. The map below shows what k-means could not: how sure the model is about each event. Colour is the event's largest responsibility; dark events are shared between clusters.

seed 2: 4 of 4 groups found, 100% of background in a cluster
events with largest responsibility below 0.9: 14%

Because each covariance can stretch and turn, an ellipse can follow a fault along its length, which a round k-means tile cannot. The cost is that every event is still assigned somewhere: the broad component takes the background, and also the sparse cluster's edges and anything else that fits no narrow ellipse.

k-means is a Gaussian mixture with its freedom removed

Restrict the mixture in three ways: every covariance is the same circle, $\Sigma_j = \sigma^2 I$; the weights are equal; and $\sigma$ shrinks towards zero. Then in the E-step the nearest centre's Gaussian is larger than every other by a factor that grows without bound, so its responsibility goes to 1 and all others to 0. The soft assignment becomes the hard nearest-centre assignment, and the M-step becomes the plain mean. That is Lloyd's algorithm: k-means is EM for a mixture of identical, infinitely narrow round Gaussians. Everything the mixture adds — shape, orientation, unequal sizes, soft membership — is a restriction of k-means removed.

The shape restriction can also be removed by degrees. scikit-learn's GaussianMixture offers four covariance_type settings: spherical (round clusters, each its own width), diag (ellipses aligned with the axes), tied (one shared ellipse) and full (any ellipse). A fault that strikes obliquely needs full:

spherical: log-likelihood = -4210, 2 of 4 groups found, 100% of background in a cluster

full     : log-likelihood = -3628, 4 of 4 groups found, 100% of background in a cluster

Optional: why each EM step raises the likelihood

Take one event and drop the index $i$. Its log-likelihood is $\log \sum_j \pi_j \mathcal{N}(x \mid \mu_j, \Sigma_j)$. Write $\theta$ for all the parameters together, and let $q_j$ be any set of probabilities over the clusters ($q_j \ge 0$, $\sum_j q_j = 1$).

Multiply and divide each term by $q_j$:

$$\log p(x \mid \theta) = \log \sum_j q_j \,\frac{\pi_j \mathcal{N}(x \mid \mu_j, \Sigma_j)}{q_j}.$$

The right side is the logarithm of an average (weights $q_j$). The logarithm is concave, so the logarithm of an average is at least the average of the logarithms (Jensen's inequality):

$$\log p(x \mid \theta) \;\ge\; \sum_j q_j \log \frac{\pi_j \mathcal{N}(x \mid \mu_j, \Sigma_j)}{q_j} \;\equiv\; \mathcal{L}(q, \theta).$$

$\mathcal{L}$ is a lower bound on the log-likelihood for every choice of $q$; it is called the evidence lower bound (ELBO). EM raises this bound in turn over each of its two arguments:

• E-step: with $\theta$ fixed, choose $q$ to make the bound as high as possible. The bound equals the log-likelihood (Jensen's inequality becomes an equality) when the ratio inside the logarithm is the same for every $j$, that is, when $q_j \propto \pi_j \mathcal{N}(x \mid \mu_j, \Sigma_j)$. Normalised, that is exactly the responsibility $\gamma_j$. After the E-step the bound touches the log-likelihood.
• M-step: with $q$ fixed, choose $\theta$ to make the bound as high as possible. Only the numerator depends on $\theta$, so this maximises $\sum_j \gamma_j \log \pi_j \mathcal{N}(x \mid \mu_j, \Sigma_j)$, a weighted Gaussian fit, whose solution is the M-step's three formulas.

The M-step raises the bound; the bound is below the log-likelihood; and before the M-step the two were equal. So the log-likelihood after the M-step is at least the bound after it, which is at least the log-likelihood before it. No EM step can lower the likelihood. Nothing stops it at a local maximum, however, which is why the seed matters.

A.4 · A cluster as a chain of neighbours: hierarchical clustering

k-means and the mixture both start from a number of clusters and a shape. Hierarchical clustering starts from neither. It builds the whole family of groupings at once:

1. Start with every event in a group of its own.
2. Find the two groups that are closest, and merge them.
3. Repeat until one group is left.

The record of the merges is a tree, the dendrogram: each merge is drawn at the height of the distance between the two groups it joined. Cutting the tree at some height gives a grouping; cutting low gives many small groups, cutting high gives few large ones.

Step 2 needs a distance between two groups, not two events, and there are several choices, each called a linkage. For groups $A$ and $B$:

• single: the distance between their closest pair of events, $\min_{a \in A, b \in B} d(a, b)$;
• complete: the distance between their farthest pair, $\max_{a, b} d(a, b)$;
• average: the mean distance over all pairs;
• Ward (Ward 1963): how much the merge would raise k-means' objective $J$, the summed squared distance of events to their group's mean. Merging $A$ and $B$, with $n_A$ and $n_B$ events and means $\mu_A$ and $\mu_B$, raises $J$ by $\frac{n_A n_B}{n_A + n_B}\lVert \mu_A - \mu_B \rVert^2$.

The linkage decides what kind of cluster the tree prefers. Single linkage joins two groups as soon as any two of their events are close, so it follows chains; Ward linkage keeps groups compact, like k-means. We build both trees for the map and cut each into five groups.

single: group sizes [1006, 23, 7, 5, 4]; 0 of 4 groups found, 100% of background in a cluster
ward  : group sizes [332, 275, 174, 146, 118]; 2 of 4 groups found, 100% of background in a cluster

Single linkage puts nearly every event in one group and leaves a few isolated background events as the other "groups". Background events are spread everywhere, so there is always a chain of them, each close to the next, from any group to any other; and the thin line joins the patch to the sparse cluster directly. This failure is called chaining. Ward linkage finds compact groups, but it inherits k-means' habit of putting every event somewhere.

The dendrogram shows where to cut: a long vertical gap between successive merges means that the next merge joins groups much farther apart than any merge before it. A.7 uses that gap to choose the number of groups.

Hierarchical clustering needs the distance between every pair of events. For this map that is trivial; for a full catalogue it is not:

   1045 events: 5.45e+05 pairs,    0.004 GB of distances
  46395 events: 1.08e+09 pairs,    8.610 GB of distances
 403475 events: 8.14e+10 pairs,  651.167 GB of distances

The two larger numbers are the sizes of the Geysers catalogues used in Part B: the focal-mechanism catalogue and the full location catalogue. So on large catalogues hierarchical clustering is run only on a sparse graph that connects each event to its nearest neighbours (scikit-learn's AgglomerativeClustering with connectivity=), which is what Part B does.

A.5 · A cluster as a dense region: DBSCAN

Look at the map again. What makes the planted groups groups is that events are crowded there, much more than in the background. DBSCAN (density-based spatial clustering of applications with noise; Ester et al. 1996) turns that into a definition. It has two parameters, a radius $\varepsilon$ (eps) and a count min_samples:

1. An event is a core event if at least min_samples events, itself included, lie within distance $\varepsilon$ of it. Core events sit where the map is crowded.
2. Two core events within $\varepsilon$ of each other are in the same cluster; following such links from core event to core event gives the whole cluster.
3. An event that is not core but lies within $\varepsilon$ of a core event is a border event, and joins that core event's cluster.
4. Every other event is noise: label $-1$.

Step 4 is new. k-means, the mixture and Ward linkage all had to put every event in a cluster; DBSCAN can say that an event belongs to none, which is what the background is. The whole algorithm, with a neighbour search from scikit-learn and the linking from scipy:

eps = 0.3 km, min_samples = 10: 4 clusters, 649 core, 30 border, 366 noise events
same labels as scikit-learn's DBSCAN for every event: True

The two parameters are one density

A core event has at least min_samples events in a circle of area $\pi \varepsilon^2$, so the two parameters together set one threshold, a density in events per square kilometre:

$$\rho_{\min} = \frac{\texttt{min\_samples}}{\pi \varepsilon^2}.$$

Doubling $\varepsilon$ and multiplying min_samples by four keeps the same threshold (in three dimensions, by eight). A region is a cluster if it is denser than $\rho_{\min}$, and that is the whole question DBSCAN asks. So we measure how dense each planted group is: for every event, the local density is the number of events within its tenth-nearest neighbour, divided by the area of that circle.

fault A        : median density     130 events/km²
fault B        : median density     108 events/km²
dense patch    : median density     474 events/km²
sparse cluster : median density      12 events/km²
background     : median density       3 events/km²
densest group / sparsest group: 38x
DBSCAN eps = 0.3 km, min_samples = 10: threshold  35.4 events/km²
DBSCAN eps = 0.6 km, min_samples = 10: threshold   8.8 events/km²
DBSCAN eps = 0.8 km, min_samples = 10: threshold   5.0 events/km²

The densest group is dozens of times denser than the sparsest, and the sparsest is only a few times denser than the background. One threshold must fall somewhere in that range. We run DBSCAN at the three thresholds just printed. Dark events are core, light events border, grey noise.

eps = 0.3 km: 4 clusters; 3 of 4 groups found, 7% of background in a cluster
eps = 0.6 km: 5 clusters; 4 of 4 groups found, 38% of background in a cluster
eps = 0.8 km: 2 clusters; 0 of 4 groups found, 52% of background in a cluster

At the high threshold (small $\varepsilon$) DBSCAN finds the faults and the patch cleanly and calls the sparse cluster noise: it is not dense enough. At the middle threshold it finds all four, at the price of admitting more background events. At the low threshold the patch, the thin line and the sparse cluster become one cluster, and the two faults join into one. So DBSCAN can work, but only inside a narrow window of $\varepsilon$, and nothing in the data tells you where that window is; on this map we know only because we planted the truth. Seismicity makes the window narrower still: a swarm can be far denser than the slow seismicity around it, and one threshold then cannot suit both.

A.6 · Every density at once: HDBSCAN

HDBSCAN (hierarchical DBSCAN; Campello, Moulavi & Sander 2013) removes the single threshold. It runs DBSCAN at every $\varepsilon$ at once, keeps the result as a tree, and then picks from the tree the clusters that last longest. Four steps:

1. Core distance. For each event, the distance to its min_samples-th nearest neighbour, $c_i$. It is small where the map is crowded: $c_i$ is the smallest $\varepsilon$ at which event $i$ would be a core event.
2. Mutual reachability distance. Replace the distance between two events by $d_{\text{mr}}(i, j) = \max\{c_i,\; c_j,\; d(i, j)\}$. Inside a dense group this changes little; an isolated event is pushed away from everything, because its own core distance is large. So a chain of background events can no longer join two groups cheaply — the cure for chaining.
3. The tree. Build the single-linkage tree on $d_{\text{mr}}$. Cutting it at height $\varepsilon$ gives the clusters DBSCAN would find with that $\varepsilon$ and min_samples. Moving the cut down from large $\varepsilon$ to small, clusters shrink and split, and events drop out of them as noise.
4. Select. A split counts only if both parts have at least min_cluster_size events; a smaller part is just events dropping out. HDBSCAN measures density as $\lambda = 1/\varepsilon$, so moving down the tree is moving up in density. Each cluster is born at the $\lambda$ where it splits off its parent; each of its events stays until it drops out or the cluster itself splits. The cluster's stability is the sum, over its events, of how far in $\lambda$ each event stays: a cluster that holds many events over a wide range of density is stable. Working up from the bottom of the tree, a parent is kept instead of its children when its stability is larger than theirs combined; otherwise the children are kept (the excess of mass rule, cluster_selection_method="eom"). The alternative, "leaf", keeps the smallest clusters at the bottom of the tree.

Steps 1–3 can be written in a few lines, and they show what the tree holds. We cut it at three heights, which is DBSCAN at three thresholds, and drop any group smaller than min_cluster_size.

cut at 0.8 km: 2 clusters of at least 30 events, 18% noise
cut at 0.4 km: 3 clusters of at least 30 events, 37% noise
cut at 0.2 km: 4 clusters of at least 30 events, 52% noise

HDBSCAN: 4 clusters; 4 of 4 groups found, 31% of background in a cluster

High in the tree (left) the two faults are still one cluster, and the patch, the thin line and the sparse cluster another. Lower down the faults separate, but the sparse cluster breaks up and drops out; lower still, the long fault breaks up too. No single cut gives all four groups. HDBSCAN's selection does not cut at one level: the sparse cluster never splits into two parts large enough to count, so its events only drop out one by one and it stays one cluster; the merged faults split into two children that each hold many events over a wide range of density, so the children are kept. The result is all four groups. The thin line goes with the sparse cluster, which it touches; the scattered background is left out.

HDBSCAN still has parameters, but min_cluster_size has a meaning a seismologist can decide without trying values: the smallest group of events you would call a cluster. Its real advantage over DBSCAN shows in A.7: where DBSCAN works only in a narrow window of $\varepsilon$, HDBSCAN gives the same answer over a wide range of min_cluster_size.

Five methods, one map

The five definitions side by side, each at the setting used above:

k-means, k = 4          : 2 of 4 groups found, 100% of background in a cluster
Gaussian mixture, k = 5 : 4 of 4 groups found, 100% of background in a cluster
Ward linkage, 5 groups  : 2 of 4 groups found, 100% of background in a cluster
DBSCAN, eps = 0.6 km    : 4 of 4 groups found, 38% of background in a cluster
HDBSCAN, min size 30    : 4 of 4 groups found, 31% of background in a cluster

Two lessons hold beyond this map.

• Can the method leave events out? k-means, the mixture and Ward linkage are partitions: every event is placed in some cluster, so all of the background lands in clusters. The mixture finds all four groups only because it spent a fifth Gaussian on the background. DBSCAN and HDBSCAN can label an event noise. A catalogue always contains background seismicity that belongs to no structure, so for seismicity this difference matters more than any parameter.
• The definition decides the answer. The map did not change between panels; only the definition of a cluster did. The density methods do well here because their definition — regions denser than their surroundings — matches how this map was built. A map made only of Gaussian blobs with no background would favour the mixture. On real data the right definition is a statement about the Earth, and it has to be argued.

A.7 · Choosing the parameter without the truth

Every method so far had one parameter that we set by looking at the planted truth: $k$ for k-means and the mixture, where to cut the tree, $\varepsilon$ for DBSCAN, min_cluster_size for HDBSCAN. On a real catalogue there is no truth to look at. Each method comes with a label-free diagnostic meant to choose its parameter from the data alone. Here, where the truth is known, we can check whether each diagnostic chooses well.

Silhouette (for k-means; Rousseeuw 1987). For one event, let $a$ be its mean distance to the other events in its own cluster, and $b$ its mean distance to the events of the nearest other cluster. Its silhouette is

$$s = \frac{b - a}{\max(a, b)} .$$

$s$ is near 1 when the event sits deep inside its cluster, near 0 when it lies between two clusters, and negative when it is closer to another cluster than its own. The mean silhouette over all events measures how well separated the clusters are; choose the $k$ that maximises it.

BIC (for the mixture). A mixture with more Gaussians always fits better, so the likelihood alone would always choose more. The Bayesian information criterion (Schwarz 1978) charges for every parameter:

$$\text{BIC} = -2 \log L + p \log N ,$$

where $\log L$ is the maximised log-likelihood, $p$ the number of fitted parameters (a 2-D full-covariance mixture has 2 for each mean, 3 for each covariance and $k - 1$ free weights) and $N$ the number of events. Choose the $k$ with the lowest BIC.

The k-distance plot (for DBSCAN). Sort the events by the distance to their min_samples-th neighbour, from small to large, and plot that distance against rank. Events inside clusters have small distances and background events large ones, so the sorted curve bends upward where clusters end and background begins. The usual rule sets $\varepsilon$ at that bend, the knee; here it is found as the point of the curve farthest from the straight line joining its ends.

Stability across the parameter (for HDBSCAN). min_cluster_size is the smallest group you would call a cluster, so it is set from what you care about rather than from a curve. What can be checked is whether the answer changes when it changes.

Each top panel below is a diagnostic; each bottom panel is how many planted groups the method found at the same setting. The dashed line is the diagnostic's choice.

k-means: silhouette chooses k = 7; 3 of 4 groups found, 100% of background in a cluster
mixture: BIC chooses k = 8; 4 of 4 groups found, 100% of background in a cluster
DBSCAN:  the knee chooses eps = 0.53 km; 4 of 4 groups found, 29% of background in a cluster
HDBSCAN (eom): 8 of 9 values of min_cluster_size from 5 to 160 find all 4 groups

Read the columns one at a time.

• Silhouette rises and falls irregularly, and its peak picks a $k$ that misses a group. The bottom panel shows the deeper problem: k-means finds all four groups at no value of $k$, so no diagnostic could choose a good one. A diagnostic chooses among the answers a method can give; it cannot repair the method's definition of a cluster.
• BIC drops sharply once the mixture has enough components for the groups and the background, then stays almost flat. Every $k$ on the flat part finds all four groups, so here BIC works. This map is close to what a mixture assumes: compact blobs plus one broad background. On the Geysers catalogue, which is not, BIC keeps falling as $k$ grows and never chooses (Part B).
• The knee of the k-distance plot (the dashed line, carried down to the eps axis) lands inside DBSCAN's narrow good window. It works because this map has one clear background density, and the knee marks where the background begins. On the Geysers catalogue, whose density varies smoothly rather than in steps, the knee puts the whole field into one cluster (Part B).
• HDBSCAN with the default selection gives the same answer over a wide range of min_cluster_size, and that stability is itself the evidence. (Its other parameter, min_samples, is held at 10 throughout; it sets how smooth the density estimate is.) With leaf selection, which keeps the smallest clusters at the bottom of the tree, small values split the faults into pieces; leaf is the setting for finding the finest dense patches, which is how Part B uses it.

The general lesson: a diagnostic answers a narrow question (how separated are the clusters? how many parameters does the fit deserve? where does the background begin?), and it chooses well only when the data match what it assumes. On a planted map you can see whether they match; on a catalogue you cannot. So use diagnostics, but also check how much the answer changes across the parameter, and prefer a parameter with a physical meaning you can defend.

One map is one random draw. Does each conclusion survive a different draw of the same structures? We rebuild the map from ten more seeds and repeat the three density checks.

seed  DBSCAN eps finding all 4 (km)  knee: groups   HDBSCAN groups at min size 15 20 30 50 80
 207                     0.50-0.65             4    4   4   4   4   4
   1                     0.55-0.60             4    3   3   3   3   3
   2                     0.55-0.60             2    3   3   3   3   3

   3                     0.50-0.70             4    4   4   4   4   4
   4                     0.50-0.60             4    4   4   4   4   4
   5                     0.60-0.65             4    4   4   4   4   4
   6                     0.55-0.65             4    4   4   4   4   4

   7                     0.60-0.70             4    4   4   4   4   4
   8                     0.55-0.70             4    3   4   4   4   4
   9                     0.55-0.65             4    3   4   4   4   4
  10                     0.50-0.55             4    4   4   4   4   4

On every map DBSCAN's good window is narrow, and the knee finds it on all but one. HDBSCAN always finds the two faults and the patch, and gives nearly the same answer across min_cluster_size; whether it counts the sparse cluster as found depends on the draw, because that cluster sits close to the 70% line. That is a first hint of what A.8 shows directly: the sparse cluster is the weakest structure on the map.

For hierarchical clustering the diagnostic is the dendrogram's largest gap between successive merges (A.4): cutting inside the largest gap separates groups that are far apart.

ward  : the largest gap gives 2 groups; 0 of 4 groups found, 100% of background in a cluster
single: the largest gap gives 1 group; 0 of 4 groups found, 100% of background in a cluster

Ward's largest gap splits the map into two halves; single linkage's gives the one chained cluster. Neither diagnostic rescues a definition that does not fit the data.

A.8 · Is it real?

Every result so far was judged against a planted truth. On a catalogue there is none, so we need tests that need no labels. The first fact to face: every method finds clusters in random data. We draw the same number of events uniformly at random over the same map — no structure at all — and cluster them.

HDBSCAN on random events: 5 clusters, 36% noise

k-means cuts the random map into four tidy tiles, as it must. HDBSCAN finds clusters too: in a uniform scatter some regions are denser than others by chance, and HDBSCAN's default selection may not return the whole map as one cluster (unless allow_single_cluster=True), so it splits those chance fluctuations into clusters. Nothing in either output warns that there is nothing there. Three tests turn this into practice.

Test 1: a null catalogue

Build catalogues that keep everything except the structure you are testing, cluster each in exactly the same way, and compare. The comparison is then between the data and what the method finds when the answer is "nothing". Which properties a null catalogue keeps is the scientific choice. Here we keep the number of events and the map, and destroy all structure. Part B keeps the large-scale shape of the seismicity and moves the events instead: by each event's own location error, to see whether the clusters survive the uncertainty of the locations, and by about half a kilometre, to see whether they are more than chance lumps.

What to compare? The number of clusters is a poor statistic, since random catalogues give clusters too. A better one asks how strong each cluster is: its density contrast, the median local density of its events divided by the map's average density (events per square kilometre over the whole map). A cluster is significant if its contrast is higher than that of every cluster found in the null catalogues.

HDBSCAN clusters per catalogue: map 4; 20 null catalogues 2 to 7
density contrast of the 91 null clusters: median 1.3, largest 2.1
map cluster 0: 319 events, density contrast  17.1 -> significant
map cluster 1: 160 events, density contrast  14.4 -> significant
map cluster 2: 136 events, density contrast   1.9 -> not distinguishable from chance
map cluster 3: 208 events, density contrast  62.0 -> significant
k-means silhouette (k = 4): map 0.50; null catalogues 0.41 ± 0.01

Counting clusters would have failed: the null catalogues give as many as the map. Density contrast separates them. The two faults and the patch are many times denser than any chance cluster. The sparse cluster is not: a uniform catalogue of the same size produces groups as dense by chance. It is real — we planted it — but these data cannot prove it. A null test protects you from false clusters, and the price is that weak real ones may not pass. Whether a cluster passes also depends on the null. A null with the density of the background around the sparse cluster, which is lower than the map's average, would produce weaker chance clusters, and the sparse cluster might then pass; but that null assumes we already know where the background is. Choosing the null is part of the science.

The same idea rescues the silhouette. k-means always returns four clusters, but the map's silhouette is well above the range of the null catalogues, so the map has more separation than random events. The silhouette alone said little; the silhouette against its null distribution is a test.

Test 2: stability

A real structure should be found again from different data of the same kind. Split the events at random into two halves, cluster each half separately, and ask whether the two answers agree. To compare them on the same events, every event of the full map is labelled by its nearest event in each half. The agreement is measured by the adjusted Rand index (ARI; Hubert & Arabie 1985):

• Take every pair of events, and ask whether the two labellings agree on it: both put the two events in the same cluster, or both put them in different clusters. The fraction of pairs on which they agree is the Rand index. It does not care what the clusters are called, so cluster 3 in one labelling may be cluster 0 in the other.
• Two unrelated labellings already agree on many pairs by chance. The adjusted index subtracts the agreement expected by chance and rescales, so that $\text{ARI} = 1$ means identical groupings and $\text{ARI} \approx 0$ means no better than chance.

We repeat the split several times and average, because one split is itself random.

                          map  random
k-means, k = 4           0.87    0.81

k-means, k = 8           0.78    0.59
HDBSCAN, min size 15     0.81    0.04

Each half has half the events, so HDBSCAN's sizes are halved too (min_cluster_size 15 and min_samples 5 instead of 30 and 10). On the map, every method is reproducible. On random events the two methods part ways. HDBSCAN's random clusters are not reproducible: each half has its own chance fluctuations. k-means with four clusters is reproducible even on random events, because four tiles of a uniform square land in nearly the same place every time. So stability alone does not prove structure: it shows that the result is not an accident of the sample, and it must be read with a null test. An unstable setting, though, is a warning on any data.

Test 3: a feature the method never saw

The strongest test asks whether the clusters are coherent in something they were not built from. Each event carries a strike, held back in A.1. If a cluster is one fault, its events should share a strike more closely than a random group of events of the same size. Strike is an angle that repeats every 180°, so we double each angle before averaging it as a direction: the length $R$ of the mean of the unit vectors at twice the strikes is 1 when all strikes agree and near 0 when they are random.

HDBSCAN cluster 0: 319 events, strike coherence R = 0.87; random groups of that size: 0.20 (largest of 500: 0.33)
HDBSCAN cluster 1: 160 events, strike coherence R = 0.87; random groups of that size: 0.21 (largest of 500: 0.37)
HDBSCAN cluster 2: 136 events, strike coherence R = 0.24; random groups of that size: 0.20 (largest of 500: 0.37)
HDBSCAN cluster 3: 208 events, strike coherence R = 0.89; random groups of that size: 0.21 (largest of 500: 0.36)

Three clusters share a strike far more closely than any random group: the two faults and the patch. The fourth, the sparse cluster, does not; its strikes were planted at random. The held-out feature agrees with Test 1 that this cluster is the weakest, and adds something position cannot: even if it is a real group, nothing says its events broke one structure. A group can be real without being one fault. Part B makes the same test on the Geysers catalogue, holding back the focal mechanism while clustering on position.

A.9 · When time joins position: what does "close" mean?

So far every event was a point on a map, and the distance between two events was a distance in kilometres. Nobody had to decide what "close" meant. A catalogue also gives each event an origin time, and much of what seismologists want to find lives in time: an aftershock sequence, the burst of events that follows a larger earthquake and decays with time, and a swarm, a burst of many events of similar size with no dominant earthquake, which Vidale and Shearer (2006) attribute mainly to pore-fluid pressure changes or slow aseismic slip rather than to one event triggering the next.

We build a synthetic catalogue with four kinds of events, so the truth is known:

• background: events at random places and random times, the steady seismicity of a region;
• a steady patch: events crowded in one place but spread evenly over the whole time span, like a persistently active volume around an injection well. Clustered in space, not in time;
• an aftershock sequence: one larger earthquake, then events close to it whose rate decays with time after it, following the Omori–Utsu law (reviewed by Utsu, Ogata & Matsu'ura 1995): the rate falls as $1/(t + c)^p$, with $t$ the time since the mainshock and $c$ and $p$ constants, $p$ usually near 1;
• a swarm: many small events over a few days, of similar magnitude, spreading outward from a point the way a pressure change diffuses, so that the distance reached grows as the square root of time.

Every event also has a magnitude, drawn from the Gutenberg–Richter relation $\log_{10} N = a - b m$ (Aki 1965 gives its standard estimate of $b$): the number of events above magnitude $m$ falls by a factor $10^{b}$ for each unit of $m$, with $b$ near 1, so small events are far more common than large ones.

background   :  1500 events, largest magnitude 3.8
steady patch :   300 events, largest magnitude 3.8
aftershocks  :   201 events, largest magnitude 4.5
swarm        :   200 events, largest magnitude 2.5

On the map the steady patch, the aftershocks and the swarm all look like dense blobs. Against time they differ: the steady patch is a horizontal band, active throughout; the aftershocks and the swarm are short vertical streaks. A method that should find sequences, and not steady activity, has to use time.

The unit problem

Put time next to position as a third feature, $(x, y, t)$, and cluster with DBSCAN. The distance between two events is now

$$d = \sqrt{\Delta x^2 + \Delta y^2 + (v\,\Delta t)^2},$$

and $v$, a speed in kilometres per day, says how many kilometres one day is worth. Nothing in the data fixes $v$: kilometres and days are different units, and any conversion is a choice. We cluster the same catalogue with the same eps and min_samples at three values of $v$.

v = 0.0001 km/day: 4 clusters; share in a cluster: steady patch 100%, aftershocks 100%, swarm 100%
v = 0.1 km/day: 2 clusters; share in a cluster: steady patch 0%, aftershocks 83%, swarm 100%
v = 10 km/day: 3 clusters; share in a cluster: steady patch 0%, aftershocks 33%, swarm 12%

• With $v$ tiny, time is ignored. The clusters are the dense places on the map, and the steady patch is one of them, although it is not a sequence at all. (The fourth cluster is a chance group of background events.)
• With one day worth a tenth of a kilometre, the steady patch disappears (its events are close in space but far apart in time), and the aftershocks and the swarm are found.
• With $v$ large, position is ignored and the clusters are moments of high activity. Only the earliest, most intense part of the aftershock sequence survives, and most of the swarm is lost.

So the answer depends on a unit conversion that the data do not determine. Seismologists have used exactly this distance: Frohlich and Davis (1990), as described by van Stiphout, Zhuang and Marsan (2012), linked earthquakes by $\sqrt{r^2 + C^2 \Delta t^2}$ with one day counted as one kilometre. This is not special to time. Any two features in different units need a conversion; a common rule is to divide each feature by a scale that means something: a location by its location error (Session 2, Location & relocation), a focal mechanism by its uncertainty as an angle (Session 3, Focal mechanism & moment tensor). Then one unit of every feature is "one uncertainty", and a distance of one is "as far apart as the errors allow". Time has no measurement error worth the name, so it has no such scale, and its conversion has to come from the physics.

Plotting has the same problem. A depth section drawn with a vertically stretched axis makes compact clusters look like columns; a map drawn with unequal east and north scales makes round clusters look like streaks. Draw distances at a common scale, as every map in this notebook does.

A.10 · A distance from physics: the nearest-neighbour method

Seismology has its own answer to the unit problem, the nearest-neighbour distance of Baiesi and Paczuski (2004), developed into a clustering method by Zaliapin et al. (2008) and Zaliapin and Ben-Zion (2013). Instead of choosing how many kilometres a day is worth, it asks a question with no free units: if seismicity were random, how many events would we expect to see between these two?

Take an earlier event $i$ of magnitude $m_i$ and a later event $j$, a time $\Delta t$ apart and a distance $r$ apart. In a random catalogue the number of events expected in a time window $\Delta t$ and within distance $r$ is proportional to $\Delta t$ times the area within $r$. That area grows as $r^{d_f}$, where $d_f$ is the fractal dimension of the epicentres: 2 if they fill the map evenly, less if they crowd onto faults. And by Gutenberg–Richter, events at least as large as $m_i$ are rarer by $10^{-b\,m_i}$, so a large event is a more surprising neighbour. Their product is the nearest-neighbour proximity

$$\eta_{ij} = \Delta t \; r^{d_f} \; 10^{-b\, m_i} \qquad (\Delta t > 0).$$

Small $\eta$ means that event $j$ came so soon and so close after $i$ that chance is an unlikely explanation. Changing the units of time or distance multiplies every $\eta$ by the same constant, so it changes no comparison: the conversion problem of A.9 is gone. What remains are $d_f$ and $b$, which are properties of the catalogue that can be measured. This catalogue was built with $b = 1$ and epicentres that fill the map, so $d_f = 2$; on a real catalogue both are estimated from the data, as Part B does.

The method then has three steps.

1. Link each event to its nearest neighbour. For every event $j$, find the earlier event $i$ with the smallest $\eta_{ij}$: its parent. Every event except the first now has one link to an earlier event, so the catalogue becomes a set of trees.
2. Keep the strong links. A link with small $\eta$ is a likely cause-and-effect pair; a link with large $\eta$ is just the nearest of many random neighbours. Cut every link with $\eta$ above a threshold $\eta_0$.
3. Read the families. The trees left after cutting are families: a first event and every event linked to it, directly or through others.

2201 events linked; log10 eta ranges from -10.4 to 1.0

Choosing the threshold with a null catalogue

The split of $\eta$ into a time part $T$ and a space part $R$ (each takes half the magnitude factor, so $\eta = T R$) gives the method's best-known picture: a 2-D histogram of $\log T$ against $\log R$. To see what chance alone produces, we build a null catalogue as in A.8, keeping everything except the structure under test. Here the structure is clustering in time, so the null keeps every position and magnitude and gives each event a random time. The steady patch stays dense in space in the null; only the bursts in time are destroyed.

threshold log10 eta0 = -4.15 (chance links at most 5% of links kept)
strong links: 16.8% of events in the catalogue, 0.7% in the null

The rule asks: among the links we keep, how many would a random catalogue also produce? Lowering $\eta_0$ keeps fewer links and fewer chance ones; we take the largest $\eta_0$ at which chance links are at most 5% of the links kept, the same rule Part B uses on the Geysers catalogue. This is not how Zaliapin and Ben-Zion chose theirs: they set $\eta_0 = 10^{-5}$ by inspecting the two populations in the histogram below, and cite a formal alternative that fits a two-component Gaussian mixture, the model of A.3, to the two populations. The null rule has the advantage of saying how many of the kept links chance alone would produce. The figure shows the two histograms and the threshold, a line of constant $\eta = TR$.

The catalogue has two populations. Most events sit in a large cloud at large $T$ and $R$, the same cloud the null has: these are background and steady events whose nearest neighbour is just the nearest of many random ones. A second, smaller population lies down and to the left, short times and short distances, and the null has nothing there: these are the clustered events. The threshold separates them.

This picture says only whether an event is clustered, not what kind of cluster it belongs to: the aftershocks and the swarm both lie in the lower-left population, and their links overlap there. Telling them apart needs the next step, the shape of each family.

Families: bursts and swarms

Now keep the strong links and follow them back to each family's first event. Which events did the method put together, and do they match what we planted?

family of 184 events: swarm 184
family of 169 events: aftershocks 169
family of  11 events: swarm 11
family of   4 events: steady patch 4
family of   2 events: background 2
background   :  0.1% of events have a strong link to a parent
steady patch :  2.3% of events have a strong link to a parent
aftershocks  : 83.6% of events have a strong link to a parent
swarm        : 96.5% of events have a strong link to a parent

The largest families are the aftershock sequence and the swarm, each with almost no events of another kind, and almost no background or steady-patch event has a strong link. The swarm is split into one large family and smaller pieces: a chain breaks wherever one link is weaker than the threshold, and every piece after the break becomes a family of its own.

The two big families differ in shape. In the aftershock family most events link straight to the mainshock: the large magnitude makes $10^{-b m_i}$ tiny, so the mainshock is the nearest neighbour of nearly every later event. The tree is a star. In the swarm no event is much larger than the others, so each event's nearest neighbour is a recent event close by, and the tree is a chain, each event linked to the one before. Zaliapin and Ben-Zion measure the shape by the average leaf depth $\langle d \rangle$: a leaf is an event with no children in its family, its depth is the number of links from it back to the family's first event, and $\langle d \rangle$ is the average depth over the family's leaves. A star has $\langle d \rangle$ near 1; a chain has a large one.

aftershock family: 169 events, average leaf depth 1.13, deepest event 3 links, largest magnitude 4.5
swarm family: 184 events, average leaf depth 13.43, deepest event 25 links, largest magnitude 2.5

Each event is drawn at its depth, the number of links back to the family's first event, and each grey line is one strong link from a parent to its child. In the aftershock family nearly every event sits at depth 1: its parent is the mainshock. In the swarm family the depth keeps growing with time, one link after another: each event's parent is a recent event close by.

On southern California, Zaliapin and Ben-Zion found that families of small-to-medium earthquakes fall into two main types by this measure: burst-like families, with small average leaf depth, and swarm-like families, with larger depth that grows with the family's size, with intermediate forms between and no sharp boundary. They suggest that swarm-like families are associated with regions of relatively high temperature or fluid content, which is one reason the measure matters for a geothermal field. Part B computes it for the Geysers.

A family with only a few events is a short tree whatever caused it, so leaf depth classifies only families with enough events; and the classification depends on the threshold, which decides where chains break.

The programs the field runs

For the general methods, scikit-learn (Pedregosa et al. 2011) implements k-means, Gaussian mixtures, agglomerative clustering, DBSCAN and HDBSCAN, all used above; the hdbscan package (McInnes, Healy & Astels 2017) adds HDBSCAN's cluster tree, soft membership and prediction for new events; scipy builds and draws dendrograms.

For earthquakes, grouping events into sequences is usually called declustering: separating a catalogue into background events and the events triggered by them, often to estimate the background rate for seismic hazard. Four approaches dominate (reviewed by van Stiphout, Zhuang & Marsan 2012):

• window methods (Gardner & Knopoff 1974): every event within a time window and a distance window, both set by the mainshock's magnitude, is an aftershock. Simple, and the window sizes are chosen by hand and vary between studies;
• link methods (Reasenberg 1985): an event joins a cluster when it falls inside an interaction zone in space and time around the cluster's events, with the zone's size set by a physical model; clusters grow by chains of links, as in single linkage (A.4);
• nearest-neighbour methods (Baiesi & Paczuski 2004; Zaliapin et al. 2008; Zaliapin & Ben-Zion 2013), as in A.10;
• stochastic declustering (Zhuang, Ogata & Vere-Jones 2002), which fits a model of triggering, the ETAS model (epidemic-type aftershock sequence), and gives each event a probability of being background, much as the mixture's responsibilities in A.3 gave each event a probability of belonging to each cluster.

For structure in space rather than in time, Ouillon & Sornette (2011) fit fault planes to hypocentres with a Gaussian mixture, and use cross-validation to choose the number of planes.

A different clustering problem is phase association: grouping the P and S arrivals picked at many stations into the earthquakes that produced them, before any event has a location. GaMMA (Zhu et al. 2022) treats it as a Gaussian mixture whose clusters are earthquakes, each with a location, origin time and magnitude, fitted by the same expectation–maximisation as A.3.

Takeaways

Machine learning

• A clustering is fixed by three choices — what counts as close, what counts as a group, and at what scale — so without labels no answer is correct in itself.
• k-means, the Gaussian mixture and Ward linkage put every event in a cluster; DBSCAN and HDBSCAN can leave background events out, which is why density methods suit seismicity.
• k-means and EM stop at a local optimum that depends on the start, so run many starts and keep the best objective.
• A label-free diagnostic chooses well only when the data match what it assumes; a result that stays the same across its parameter, as HDBSCAN's does, is stronger evidence.
• Every method finds clusters in random data: a cluster is credible when it beats every cluster in a null catalogue, is reproduced by random halves of the data, and is coherent in a feature the method never saw.

Seismology

• Position alone cannot tell steady activity from a sequence; time is needed, and a distance mixing days and kilometres depends on a conversion the data do not fix.
• The nearest-neighbour proximity $\eta = \Delta t\, r^{d_f}\, 10^{-b m}$ measures how unlikely two events are to be neighbours by chance, which removes that conversion.
• An aftershock sequence forms a star-shaped family around its largest event and a swarm a chain of similar events, and the average leaf depth measures the difference.

References

• Aki, K. (1965). Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits. Bull. Earthq. Res. Inst. 43, 237–239. 10.15083/0000033631
• Baiesi, M., & Paczuski, M. (2004). Scale-free networks of earthquakes and aftershocks. Phys. Rev. E 69, 066106. 10.1103/PhysRevE.69.066106
• Campello, R. J. G. B., Moulavi, D., & Sander, J. (2013). Density-based clustering based on hierarchical density estimates. PAKDD 2013, LNCS 7819, 160–172. 10.1007/978-3-642-37456-2_14
• Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. B 39, 1–22. 10.1111/j.2517-6161.1977.tb01600.x
• Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. (1996). A density-based algorithm for discovering clusters in large spatial databases with noise. KDD-96, 226–231. PDF
• Frohlich, C., & Davis, S. D. (1990). Single-link cluster analysis as a method to evaluate spatial and temporal properties of earthquake catalogues. Geophys. J. Int. 100, 19–32 (cited as described by van Stiphout et al. 2012).
• Gardner, J. K., & Knopoff, L. (1974). Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian? BSSA 64, 1363–1367. 10.1785/BSSA0640051363
• Hubert, L., & Arabie, P. (1985). Comparing partitions. J. Classification 2, 193–218. 10.1007/BF01908075
• Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Trans. Inf. Theory 28, 129–137. 10.1109/TIT.1982.1056489
• MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. Proc. Fifth Berkeley Symp. Math. Stat. Prob. 1, 281–297.
• McInnes, L., Healy, J., & Astels, S. (2017). hdbscan: Hierarchical density based clustering. J. Open Source Softw. 2, 205. 10.21105/joss.00205
• Ouillon, G., & Sornette, D. (2011). Segmentation of fault networks determined from spatial clustering of earthquakes. J. Geophys. Res. 116, B02306. 10.1029/2010JB007752
• Pedregosa, F., et al. (2011). Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830. JMLR
• Reasenberg, P. (1985). Second-order moment of central California seismicity, 1969–1982. J. Geophys. Res. 90, 5479–5495. 10.1029/JB090iB07p05479
• Rousseeuw, P. J. (1987). Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 20, 53–65. 10.1016/0377-0427(87)90125-7
• Schwarz, G. (1978). Estimating the dimension of a model. Ann. Stat. 6, 461–464. 10.1214/aos/1176344136
• Utsu, T., Ogata, Y., & Matsu'ura, R. S. (1995). The centenary of the Omori formula for a decay law of aftershock activity. J. Phys. Earth 43, 1–33. 10.4294/jpe1952.43.1
• van Stiphout, T., Zhuang, J., & Marsan, D. (2012). Seismicity declustering. Community Online Resource for Statistical Seismicity Analysis. 10.5078/corssa-52382934
• Vidale, J. E., & Shearer, P. M. (2006). A survey of 71 earthquake bursts across southern California: exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers. J. Geophys. Res. 111(B5). 10.1029/2005JB004034
• Ward, J. H. (1963). Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 58, 236–244. 10.1080/01621459.1963.10500845
• Zaliapin, I., & Ben-Zion, Y. (2013). Earthquake clusters in southern California I: Identification and stability. J. Geophys. Res. Solid Earth 118, 2847–2864. 10.1002/jgrb.50179
• Zaliapin, I., & Ben-Zion, Y. (2013). Earthquake clusters in southern California II: Classification and relation to physical properties of the crust. J. Geophys. Res. Solid Earth 118, 2865–2877. 10.1002/jgrb.50178
• Zaliapin, I., Gabrielov, A., Keilis-Borok, V., & Wong, H. (2008). Clustering analysis of seismicity and aftershock identification. Phys. Rev. Lett. 101, 018501. 10.1103/PhysRevLett.101.018501
• Zhuang, J., Ogata, Y., & Vere-Jones, D. (2002). Stochastic declustering of space-time earthquake occurrences. J. Am. Stat. Assoc. 97, 369–380. 10.1198/016214502760046925
• Zhu, W., McBrearty, I. W., Mousavi, S. M., Ellsworth, W. L., & Beroza, G. C. (2022). Earthquake phase association using a Bayesian Gaussian mixture model. J. Geophys. Res. Solid Earth 127, e2021JB023249. 10.1029/2021JB023249