Flow velocity of the River Thames

Use this notebook to carry out the analysis, inspect the results, and prepare a concise set of handwritten notes. In particular, make sure you complete the goodness-of-fit assessment in Part 4, the dependence analysis in Part 6, and the two raw contingency tables in Part 7 before laptops close.

What is the propagated uncertainty? How fast does the River Thames flow at Kingston?

You have paired observations of water depth \(h\) [m] and flow rate (discharge) \(q\) [m\(^3\)/s] from a measuring station on the River Thames at Kingston, UK. The observations were collected every 15 minutes over one month.

Discharge \(q\) is related to mean flow velocity \(u\) by

\[q = uS,\]

where \(S\) is the cross-sectional area. For this simplified analysis, approximate the cross-section as rectangular with an effective river width of

\[w=95\ \mathrm{m},\]

so that

\[S \approx wh, \qquad u=\frac{q}{wh}=\frac{q}{95h}.\]

Your goals are to:

1. explore and model the marginal distributions of \(h\) and \(q\);

2. assess the fitted distributions graphically;

3. propagate the fitted distributions to \(u\), initially assuming \(h\) and \(q\) are independent;

4. compare the dependence in the observations with that in the independent simulation;

5. prepare two raw contingency tables from the course questionnaire for the paper-only part of the GA.

The Thames data originate from the UK Department for Environment, Food & Rural Affairs’ Data Services platform.

import os
from urllib.request import urlretrieve

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from IPython.display import display

plt.rcParams.update({'font.size': 14})

Part 1. Explore the data

def findfile(fname):
    """Download a MUDE source file only if it is not already in the repository."""
    if not os.path.isfile(fname):
        print(f'Downloading {fname}...')
        urlretrieve(
            'https://github.com/TUDelft-MUDE/source-files/raw/main/file/GA1.4/' + fname,
            fname
        )

findfile('dataset_thames.csv')
h, q = np.genfromtxt('dataset_thames.csv', delimiter=',', unpack=True, skip_header=True)

Downloading dataset_thames.csv...

Task 1 — explore the observations:

Plot the observations of \(h\) and \(q\) and summarize each variable numerically. At minimum, determine its mean, standard deviation, range, and skewness.

Compare the relative variability of the two variables using the coefficient of variation,

\[CV = \frac{s}{\bar{x}}.\]

What features of the two empirical distributions would matter when selecting a parametric model later?

fig, ax = plt.subplots(2, 1, figsize=(10, 7), layout='constrained')
ax[0].plot(h, 'k')
ax[0].set_ylabel('Water depth, $h$ [m]')
ax[0].set_xlabel('Observation')
ax[0].grid()
ax[1].plot(q, 'k')
ax[1].set_ylabel('Flow rate, $q$ [m$^3$/s]')
ax[1].set_xlabel('Observation')
ax[1].grid()

for name, x in [('h', h), ('q', q)]:
    desc = stats.describe(x)
    s = np.std(x, ddof=1)
    cv = s / np.mean(x)
    print(f'{name}:')
    print(f'  n = {len(x)}')
    print(f'  mean = {np.mean(x):.3f}')
    print(f'  std = {s:.3f}')
    print(f'  min/max = {np.min(x):.3f}, {np.max(x):.3f}')
    print(f'  skewness = {desc.skewness:.3f}')
    print(f'  CV = {cv:.3f}')

h:
  n = 2296
  mean = 3.402
  std = 0.023
  min/max = 3.338, 3.483
  skewness = 0.692
  CV = 0.007
q:
  n = 2296
  mean = 8.657
  std = 2.592
  min/max = 0.784, 19.580
  skewness = 1.145
  CV = 0.299

Solution

The water depth varies only slightly during the measurement period. Its mean is about \(3.402\) m and its standard deviation is about \(0.023\) m, giving \(CV(h)\approx0.007\). In contrast, the mean discharge is about \(8.657\) m\(^3\)/s with a standard deviation of about \(2.592\) m\(^3\)/s, giving \(CV(q)\approx0.299\). Thus, discharge shows much greater variability relative to its typical magnitude.

Both variables are positively skewed, but the effect is much stronger for \(q\) (skewness \(\approx1.15\)) than for \(h\) (skewness \(\approx0.69\)). The distribution of \(h\) is therefore concentrated in a narrow range with a modest right tail, whereas \(q\) has a much more pronounced right tail.

Part 2. Empirical distribution functions

The function below is given — you used the same idea in the workshop.

def ecdf(var):
    x = np.sort(var)
    n = x.size
    y = np.arange(1, n + 1) / (n + 1)
    return [y, x]

Task 2 — visualize the empirical distributions:

For both \(h\) and \(q\), make a density-normalized histogram and an empirical CDF. Use the plots together with your statistics from Task 1 to characterize the shape of each distribution.

fig, axes = plt.subplots(2, 2, figsize=(12, 9), layout='constrained')

axes[0, 0].hist(h, bins=20, density=True, edgecolor='k', linewidth=0.2)
axes[0, 0].set_xlabel('Water depth, $h$ [m]')
axes[0, 0].set_ylabel('Density')
axes[0, 0].set_title('$h$: histogram')
axes[0, 0].grid()

axes[0, 1].hist(q, bins=20, density=True, edgecolor='k', linewidth=0.2)
axes[0, 1].set_xlabel('Flow rate, $q$ [m$^3$/s]')
axes[0, 1].set_ylabel('Density')
axes[0, 1].set_title('$q$: histogram')
axes[0, 1].grid()

F_h, h_sorted = ecdf(h)
F_q, q_sorted = ecdf(q)
axes[1, 0].step(h_sorted, F_h, where='post')
axes[1, 0].set_xlabel('Water depth, $h$ [m]')
axes[1, 0].set_ylabel(r'$P(H \leq h)$')
axes[1, 0].set_title('$h$: empirical CDF')
axes[1, 0].grid()

axes[1, 1].step(q_sorted, F_q, where='post')
axes[1, 1].set_xlabel('Flow rate, $q$ [m$^3$/s]')
axes[1, 1].set_ylabel(r'$P(Q \leq q)$')
axes[1, 1].set_title('$q$: empirical CDF')
axes[1, 1].grid()

Part 3. Fitting a distribution

Task 3 — choose and fit parametric models:

Based on Tasks 1 and 2, choose one candidate distribution for each variable:

• for \(h\): Uniform or Lognormal;

• for \(q\): Gaussian or Gumbel.

Briefly justify each choice, then fit the selected distributions using scipy.stats (MLE). Store the fitted parameters as params_h and params_q.

Your justification here.

params_h = stats.lognorm.fit(h)
params_q = stats.gumbel_r.fit(q)

print('Lognormal params for h (shape, loc, scale):', params_h)
print('Gumbel params for q (loc, scale):', params_q)

Lognormal params for h (shape, loc, scale): (np.float64(0.1921240552490838), 3.2835837648621853, np.float64(0.11666407834912991))
Gumbel params for q (loc, scale): (np.float64(7.498689888072776), 2.1212576118774233)

Solution

For \(h\), the Lognormal distribution is the more plausible of the two candidates. The observed depths are strongly concentrated around approximately \(3.4\) m rather than being spread uniformly between their minimum and maximum values. The observations also have a modest positive skewness, which a Lognormal distribution can represent whereas a Uniform distribution cannot.

For \(q\), the Gumbel distribution is the more plausible candidate. The discharge observations have a pronounced positive skewness and a clear right tail. A Gaussian distribution is symmetric, whereas the right-skewed Gumbel distribution can reproduce right tails.

Part 4. Assessing goodness of fit

A fitted distribution is a model, not a guarantee that the observations actually follow it. Use a graphical goodness-of-fit check to compare each fitted distribution with the empirical distribution. Inspect the agreement across the observed range, paying particular attention to regions where the curves diverge. No formal hypothesis test is required here.

Task 4 — graphical goodness of fit:

For both \(h\) and \(q\), compare the empirical exceedance probability with the fitted model on a logarithmic probability axis.

For each variable, decide whether the fitted distribution represents the observations reasonably well. While the plots are visible, write a short note by hand recording your overall judgement and the main evidence from the plot that supports it. Make a qualitative sketch of these plots for later.

fig, axes = plt.subplots(1, 2, figsize=(14, 5), layout='constrained')

F_h, h_sorted = ecdf(h)
axes[0].step(h_sorted, 1 - F_h, color='k', label='$h$ (observed)')
axes[0].plot(h_sorted, stats.lognorm.sf(h_sorted, *params_h), '--', color='grey', label='Lognormal')
axes[0].set_xlabel('Water depth, $h$ [m]')
axes[0].set_ylabel('$P(H > h)$')
axes[0].set_title('$h$')
axes[0].set_yscale('log')
axes[0].legend()
axes[0].grid()

F_q, q_sorted = ecdf(q)
axes[1].step(q_sorted, 1 - F_q, color='k', label='$q$ (observed)')
axes[1].plot(q_sorted, stats.gumbel_r.sf(q_sorted, *params_q), '--', color='grey', label='Gumbel')
axes[1].set_xlabel('Flow rate, $q$ [m$^3$/s]')
axes[1].set_ylabel('$P(Q > q)$')
axes[1].set_title('$q$')
axes[1].set_yscale('log')
axes[1].legend()
axes[1].grid()

Solution

For \(h\), the fitted Lognormal distribution follows the empirical exceedance curve fairly closely through most of the observed range. The agreement is not exact, however: the empirical curve is less smooth, and under- or over-estimates the exceedance probabilities in different parts of the graph. In the right tail, in particular, we observe significant deviations from the fitted distribution; this is in part because a log plot amplifies mismatches for small \(\log y\) values.

For \(q\), the fitted Gumbel distribution also captures the overall trend of the empirical exceedance curve quite well. As for \(h\), differences become visible toward some of the largest observed discharges.

Overall, both distributions provide acceptable approximations to the bulk of their respective datasets, but neither model reproduces every part of the empirical distribution exactly.

Part 5. Propagating the uncertainty

Using the fitted distributions, propagate the uncertainty from \(h\) and \(q\) to flow velocity \(u\), assuming that \(h\) and \(q\) are independent. You will examine that assumption in Section 6.

Recall

\[u=\frac{q}{95h}.\]

Task 5 — Monte Carlo propagation:

1. Draw 10,000 independent random samples from the fitted distributions for \(h\) and \(q\).

2. Compute the corresponding simulated velocities \(u_{sim}=q_{sim}/(95h_{sim})\).

3. Compute the observed velocities \(u_{obs}=q/(95h)\) from the paired measurements.

4. Compare the observed and simulated velocity distributions using density-normalized histograms and exceedance plots.

5. Report the mean and standard deviation of both velocity distributions.

Use variable names rs_h, rs_q, rs_u, and u so they can be reused below.

rng_mc = np.random.default_rng(42)
rs_h = stats.lognorm.rvs(*params_h, size=10000, random_state=rng_mc)
rs_q = stats.gumbel_r.rvs(*params_q, size=10000, random_state=rng_mc)

w = 95.0
rs_u = rs_q / (w * rs_h)
u = q / (w * h)

fig, axes = plt.subplots(1, 2, figsize=(13, 5), layout='constrained')
axes[0].hist(rs_u, bins=25, density=True, alpha=0.6, edgecolor='k', linewidth=0.2, label='Simulated')
axes[0].hist(u, bins=25, density=True, alpha=0.6, edgecolor='k', linewidth=0.2, label='Observed')
axes[0].set_xlabel('Flow velocity, $u$ [m/s]')
axes[0].set_ylabel('Density')
axes[0].set_title('Velocity distribution')
axes[0].legend()
axes[0].grid()

F_rs_u, rs_u_sorted = ecdf(rs_u)
F_u, u_sorted = ecdf(u)
axes[1].step(rs_u_sorted, 1 - F_rs_u, label='Simulated')
axes[1].step(u_sorted, 1 - F_u, label='Observed')
axes[1].set_xlabel('Flow velocity, $u$ [m/s]')
axes[1].set_ylabel('$P(U > u)$')
axes[1].set_title('Velocity exceedance')
axes[1].set_yscale('log')
axes[1].legend()
axes[1].grid()

print('Observed u:  mean=%.4f, std=%.4f' % (np.mean(u), np.std(u, ddof=1)))
print('Simulated u: mean=%.4f, std=%.4f' % (np.mean(rs_u), np.std(rs_u, ddof=1)))

Observed u:  mean=0.0267, std=0.0078
Simulated u: mean=0.0271, std=0.0085

Part 6. Scatter plot and correlation

Task 6 — compare dependence:

Make a scatter plot of \(h\) against \(q\) for:

• the paired observations \((h_i,q_i)\); and

• the independently generated simulations \((h_{sim},q_{sim})\).

Compute the Pearson correlation coefficient for each set and store the results as rho_obs and rho_sim.

Copy both correlation coefficients into your handwritten reference notes.

fig, ax = plt.subplots(figsize=(7, 7))
ax.scatter(rs_h, rs_q, s=18, alpha=0.25, label='Independent simulations')
ax.scatter(h, q, s=18, marker='x', label='Paired observations')
ax.set_xlabel('Water depth, $h$ [m]')
ax.set_ylabel('Flow rate, $q$ [m$^3$/s]')
ax.legend()
ax.grid()

rho_obs = np.corrcoef(h, q)[0, 1]
rho_sim = np.corrcoef(rs_h, rs_q)[0, 1]
print('Correlation, observations: rho=%.3f' % rho_obs)
print('Correlation, simulations:  rho=%.3f' % rho_sim)

Correlation, observations: rho=0.753
Correlation, simulations:  rho=-0.002

Part 7. Contingency tables: prepare two raw tables

This section uses a different dataset: this year’s course questionnaire (dataset_class_survey_contingency.csv). Very rare response categories were removed or grouped before release.

For the setup phase, your job is only to construct and copy raw count tables.

findfile('dataset_class_survey_contingency.csv')
survey = pd.read_csv('dataset_class_survey_contingency.csv')
survey.head()

Downloading dataset_class_survey_contingency.csv...

	programme	gender	world_region	cats_or_dogs	rock_paper_scissors_round_1	rock_paper_scissors_round_2	rock_paper_scissors_round_3	vacation_style	modelling_uncertainty_or_data	shares_apartment
0	Civil Engineering	Female	Europe	Cats	Paper	Scissors	Rock	Sport holiday	*Incoherent Screaming*	I live alone
1	Earth, Climate, and Technology	Male	Europe	Dogs	Rock	Paper	Rock	Mountain hike	*Incoherent Screaming*	I live with others
2	Earth, Climate, and Technology	Female	Europe	Both?	Scissors	Rock	Paper	Nature park	*Incoherent Screaming*	I live with others
3	Civil Engineering	Male	Europe	Dogs	Scissors	Scissors	Scissors	City trip	Modelling	I live with others
4	Civil Engineering	Male	Europe	Dogs	Rock	Rock	Rock	Beach holiday	Modelling	I live with others

Task 7.1 — worked example:

The cell below shows the mechanics of constructing a contingency table. Read it and make sure everyone in the group understands what the rows, columns, cells, and totals represent.

Do not use this worked example as one of your own two tables below.

# WORKED EXAMPLE -- illustrates the mechanic only
example_counts = pd.crosstab(
    survey['programme'],
    survey['cats_or_dogs'],
    margins=True
)
example_counts

cats_or_dogs	Both?	Cats	Dogs	Neither! >:(	All
programme
Civil Engineering	16	23	30	7	76
Earth, Climate, and Technology	4	9	6	2	21
Environmental Engineering	9	4	3	1	17
All	29	36	39	10	114

Task 7.2 — choose and prepare two tables:

Choose two different pairs of categorical variables from the questionnaire that you find worth comparing. Construct a raw contingency table for each pair, including row and column totals.

Make sure the two tables use different variable pairs. Copy both complete raw count tables, including labels and totals, accurately into your handwritten reference notes. You will decide how to use one of these tables only after the written questions are handed out.

# Example solution choices -- your own pairs may be different
pair1 = ('programme', 'vacation_style')
pair2 = ('world_region', 'cats_or_dogs')

counts1 = pd.crosstab(survey[pair1[0]], survey[pair1[1]], margins=True)
counts2 = pd.crosstab(survey[pair2[0]], survey[pair2[1]], margins=True)

display(counts1)
display(counts2)

vacation_style	Antarctica cruise	Beach holiday	City trip	Mountain hike	Nature park	Road trip	Sport holiday	All
programme
Civil Engineering	4	17	12	12	7	19	5	76
Earth, Climate, and Technology	1	3	5	5	3	4	0	21
Environmental Engineering	2	2	3	2	5	2	1	17
All	7	22	20	19	15	25	6	114

cats_or_dogs	Both?	Cats	Dogs	Neither! >:(	All
world_region
Americas	4	3	8	0	15
Asia	6	4	1	1	12
Europe	19	29	30	9	87
All	29	36	39	10	114

Part 8. Prepare your handwritten reference notes — do this before 11:30!

Once laptops close, you may use only the handwritten notes your group prepared during the setup.

The numerical cell below collects the values you will need. Your goodness-of-fit judgements and contingency tables must come from your own work above.

print('='*50)
print('HANDWRITTEN REFERENCE NOTES -- COPY THESE RESULTS BEFORE 11:30')
print('='*50)
print()

print('--- Distribution fitting results ---')
print('q empirical moments: mean=%.3f, std=%.3f, skewness=%.3f' % (np.mean(q), np.std(q, ddof=1), stats.describe(q).skewness))
print('q MLE fitted parameters: '+"%.3f "*len(params_q) % params_q)
print('h empirical moments: mean=%.3f, std=%.3f, skewness=%.3f' % (np.mean(h), np.std(h, ddof=1), stats.describe(h).skewness))
print('h MLE fitted parameters: '+"%.3f "*len(params_h) % params_h)
print()

print('--- Dependence results ---')
print('rho_obs = %.3f' % rho_obs)
print('rho_sim = %.3f' % rho_sim)
print()

print('Also copy by hand:')
print('  1. Your short GOF judgement + supporting graphical evidence (qualitatively) for h')
print('  2. Your short GOF judgement + supporting graphical evidence (qualitatively) for q')
print('  3. A brief justification of your distribution choices for h and q')
print('  4. BOTH complete raw contingency tables from Section 7')

==================================================
HANDWRITTEN REFERENCE NOTES -- COPY THESE RESULTS BEFORE 11:30
==================================================

--- Distribution fitting results ---
q empirical moments: mean=8.657, std=2.592, skewness=1.145
q MLE fitted parameters: 7.499 2.121 
h empirical moments: mean=3.402, std=0.023, skewness=0.692
h MLE fitted parameters: 0.192 3.284 0.117 

--- Dependence results ---
rho_obs = 0.753
rho_sim = -0.002

Also copy by hand:
  1. Your short GOF judgement + supporting graphical evidence (qualitatively) for h
  2. Your short GOF judgement + supporting graphical evidence (qualitatively) for q
  3. A brief justification of your distribution choices for h and q
  4. BOTH complete raw contingency tables from Section 7

By Max Ramgraber, Patricia Mares Nasarre and Robert Lanzafame, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.