import numpy as np
from corehydropy import (
ml_gaussian_mixture,
ml_glm,
ml_jenks_breaks,
ml_kmeans,
ml_knn,
ml_naive_bayes,
ml_random_forest,
)29. Machine learning
The Numerics MachineLearning layer, which had no R or Python binding until this release: three unsupervised methods (k-means, Gaussian mixture models, Jenks natural breaks) and five supervised ones (decision trees, random forests, k-nearest neighbors, Gaussian naive Bayes, and generalized linear models). This page works through all eight on one regional flood-frequency problem – the kind of regionalization and prediction-in-ungauged-basins task these methods actually get used for in hydrology.
What you’ll learn
- How to group catchments into homogeneous regions with
ml_kmeans()andml_gaussian_mixture(), and what the mixture model tells you that k-means cannot. - How
ml_jenks_breaks()cuts a continuous variable into hazard classes, and why it will happily give you a class of one. - How to fit a regional flood-frequency regression with
ml_glm(), and a Poisson count model with the same verb and a different link. - How
ml_random_forest()compares against that regression out of sample – including a case where the parametric model wins, and why. - How
ml_naive_bayes()andml_knn()classify catchments above and below a flood threshold. - Why every fit on this page is reproducible, and what the seeds are actually doing.
Setup
A regional dataset
Thirty-six gauged catchments, each with a drainage area (km2), mean annual precipitation (mm), mean basin slope (m/m), and mean annual flood (m3/s). The values are illustrative rather than observed – they are generated from a power-law regional relation with a deterministic wobble, so the page’s numbers are stable and the R and Python versions define byte-identical inputs – but the shape and the ranges are those of a real regional study.
area = np.array([
14.1, 11.4, 10.0, 13.1, 13.7, 10.5,
10.7, 13.9, 12.9, 9.9, 11.7, 14.2,
165.0, 133.6, 116.7, 152.5, 160.1, 122.4,
124.4, 161.7, 150.0, 115.8, 136.2, 165.2,
1060.6, 858.6, 750.0, 980.0, 1029.4, 786.6,
799.9, 1039.2, 964.3, 744.3, 875.9, 1062.0,
])
precip = np.array([
606, 634, 682, 590, 665, 657,
593, 686, 626, 613, 690, 600,
928, 971, 1044, 904, 1018, 1006,
908, 1051, 958, 938, 1057, 919,
1373, 1437, 1545, 1337, 1506, 1488,
1343, 1554, 1418, 1388, 1564, 1359,
], dtype=float)
slope = np.array([
0.0911, 0.0821, 0.0677, 0.0589, 0.0622, 0.0753,
0.0881, 0.0910, 0.0818, 0.0675, 0.0588, 0.0624,
0.0389, 0.0350, 0.0289, 0.0251, 0.0266, 0.0321,
0.0376, 0.0388, 0.0349, 0.0288, 0.0251, 0.0266,
0.0134, 0.0120, 0.0099, 0.0086, 0.0091, 0.0110,
0.0129, 0.0133, 0.0120, 0.0099, 0.0086, 0.0092,
])
flood = np.array([
2.18, 1.90, 1.79, 1.69, 2.13, 1.80,
1.73, 2.47, 2.09, 1.50, 1.98, 1.83,
15.67, 15.17, 12.87, 13.48, 15.24, 14.37,
12.41, 19.79, 14.91, 12.05, 14.15, 14.68,
68.03, 66.88, 55.60, 59.71, 65.53, 63.93,
53.11, 88.39, 63.52, 54.04, 60.00, 66.11,
])
len(area)36Basin characteristics span orders of magnitude, so every method on this page works in log10 space. That is not a detail: k-means and k-nearest neighbors measure Euclidean distance, so an untransformed area column would swamp a slope column outright.
X = np.column_stack([np.log10(area), np.log10(precip), np.log10(slope)])
COLS = ["log_area", "log_precip", "log_slope"]
np.round(X[:3], 4)array([[ 1.1492, 2.7825, -1.0405],
[ 1.0569, 2.8021, -1.0857],
[ 1. , 2.8338, -1.1694]])Regionalization: k-means
Regionalization is the task of grouping catchments that behave alike, so a flood estimate can be transferred from gauged to ungauged basins within a group. ml_kmeans() does it by minimizing within-cluster distance.
km = ml_kmeans(X, k=3, seed=20260827)
np.bincount(km["labels"])array([12, 12, 12])labels are 0-based in both Python and R, matching the library’s own indexing. Here the three clusters recover the three catchment size classes exactly, twelve each:
np.round(km["means"], 4)array([[ 2.9565, 3.1585, -1.9713],
[ 1.0817, 2.8033, -1.137 ],
[ 2.1484, 2.9884, -1.5069]])The centroids are in log10 space, so the third cluster’s log_area of 2.96 is a catchment of about 900 km2.
km["iterations"]2Regionalization: a Gaussian mixture
ml_gaussian_mixture() fits the same kind of grouping but carries a full covariance matrix per component rather than only a centre, so it can represent clusters that are elongated or correlated rather than spherical. That matters here because area, precipitation and slope co-vary strongly.
gmm = ml_gaussian_mixture(X, k=3, seed=20260827)
np.round(gmm["weights"], 6)array([0.333333, 0.333333, 0.333333])Equal mixing weights, and the same three groups. What the mixture adds is the covariance:
np.round(gmm["sigmas"][2], 5)array([[ 0.00328, -0.00019, 0.0002 ],
[-0.00019, 0.0006 , -0.00016],
[ 0.0002 , -0.00016, 0.00496]])The off-diagonal terms are the useful part – within the largest catchments, log_area and log_slope move together (steeper basins in this set are smaller), which k-means has no way to express.
log_likelihood follows the library’s own definition, which omits the multivariate-normal normalizing constant. It is short of the true mixture log-likelihood by X.shape[0] * X.shape[1] / 2 * log(2 * pi). That constant cancels when comparing two fits of the same data, which is what the library uses it for, but do not feed this value straight to AIC.
Hazard classes: Jenks natural breaks
ml_jenks_breaks() cuts a continuous variable into classes that minimize within-class variance – the standard way to build a legend for a hazard map.
jn = ml_jenks_breaks(flood, n_clusters=4)
np.round(jn["breaks"], 2)array([ 2.47, 19.79, 68.03, 88.39])# columns: start_index, end_index, count, min, max, sum, average, variance
np.round(jn["clusters"][:, [2, 3, 4]], 2)array([[12. , 1.5 , 2.47],
[12. , 12.05, 19.79],
[11. , 53.11, 68.03],
[ 1. , 88.39, 88.39]])Note the fourth class: a single catchment. Jenks minimizes within-class variance and nothing else, so a genuine outlier (the 88.4 m3/s catchment) gets a class to itself rather than being absorbed. That is correct behaviour and usually a signal to look at that gauge, not to change n_clusters.
round(jn["gvf"], 6)0.987557The goodness-of-variance-fit measure runs from 0 to 1; 0.988 says four classes explain nearly all the variance in the flood column.
Regional regression: ml_glm()
The classic regional flood-frequency model is a power law in the basin characteristics, which is a linear model in log space. Split the catchments into a training set and a held-out test set of nine:
train = np.array([i - 1 for i in list(range(1, 10)) + list(range(13, 22)) + list(range(25, 34))])
test = np.array([i for i in range(len(flood)) if i not in set(train.tolist())])
fit = ml_glm(X[train], np.log10(flood[train]))
np.round(fit["coefficients"], 5)array([ 0.01087, 0.92989, -0.17801, 0.2044 ])Read as a power law, the fitted exponents are 0.930 on area, -0.178 on precipitation and 0.204 on slope. Their standard errors say how much to trust each:
np.round(fit["standard_errors"], 5)array([0.8775 , 0.06683, 0.34598, 0.10998])np.round(fit["p_values"], 5)array([0.99012, 0. , 0.60691, 0.0631 ])Only the area exponent is well determined – unsurprising, since area spans two orders of magnitude here while precipitation spans less than one.
A Poisson model with the same verb
ml_glm()’s link argument selects the response family as well as the transform, so a count response needs only a different link. Here the count is the number of days per year each catchment exceeds its local flood threshold:
exceedances = np.array([
1, 2, 3, 1, 2, 2, 1, 3, 2, 1, 3, 1,
4, 5, 7, 4, 6, 6, 4, 7, 5, 4, 7, 4,
9, 11, 14, 8, 12, 12, 8, 14, 10, 9, 14, 9,
], dtype=float)
Xp = np.column_stack([np.log10(precip), np.log10(area)])
pois = ml_glm(Xp, exceedances, link="log")
np.round(pois["coefficients"], 5)array([ 0.0132 , -0.12084, 0.92862])round(pois["aic"], 3)142.878The fitted counts track the observed ones closely:
design = np.column_stack([np.ones(len(Xp)), Xp])
fitted_counts = np.exp(design @ pois["coefficients"])
round(float(np.corrcoef(fitted_counts, exceedances)[0, 1]), 6)0.911463The other three links, "logit", "probit" and "complementary_log_log", take a 0/1 response and fit the corresponding binomial model.
Random forest, and when it loses
ml_random_forest() grows many trees on bootstrap resamples and reports the spread of their predictions as an interval, which is what makes it attractive when the functional form is unknown.
rf = ml_random_forest(X[train], np.log10(flood[train]), newdata=X[test],
seed=20260827, number_of_trees=200)
# Columns are lower, median, upper, mean.
rf.shape(9, 4)np.round(rf[:3], 4)array([[0.2529, 0.2553, 0.3201, 0.2699],
[0.2279, 0.3284, 1.1951, 0.5425],
[0.2776, 0.3385, 1.1951, 0.6306]])Compare both models out of sample:
def rmse(a, b):
return float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) ** 2)))
lm_pred = np.column_stack([np.ones(len(test)), X[test]]) @ fit["coefficients"]
rf_median = rf[:, 1]
{"glm": rmse(lm_pred, np.log10(flood[test])),
"random_forest": rmse(rf_median, np.log10(flood[test]))}{'glm': 0.04068522131045717, 'random_forest': 0.06857045960049413}The regression wins, by a factor of about 1.7 in RMSE. That is the right answer, and it is worth being explicit about why, because the usual demonstration runs the other way:
- The data really is a power law, so a model that assumes one is using information the forest has to rediscover from 27 training points.
- A tree predicts by averaging the training responses in a leaf, so it cannot extrapolate. Every prediction is inside the range of the training data, which is exactly wrong for a regional relation you want to apply to an ungauged basin outside the gauged range.
Reach for the forest when you suspect interactions or thresholds you cannot write down, not as a default. The same trade-off is why a bare ml_decision_tree() is rarely the end of an analysis: at its defaults a regression tree recurses until every leaf holds one training observation, so it memorizes the training set.
Classification: above or below a threshold
Recast the problem as a yes-or-no question – will this catchment’s mean annual flood exceed 20 m3/s? – and two classifiers answer it.
big = (flood > 20).astype(float)
nb = ml_naive_bayes(X[train], big[train], newdata=X[test])
knn_pred = ml_knn(X[train], big[train], newdata=X[test], k=3, regression=False)
np.column_stack([big[test], nb["prediction"], knn_pred])array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[1., 1., 1.],
[1., 1., 1.],
[1., 1., 1.]])Both classify all nine held-out catchments correctly. Naive Bayes also reports the per-class feature distributions it fitted, which is often more informative than the predictions:
np.round(nb["means"], 4)array([[ 1.6183, 2.8965, -1.3003],
[ 2.9598, 3.1591, -1.9499]])np.round(nb["priors"], 4)array([0.6667, 0.3333])classes follows the order the labels first appear in the training response, not sorted order, and every row of means, standard_deviations and priors is indexed to match.
ml_knn() can also report which training catchments a prediction came from, which is exactly the “donor catchment” idea regional hydrology uses directly:
neighbours = ml_knn(X[train], big[train], newdata=X[test], k=3, what="neighbors")
neighbours[:3]array([[2, 5, 1],
[4, 3, 2],
[3, 4, 8]])These are 0-based indices into the training rows, like the cluster labels.
Reproducibility
Every method on this page that draws a random number takes a seed, and a seeded fit is reproducible: repeated calls agree bit for bit, and so do the R and Python versions of this page, because the entire computation happens in the shared compiled core with the same bit-exact Mersenne Twister.
again = ml_random_forest(X[train], np.log10(flood[train]), newdata=X[test],
seed=20260827, number_of_trees=200)
bool(np.array_equal(again, rf))TrueLeave seed at its default and the generator is seeded from the clock instead, which is the right choice when you want an honest sense of run-to-run variability and the wrong one when you want a result you can publish.
Reproduction check
Every value below is deterministic given the fixed inputs and seeds above, checked at 1e-15 relative tolerance against this page’s own executed values. The k-means and naive-Bayes paths are additionally pinned against the upstream C# oracles in fixtures/ml/machine_learning.json, which the dotnet gate replays against the real library; those fixtures use the iris measurements the C# tests use, so what is checked here is that this page’s own numbers do not drift.
def near(x, literal, tol=1e-15):
return abs(x / literal - 1) < tol
# Clustering recovers the three size classes exactly.
assert np.bincount(km["labels"]).tolist() == [12, 12, 12]
assert km["iterations"] == 2
assert near(km["means"][0, 0], 2.9565201679743947)
assert near(km["means"][0, 2], -1.9713001801472083)
# The mixture puts equal weight on the three components and its covariances are symmetric.
assert np.allclose(gmm["weights"], 1 / 3, atol=1e-12)
assert near(gmm["log_likelihood"], 238.87583471023325)
assert np.allclose(gmm["sigmas"][2], gmm["sigmas"][2].T, atol=1e-12)
# Jenks, including the single-member top class.
assert near(jn["breaks"][3], 88.390000000000001)
assert jn["clusters"][:, 2].tolist() == [12.0, 12.0, 11.0, 1.0]
assert near(jn["gvf"], 0.98755725105630288)
# The regional regression.
assert near(fit["coefficients"][1], 0.92989165483884118)
assert near(fit["standard_errors"][1], 0.066834707598246085)
assert near(fit["sigma"], 0.04726688233226619)
assert fit["n"] == 27 and fit["df"] == 23
# The Poisson model.
assert near(pois["coefficients"][2], 0.92861744959233672)
assert near(pois["aic"], 142.87786000842539)
assert near(float(np.corrcoef(fitted_counts, exceedances)[0, 1]), 0.91146273024268387)
# The forest is seeded-reproducible, its interval is ordered, and the regression beats it here.
assert (rf[:, 0] <= rf[:, 1]).all() and (rf[:, 1] <= rf[:, 2]).all()
assert near(rmse(lm_pred, np.log10(flood[test])), 0.040685221310457167)
assert near(rmse(rf_median, np.log10(flood[test])), 0.068570459600494127)
assert rmse(lm_pred, np.log10(flood[test])) < rmse(rf_median, np.log10(flood[test]))
# Both classifiers get all nine held-out catchments right.
assert nb["prediction"].tolist() == big[test].tolist()
assert np.asarray(knn_pred).tolist() == big[test].tolist()
assert nb["classes"].tolist() == [0.0, 1.0]
assert neighbours.min() >= 0 and neighbours.max() < 27
assert neighbours.shape == (9, 3)
print("All reproduction checks passed.")All reproduction checks passed.