Flow velocity of the River Thames

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)

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?

### YOUR CODE HERE

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.

### YOUR CODE HERE

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.

### YOUR CODE HERE

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.

### YOUR CODE HERE

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.

### YOUR CODE HERE

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.

### YOUR CODE HERE

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()

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

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.

### YOUR CODE HERE

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')

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.