28. Hypothesis tests and paired data

Language: Python (Jupyter notebook) - R version

Two surfaces that had no R or Python binding until this release: the twelve Numerics HypothesisTests statics, reached directly through hypothesis_test() or, for a flood-frequency record, through the DataFrame facades analysis_data_hypothesis_test()/analysis_data_statistics(); and the RMC.BestFit “Paired Data” subsystem – an x-y curve container with interpolation, area, three curve-simplification algorithms, and an uncertain twin whose y-coordinate is a whole distribution per point. This page works through both, on a real flood record and a real reservoir curve.

What you’ll learn

What four hypothesis tests answer about a flood-frequency record – normality, independence, and homogeneity – and how to read a p-value from each.
The DataFrame path through analysis_data(): the same tests over a censored-data-aware container, the use_log10 switch, and the two-sample split_index split.
correlation()’s matrix form over more than two series at once.
A reservoir stage-storage curve interpolated in linear and log space, and its trapezoidal area.
How the three curve-simplification algorithms differ on the same curve, including one that drops the curve’s last point without warning.
Sampling an uncertain curve at its mean versus its median.

Setup

import numpy as np

import corehydropy as ch

A flood-frequency record

Sixty-nine years of annual peak discharge (m³/s) for the Harricana River, Quebec (Bobee & Ashkar 1991, Multivariate Flood Frequency Analysis, Table 1.2):

harricana = [
    122, 244, 214, 173, 229, 156, 212, 263, 146, 183, 161, 205, 135, 331, 225,
    174, 98.8, 149, 238, 262, 132, 235, 216, 240, 230, 192, 195, 172, 173, 172,
    153, 142, 317, 161, 201, 204, 194, 164, 183, 161, 167, 179, 185, 117, 192,
    337, 125, 166, 99.1, 202, 230, 158, 262, 154, 164, 182, 164, 183, 171, 250,
    184, 205, 237, 177, 239, 187, 180, 173, 174,
]
len(harricana)

Four questions a flood-frequency analysis usually asks of a record before fitting anything to it: is it normally distributed, is it independent from year to year, and is it homogeneous (no jump, no trend)? hypothesis_test() answers each with a p-value.

{
    "jarque_bera": ch.hypothesis_test(harricana, method="jarque_bera")["p_value"],
    "wald_wolfowitz": ch.hypothesis_test(harricana, method="wald_wolfowitz")["p_value"],
    "ljung_box": ch.hypothesis_test(harricana, method="ljung_box")["p_value"],
    "mann_kendall": ch.hypothesis_test(harricana, method="mann_kendall")["p_value"],
}

{'jarque_bera': 0.0010565797084869377,
 'wald_wolfowitz': 0.24359142152435065,
 'ljung_box': 0.8204107252215709,
 'mann_kendall': 0.775729412824826}

Jarque-Bera tests normality, and its p-value here, 0.001, rejects it firmly – unsurprising for a right-skewed peak-flow record (skewness 0.86 below), which is exactly why flood-frequency work fits a skewed family instead of a Normal one. Wald-Wolfowitz (a runs test) and Ljung-Box (serial correlation up to a default lag) both test independence between consecutive years; neither rejects it. Mann-Kendall tests for a monotonic trend – homogeneity in the “no trend” sense – and its own large p-value agrees: no trend detected over these 69 years.

Mann-Whitney tests the other kind of homogeneity, a jump: are the first and second halves of the record drawn from the same distribution? It is the one two-sample test above, so it needs both samples split out by hand when called directly (analysis_data_hypothesis_test() below does the split for you, from a split_index):

first_half = harricana[0:34]
second_half = harricana[34:69]
ch.hypothesis_test(first_half, second_half, method="mann_whitney")

{'p_value': 0.687614027389208}

The `DataFrame` path

analysis_data() wraps a record the way an actual flood-frequency analysis would – with room for censored observations the plain list above has none of. Its own hypothesis-test and summary-statistics facades read the same twelve HypothesisTests statics, plus a use_log10 switch and a split-index-based two-sample split:

d = ch.analysis_data(harricana)
print(d)

print(ch.analysis_data_hypothesis_test(d, "mann_kendall"))
ch.analysis_data_hypothesis_test(d, "mann_whitney", split_index=34)

<AnalysisData 69 exact>
{'mann_kendall': 0.775729412824826}

{'mann_whitney': 0.687614027389208}

The split_index split is on the record’s INDEX (0-based, so split_index=34 puts the observations with index below 34 – the first 34 – in the first sample and the rest in the second) – matching the manual harricana[0:34]/harricana[34:69] split above exactly.

use_log10=True reruns the same test against log-space values. For a rank-based test like Mann-Kendall this changes nothing – log10 is monotonic, so it cannot change which of two years’ peaks is larger, and Mann-Kendall depends only on those pairwise orderings:

{
    "real_space": ch.analysis_data_hypothesis_test(d, "mann_kendall")["mann_kendall"],
    "log_space": ch.analysis_data_hypothesis_test(d, "mann_kendall", use_log10=True)["mann_kendall"],
}

{'real_space': 0.775729412824826, 'log_space': 0.775729412824826}

Jarque-Bera is not rank-based, so use_log10 changes its answer outright – this is the log-normal version of the normality question asked above, and it comes back with a very different verdict:

{
    "real_space": ch.analysis_data_hypothesis_test(d, "jarque_bera")["jarque_bera"],
    "log_space": ch.analysis_data_hypothesis_test(d, "jarque_bera", use_log10=True)["jarque_bera"],
}

{'real_space': 0.0010565797084869377, 'log_space': 0.5756256964701716}

Log10-transformed, the record no longer rejects normality – the log-Normal (or log-Pearson) family a flood-frequency analysis usually reaches for is a reasonable fit exactly because this record’s SHAPE, not just its scale, is closer to Normal after the transform.

analysis_data_statistics() reports the twenty summary-statistics keys DataFrame computes – moments, percentiles, both in real space and log space:

s = ch.analysis_data_statistics(d)
{k: s["value"][k] for k in ["Record Length", "Mean", "Std Dev", "Skewness", "Kurtosis", "50%"]}

{'Record Length': 69.0,
 'Mean': 191.3173913043478,
 'Std Dev': 47.96161135411182,
 'Skewness': 0.8605451107460523,
 'Kurtosis': 4.34348681301944,
 '50%': 183.0}

A correlation matrix

correlation() takes a 2D array instead of two vectors and returns every pairwise coefficient at once. Ten years of peak flow at three tributary gauges:

gauge_a = [122, 244, 214, 173, 229, 156, 212, 263, 146, 183]
gauge_b = [98, 210, 187, 145, 201, 130, 178, 240, 120, 155]
gauge_c = [310, 180, 405, 220, 350, 190, 300, 410, 175, 260]

m = ch.correlation(np.column_stack([gauge_a, gauge_b, gauge_c]))
m

array([[1.        , 0.99696338, 0.46537499],
       [0.99696338, 1.        , 0.50432343],
       [0.46537499, 0.50432343, 1.        ]])

Gauges A and B track closely (a shared upstream driver, high peak years at one are high peak years at the other); gauge C answers to something else – its correlation with either is much weaker.

A reservoir stage-storage curve

The rest of this page is the Paired Data subsystem: OrderedPairedData, an x-y curve that keeps itself sorted and validated, with interpolation, trapezoidal area, and three ways to reduce its point count. The curve below is upstream’s own 15-point reservoir fixture – storage in acre-feet against pool elevation in feet:

storage = [230408, 288010, 345611, 403213, 460815, 518417, 576019, 633612,
           691223, 748825, 806427, 864029, 921631, 1036834, 1152038]
elevation = [1519.7, 1520.5, 1520.9, 1521.7, 1523.5, 1525.9, 1528.4, 1530.9,
             1533.2, 1534.7, 1535.9, 1538, 1541.3, 1547.7, 1552.7]

curve_interpolate() reads the elevation at a given storage, linearly by default:

target = 500000
lin = ch.curve_interpolate(storage, elevation, xout=[target])[0]
log_space = ch.curve_interpolate(storage, elevation, xout=[target], x_transform="logarithmic")[0]
{"linear": lin, "log_x": log_space}

{'linear': np.float64(1525.1326516440402),
 'log_x': np.float64(1525.1629478762325)}

The two differ by a small fraction of a foot here – the curve is nearly linear over this stretch, so which space you interpolate in barely matters; it would matter much more over a stretch where storage grows convexly with elevation.

curve_area() is the trapezoidal-rule area between the curve and the x-axis:

ch.curve_area(storage, elevation)

1413175623.4

Three ways to simplify the same curve

curve_simplify() reduces a curve’s point count three ways – Ramer-Douglas-Peucker, Visvalingam- Whyatt, and Lang – each trading a tolerance (or a target count) for fewer points:

rdp = ch.curve_simplify(storage, elevation, method="rdp", tolerance=0.5)
vis = ch.curve_simplify(storage, elevation, method="visvalingam", num_to_keep=7)
lang = ch.curve_simplify(storage, elevation, method="lang", tolerance=0.5, look_ahead=2)

for name, arr in [("rdp", rdp), ("visvalingam", vis), ("lang", lang)]:
    print(f"{name:12s} retained_points={arr.shape[0]:2d}  last_x_matches={arr[-1, 0] == storage[-1]}")

rdp          retained_points= 7  last_x_matches=True
visvalingam  retained_points= 7  last_x_matches=True
lang         retained_points= 8  last_x_matches=False

Ramer-Douglas-Peucker and Visvalingam-Whyatt both always keep the curve’s first and last point. Lang does not, and this is not merely a theoretical concern raised elsewhere on this site – at this tolerance and look-ahead it drops the reservoir curve’s own last point right here, the same defect docs/upstream-csharp-issues.md documents against the sin curve upstream’s own Test_LangSimplify uses. last_x_matches is False for Lang and True for the other two.

Sampling an uncertain curve

UncertainOrderedPairedData is the same idea with a distribution at every y instead of a number – a stage-storage curve, say, whose elevation at each storage level is only known within a range. uncertain_curve_sample() collapses it to a plain curve at a chosen quantile, or at the mean:

storage_pts = storage[0:5]
elev_dists = [ch.Distribution("Triangular", [e - 0.6, e, e + 1.2]) for e in elevation[0:5]]

median_curve = ch.uncertain_curve_sample(storage_pts, elev_dists, probability=0.5)
mean_curve = ch.uncertain_curve_sample(storage_pts, elev_dists)

np.column_stack([storage_pts, median_curve[:, 1], mean_curve[:, 1]])

array([[230408.        ,   1519.86076952,   1519.9       ],
       [288010.        ,   1520.66076952,   1520.7       ],
       [345611.        ,   1521.06076952,   1521.1       ],
       [403213.        ,   1521.86076952,   1521.9       ],
       [460815.        ,   1523.66076952,   1523.7       ]])

The Triangular distributions here are not symmetric (each is skewed further above its mode than below), so the mean and median curves differ by a consistent small offset – exactly what you would expect from a right-skewed uncertainty band at every point, and exactly what collapsing to a single number by mean instead of by median would hide.

Key takeaways

hypothesis_test() reaches all twelve HypothesisTests statics directly; analysis_data_*() reaches nine of them through a DataFrame, with the censored-data and use_log10/split_index machinery a raw vector does not have.
correlation() without a second argument treats its input as a table and returns the whole pairwise matrix in one call.
curve_interpolate()/curve_area() work in linear or log-transformed space on the same curve.
Ramer-Douglas-Peucker and Visvalingam-Whyatt always keep a curve’s endpoints; Lang does not – an upstream defect this port reproduces on purpose rather than silently fixing.
uncertain_curve_sample() collapses a distribution-valued curve to a plain one at any quantile, including the mean and the median, which need not agree.

Reproduction check

All values below are deterministic given the fixed inputs above; checked at 1e-15 relative tolerance except the Mann-Kendall p-value, which is additionally pinned against the upstream C# test literal (ExactDataHypothesisTests.Test_MannKendall, fixtures/data/data_frame_facades.json) at that test’s own tolerance.

def near(x, literal, tol=1e-15):
    return abs(x / literal - 1) < tol

mk = ch.analysis_data_hypothesis_test(d, "mann_kendall")["mann_kendall"]
mw = ch.analysis_data_hypothesis_test(d, "mann_whitney", split_index=34)["mann_whitney"]
jb_log = ch.analysis_data_hypothesis_test(d, "jarque_bera", use_log10=True)["jarque_bera"]
mw_manual = ch.hypothesis_test(first_half, second_half, method="mann_whitney")["p_value"]

checks = [
    # Pinned against the upstream C# oracle.
    abs(mk - 0.7757) < 1e-3,

    # This page's own deterministic numbers, checked bit-for-bit against themselves.
    near(mk, ch.hypothesis_test(harricana, method="mann_kendall")["p_value"]),
    near(mw, mw_manual),

    # The correlation matrix, pinned to this page's own executed values (not merely checked for
    # symmetry, which is guaranteed by construction and proves nothing about the actual numbers).
    near(m[0, 1], 0.996963377946883),
    near(m[0, 2], 0.4653749875464116),
    near(m[1, 2], 0.5043234334007348),

    near(ch.curve_area(storage, elevation), 1413175623, tol=1e-9),

    # The two interpolations, pinned to this page's own executed values (not merely recomputed
    # and compared against themselves, which is vacuous).
    near(lin, 1525.1326516440402),
    near(log_space, 1525.1629478762325),

    rdp.shape[0] == 7,
    vis.shape[0] == 7,
    lang.shape[0] == 8,
    bool(rdp[-1, 0] == storage[-1]),
    bool(vis[-1, 0] == storage[-1]),
    not bool(lang[-1, 0] == storage[-1]),  # the LangSimplify drop, reproduced on purpose
    jb_log > 0.5,                          # log-space does not reject normality
    near(mean_curve[0, 1], 1519.8999999999999),
    near(median_curve[0, 1], 1519.8607695154587),
    near(mean_curve[0, 1] - median_curve[0, 1], 0.039230484541121768),
]
assert all(checks)
print("All reproduction checks passed.")

All reproduction checks passed.

What you’ll learn

Setup

A flood-frequency record

The DataFrame path

A correlation matrix

A reservoir stage-storage curve

Three ways to simplify the same curve

Sampling an uncertain curve

Key takeaways

Reproduction check

The `DataFrame` path