Model selection

Perhaps we do not know the number of segments in a rating or that a power law is better than another type of model. In these situations, we want to compare the performance of one model against another, then select whichever models gives the best fit for the data. More specifically, we want to know the generalization performance of each model, which refers to their performance on new (out-of-sample) data.

One way to estimate generalization performance is cross validation, but it can be costly. In cross validation, a portion of the data is held out. If you don’t have a lot of data to begin with, holding out a portion can substantially degrade model fit.

An alternative approach is to use information criteria, which are measures that use in-sample data to estimate out-of-sample performance.

In this tutorial, we demonstrate how to select the number of segments in the Green River rating. We fit the rating with 1, 2, 3, and 4 segments, then use information criteria to select the number of segments that gives the best fit.

# Load the data
from ratingcurve import data

%load_ext autoreload
%autoreload 2
%xmode minimal
# suppress warnings and errors

import arviz as az
import pymc as pm

from ratingcurve.ratings import PowerLawRating

df = data.load('green channel')

Exception reporting mode: Minimal

Fit the data

Fit the data to ratings with 1 to 4 segments.

%%capture
# Output supressed, this will print "Finished" after running each of the four models

segments = [1, 2, 3, 4]
traces = []
for segment in segments:
    powerrating = PowerLawRating(segments=segment,
                                 prior={'distribution': 'uniform'})

    trace = powerrating.fit(q=df['q'],
                            h=df['stage'],
                            q_sigma=df['q_sigma'],
                            n=100_000)

    traces.append(pm.compute_log_likelihood(trace, model=powerrating.model))  # Add arg to compute log likelihood

Convergence achieved at 43900

Interrupted at 43,899 [43%]: Average Loss = 220.95

Sampling: [a, b, hs_, model_q, sigma]

Sampling: [model_q]

Finished [100%]: Average Loss = -52.049

Sampling: [a, b, hs_, model_q, sigma]

Sampling: [model_q]

Finished [100%]: Average Loss = -48.728

Sampling: [a, b, hs_, model_q, sigma]

Sampling: [model_q]

Finished [100%]: Average Loss = -44.67

Sampling: [a, b, hs_, model_q, sigma]

Sampling: [model_q]

Now use arviz.compare to format the output.

# this model will generate warnings about the LOO
import warnings

warnings.filterwarnings('ignore')

compare_dict = {f'{i} segment': traces[i - 1] for i in segments}
az.compare(compare_dict)

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight
2 segment	0	0.0	0.0	NaN		3 k̂ > 0.70	5.8	73.0	3.2	0.74
4 segment	1	-10.0	9.6	0.90	N < 100	9 k̂ > 0.70	18.2	60.0	9.7	0.00
3 segment	2	-10.0	13.0	0.85	N < 100	4 k̂ > 0.70	18.4	60.0	13.0	0.26
1 segment	3	-19.0	4.3	1.00	N < 100		3.6	54.0	3.3	0.00

As expected, the 2-segment model ranked highest.

Residual analysis

In practice, it can be helpful to plot to rating error (the deviations between the rating fit and the discharge observations). Here is a demonstration of how.

segments = 2

powerrating = PowerLawRating(segments=2,
                             prior={'distribution': 'uniform'})

trace = powerrating.fit(q=df['q'],
                        h=df['stage'],
                        q_sigma=df['q_sigma'])

Finished [100%]: Average Loss = -52.185

Sampling: [a, b, hs_, model_q, sigma]

Sampling: [model_q]

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, figsize=(7, 7), sharey=True)

powerrating.plot(ax[0])
powerrating.plot_residuals(ax[1])


plt.subplots_adjust(wspace=0.1)
ax[1].set_ylabel('')

Sampling: [model_q]

Sampling: [model_q]

Text(0, 0.5, '')

%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,xarray

Last updated: Tue, 01 Sep 2026

Python implementation: CPython
Python version       : 3.12.14
IPython version      : 9.17.1

pytensor: 3.3.0
xarray  : 2026.7.0

arviz      : 1.3.0
matplotlib : 3.11.1
pymc       : 6.3.1
ratingcurve: 0.1.dev1+ga5a0fe6ba

Watermark: 2.6.0