Fitting probability distributions

Fitting probability distributions#

Fitting probability distributions to data is an important task in any discipline of science and engineering, as these distributions can be used to derive quantitative statements about risks and frequencies of uncertain properties. In this workshop, you will first work with a dataset of your choosing to estimate a distribution, and evaluate the quality of your fit. Then, you will build and analyze a contingency table that describes the relationship between two variables.

The goal of this project is:

  1. Choose a reasonable distribution function for your chosen dataset, analyzing the statistics of the observations.

  2. Fit the chosen distribution by moments.

  3. Assess the fit computing probabilities analytically.

  4. Assess the fit using goodness-of-fit techniques and computer code.

  5. Build a contingency table from two variables and reason about their dependence.

The workshop is divided into 4 parts: 1) data analysis, 2) pen-and-paper practice, 3) programming, and 4) contingency tables.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import pandas as pd
from scipy import stats

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

Please choose one of the following datasets:

  1. Observations of the compressive strength of concrete. The compressive strength of concrete is key for the safety of infrastructures and buildings. However, a lot of boundary conditions influence the final resistance of the concrete, such as the cement content, the environmental temperature, or the age of the concrete.

  2. ERA5 predictions of the hourly temperature at 2m height during the summer months (June, July, August) from 2005 to 2025. ERA5 predictions are re-analysis data, which are observation-corrected weather model predictions.

import os
from urllib.request import urlretrieve

def findfile(fname):
    if not os.path.isfile(fname):
        print(f'Downloading {fname}...')
        urlretrieve('https://github.com/TUDelft-MUDE/source-files/raw/main/file/'+fname, fname)

findfile('dataset_concrete.csv')
findfile('dataset_temperature.csv')

# Please choose one of the datasets below
viable_datasets = ['concrete', 'temperature']
dataset = 'concrete'  # Choose one dataset from the list above

assert dataset in viable_datasets, 'Dataset must be in {}.'.format(str(viable_datasets))

data = np.genfromtxt('dataset_{}.csv'.format(dataset), skip_header=1)

if dataset == 'concrete':
    label = 'Concrete compressive strength [MPa]'
    number_bins = 10
elif dataset == 'temperature':
    label = 'Summer air temperature in Delft [°C]'
    number_bins = 20

Now, let us clean and plot the data.

data = data[~np.isnan(data)]

plt.figure(figsize=(10, 6))
gs = GridSpec(nrows=1, ncols=2, width_ratios=[1, 0.2])

plt.subplot(gs[0, 0])
plt.plot(data, 'ok')
plt.xlabel('observation number'); plt.ylabel(label)
ylims = plt.gca().get_ylim()

plt.subplot(gs[0, 1])
plt.hist(data, orientation='horizontal', color='lightblue', rwidth=0.9, bins=number_bins)
plt.xlabel('frequency'); plt.gca().set_ylim(ylims)

In the figure above, you can see all the observations of your chosen dataset. There is no clear pattern in the observations. Let’s see how the statistics look!

df_describe = pd.DataFrame(data)
df_describe.describe()

Task 1.1

Using ONLY the statistics calculated above:

  • Choose an appropriate distribution to model the data between the following: (1) Gumbel, (2) Uniform, and (3) Gaussian.

  • Justify your choice.

Part 2: Use pen and paper!#

Once you have selected the appropriate distribution, fit it by moments manually and check the fit by computing some probabilities analytically. Remember that you have all the information you need in the textbook. Do not use any computer code for this section — you have to do it with pen and paper. You can use the notebook as a calculator.

Task 2.1:

Fit the selected distribution by moments.

We can now check the fit by computing manually some probabilities from the fitted distribution and comparing them with the empirical ones.

Task 2.2:

Check the fit of the distribution:

  • Use the min, 25%, 50%, 75%, and max values. What are the non-exceedance probabilities (from the empirical distribution) that correspond to those values?

  • Compute the values of the random variable corresponding to those probabilities using the fitted distribution.

  • Compare the obtained values with the empirical ones and assess the fit.

Tip: to compute the predicted minimum and maximum, use the non-exceedance probability of the first and last rank sample in an empirical CDF (\(1/(n+1)\) and \(n/(n+1)\)).

Your answer here — fill in the table below.

Minimum

P25%

P50%

P75%

Maximum

Non-exceedance prob.

Empirical

Predicted

Part 3: Let’s do it with Python!#

Now, let’s assess performance with further goodness-of-fit techniques and see whether they are consistent with the analysis done by hand. Note that you have the pseudo-code for the empirical CDF in the reader.

Task 3.1:

Prepare a function to compute the empirical cumulative distribution function.

def ecdf(observations):
    """Return [non_exceedance_probabilities, sorted_values]."""
    ### YOUR CODE HERE
    return [non_exceedance_probabilities, sorted_values]

Task 3.2:

Transform the parameters of the selected distribution you fitted by moments to loc-scale-shape.

Hint: the distributions are listed in the online textbook, but check that the book’s formulation matches the Python package’s (see the Scipy.stats documentation).

Your answer here.

Task 3.3:

Assess the goodness of fit of the distribution you fitted by moments by:

  • Visually comparing the empirical and fitted PDF.

  • Using the exceedance plot in log-scale.

  • Using the QQplot.

  • Interpreting them — do you reach a conclusion similar to Part 2?

# Visually compare the PDFs
### YOUR CODE HERE

# Visually compare the exceedance probabilities
### YOUR CODE HERE

# Visually compare the QQ plots
### YOUR CODE HERE

Your answer here.

End of Part 3. So far, everything you’ve fitted and assessed has been about a single variable at a time. In Part 4, we ask a different kind of question: how do two variables relate to each other? We’ll answer it using tables built directly from data, connecting to the discrete joint/marginal/conditional probability material from this week’s lecture.

Case: height and shoe size in this year’s MUDE class#

In Parts 1-3, you fitted a continuous parametric distribution to a single variable. Here, we look at two continuous variables at once — height and shoe size, from this year’s course questionnaire — and ask a new kind of question: not “what distribution fits each one on its own?”, but “how are they related?”. We’ll answer that using the discrete joint / marginal / conditional probability tools from Monday’s lecture, applied to height and shoe size after grouping each into a few categories.

The goal of this part is:

  1. Look at the raw relationship between two continuous variables.

  2. Turn each variable into a small number of categories.

  3. Build a contingency table, and compute joint, marginal, and conditional probabilities from it.

  4. Use those probabilities to informally check whether the two variables look independent.

A note on the data: dataset_class_survey_2026.csv only includes responses from students who agreed to let us reuse their answers in future years (112 of 127 respondents). No names, emails, or timestamps are included.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

findfile('dataset_class_survey_2026.csv')

survey = pd.read_csv('dataset_class_survey_2026.csv')
survey[['height_cm', 'shoe_size']].describe()

Task 4.1:

Make a scatter plot of shoe_size (y-axis) against height_cm (x-axis). Do the two variables look related?

### YOUR CODE HERE

Task 4.2 — clean the data first:

The questionnaire asked respondents to report their shoe size using the EU sizing system. Inspect the smallest values of shoe_size. A few responses appear inconsistent with the requested scale and may have been entered using a different convention, such as foot length in centimetres.

Identify those observations and decide how to handle them. State your cleaning rule explicitly before continuing.

### YOUR CODE HERE
### inspect the smallest values of survey['shoe_size'] and clean the dataframe

We want to turn both continuous variables into three categories: Small, Medium, and Large.

For each variable, we divide the range from its minimum to its maximum into three intervals of equal width. Then we use pd.cut to assign each observation to one of the three categories. Finally, we store the resulting categories as new columns height_cat and shoe_cat in survey.

# Our labels for both categories
labels = ['Small', 'Medium', 'Large']

# Three bins have four edges: | 1 | 2 | 3 |
height_edges = np.linspace(
    survey['height_cm'].min(),
    survey['height_cm'].max(),
    4
)

# Three bins have four edges: | 1 | 2 | 3 |
shoe_edges = np.linspace(
    survey['shoe_size'].min(),
    survey['shoe_size'].max(),
    4
)

# Divide height into the bin categories
survey['height_cat'] = pd.cut(
    survey['height_cm'],
    bins=height_edges,
    labels=labels,
    include_lowest=True
)

# Divide shoe size into the bin categories
survey['shoe_cat'] = pd.cut(
    survey['shoe_size'],
    bins=shoe_edges,
    labels=labels,
    include_lowest=True
)

print('Height edges:', np.round(height_edges, 1))
print('Shoe-size edges:', np.round(shoe_edges, 1))

survey[['height_cm', 'height_cat', 'shoe_size', 'shoe_cat']].head()

Task 4.3: Before building the full contingency table, calculate one entry manually.

How many students have both:

  • height_cat == 'Small', and

  • shoe_cat == 'Small'?

Use Boolean indexing to count the observations satisfying both conditions. tore the result as n_small_small.

n_small_small = ### YOUR CODE HERE
print(n_small_small)

Task 4.4:

You just counted one combination of categories manually. Pandas can repeat this operation for every combination at once using .crosstab, which creates a full contingency table for you! The first argument defines the row categories and the second defines the column categories.

We use the categorized columns in survey to construct the full contingency table, with height categories as rows and shoe-size categories as columns.

Verify that the Small / Small entry agrees with your result from Task 4.4.

Then convert the counts into joint probabilities by dividing the entire table by the total number of observations:

\[P(A,B)=\frac{\mathrm{count}(A,B)}{n}.\]
# Count every combination of height and shoe-size category
counts = pd.crosstab(
    survey['height_cat'],
    survey['shoe_cat']
)

print("Raw data table:")
print(counts)

# Check one entry against the count you calculated manually in Task 4.4
assert counts.loc['Small', 'Small'] == n_small_small

# Convert counts to joint probabilities
joint = ### YOUR CODE HERE

print("Joint probabilities:")
joint.round(3)

Task 4.5:

The joint probability table describes height and shoe size jointly. Now compute the two marginal probability distributions:

  • the probability distribution of height category, regardless of shoe size;

  • the probability distribution of shoe-size category, regardless of height.

Compute each marginal distribution by summing the appropriate entries of the joint probability table. Check that each marginal distribution sums to 1.

Finally, think back to how we created the categories. We used intervals of equal width. If instead we had divided each variable into three groups using its empirical quantiles, approximately what would you expect the three marginal probabilities to be?

marginal_height = ### YOUR CODE HERE
marginal_shoe = ### YOUR CODE HERE

Task 4.6:

So far, the marginal distribution of shoe size described the whole class.

Now consider only students in the Small height category. Among these students, what is the probability distribution of shoe-size category?

Then repeat the calculation for students in the Large height category.

Use the joint and marginal probabilities from the previous tasks. Recall that conditioning means restricting the population we are considering and renormalizing the corresponding probabilities.

shoe_given_small_height = ### YOUR CODE HERE
shoe_given_large_height = ### YOUR CODE HERE

print('Shoe size | Small height')
print(shoe_given_small_height.round(3))

print('\nShoe size | Large height')
print(shoe_given_large_height.round(3))

Task 4.7:

Compare the marginal distribution of shoe size from Task 4.5 with the two conditional distributions of shoe size from Task 4.6.

If height category and shoe-size category were independent, what relationship would you expect between the marginal and conditional distributions?

Does that appear to be the case for these data? Support your conclusion with specific probabilities from your results.

Task 4.8 (optional / if you have time):

Try the same recipe (scatter → clean → bin into equal-width categories → contingency table → conditional probabilities) on another pair of variables in dataset_class_survey_2026.csv:

  • two more continuous variables, e.g. rent_eur vs. room_size_m2 (also needs cleaning — check for implausible values first!)

  • or two categorical variables, e.g. hair_color vs. eye_color (the original lecture example, for comparison), or cats_or_dogs vs. degree (equal-width binning doesn’t apply to these — they’re already categories) Which pairs show a strong dependence, and which look closer to independent?

### YOUR CODE HERE

End of Part 4.

You’ve now practiced the discrete counterpart of the fitting workflow from Parts 1-3: instead of fitting a continuous parametric distribution to one variable, you grouped two continuous variables into equal-width categories and built an empirical joint distribution directly from the counts, then used conditional probabilities to reason about dependence. On Friday, you’ll apply the same recipe — scatter plot, then bin into equal-width Low/Medium/High categories, then contingency table — to the two continuous variables in the group assignment.

By Max Ramgraber, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.