Assignment#
Part 1: Introduction and set up#
In this assignment, you will investigate the joint variability of three soil properties. Soil properties vary naturally from sample to sample due to differences in soil composition, depositional history, stress history, and water content. These properties may also be statistically dependent. Therefore, considering them jointly can provide information that is not captured when each property is analysed separately.
The dataset Soil Properties contains laboratory soil-characterisation data. In this exercise, we focus on three random variables:
Unit weight, \(\gamma\) (\(\mathrm{kN/m^3}\))
Water content, \(w\) (%)
Plasticity index, \(PI\) (%)
import numpy as np
from scipy.stats import multivariate_normal
from scipy.stats import norm
import matplotlib.pyplot as plt
from statistics import mean, stdev
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)
fname = 'Soil_Properties.csv'
findfile(fname)
# Load the three soil-property columns from the input file
raw = np.genfromtxt(fname, delimiter=',', names=True, dtype=float, encoding=None)
gamma = np.asarray(raw['gamma'], dtype=float)
w = np.asarray(raw['w'], dtype=float)
PI = np.asarray(raw['PI'], dtype=float)
data = np.column_stack((gamma, w, PI))
print(f'Number of complete observations = {len(data)}')
data
Part 2: Data exploration#
Task 2.1:
Plot the data. Examine the histogram of each soil property and the scatter plots for each pair of properties.
What can you conclude about:
the shape and variability of the marginal distributions; and
the direction and approximate strength of dependence between the soil properties?
variable_names = [
r'$\gamma$ (kN/m$^3$)',
r'$w$ (%)',
r'$PI$ (%)'
]
fig, axs = plt.subplots(3, 3, figsize=(8, 6), layout='constrained')
for i in range(data.shape[1]):
for j in range(data.shape[1]):
if i == j:
axs[i, j].grid()
axs[i, j].hist(
data[:, i],
bins=20,
density=True,
alpha=0.5,
edgecolor='black'
)
axs[i, j].set_xlabel(variable_names[i])
axs[i, j].set_ylabel('PDF')
elif i < j:
axs[i, j].grid()
axs[i, j].scatter(
data[:, i],
data[:, j],
alpha=0.5,
edgecolor='black'
)
axs[i, j].set_xlabel(variable_names[i])
axs[i, j].set_ylabel(variable_names[j])
else:
axs[i, j].axis('off')
plt.show()
Task 2.2:
Quantify the strength of the linear dependence between each pair using Pearson’s correlation coefficient. What can you conclude?
def calculate_covariance(X1, X2):
mean_x1 = mean(X1)
mean_x2 = mean(X2)
# Calculate deviations from the mean
diff_x1 = ### YOUR CODE HERE ###
diff_x2 = ### YOUR CODE HERE ###
# Multiply paired deviations
product = ### YOUR CODE HERE ###
# Sample covariance
covariance = ### YOUR CODE HERE ###
return covariance
def pearson_correlation(X1, X2):
covariance = ### YOUR CODE HERE ###
correl_coeff = ### YOUR CODE HERE ###
return correl_coeff
print("Pearson's correlation coefficient between γ and w is",
np.round(pearson_correlation(data[:,0], data[:,1]), 3))
print("Pearson's correlation coefficient between γ and PI is",
np.round(pearson_correlation(data[:,0], data[:,2]), 3))
print("Pearson's correlation coefficient between w and PI is",
np.round(pearson_correlation(data[:,1], data[:,2]), 3))
Part 3: Multivariate Gaussian distribution#
Task 3.1:
Model the joint distribution of \(\gamma\), \(w\) and \(PI\) using a multivariate Gaussian distribution. Follow the steps:
Define the vector of means.
Define the covariance matrix.
Draw 500 samples from the fitted trivariate Gaussian distribution and compare them with the observations.
Do you see differences between the observations and the Gaussian samples? Compare both the univariate marginal distributions and the bivariate scatter plots.
# ---------------------------------------------------------
# Define the vector of means
# ---------------------------------------------------------
mu1 = ### YOUR CODE HERE ###
mu2 = ### YOUR CODE HERE ###
mu3 = ### YOUR CODE HERE ###
mu = ### YOUR CODE HERE ###
# ---------------------------------------------------------
# Define the covariance matrix
# ---------------------------------------------------------
s1 = ### YOUR CODE HERE ###
s2 = ### YOUR CODE HERE ###
s3 = ### YOUR CODE HERE ###
covariance_12 = ### YOUR CODE HERE ###
covariance_13 = ### YOUR CODE HERE ###
covariance_23 = ### YOUR CODE HERE ###
sigma = ### YOUR CODE HERE ###
# ---------------------------------------------------------
# Draw 500 samples from the fitted trivariate Gaussian
# ---------------------------------------------------------
samples = ### YOUR CODE HERE ###
# ---------------------------------------------------------
# Compare samples with observations
# ---------------------------------------------------------
fig, axs = plt.subplots(
3, 3,
figsize=(8, 6),
layout='constrained'
)
for i in range(data.shape[1]):
for j in range(data.shape[1]):
# Marginal distributions
if i == j:
axs[i, j].grid()
axs[i, j].hist(
data[:, i],
bins=20,
density=True,
alpha=0.9,
edgecolor='black',
facecolor='none',
label='Observations'
)
axs[i, j].hist(
samples[:, i],
bins=20,
density=True,
alpha=0.5,
edgecolor='black',
label='Gaussian samples'
)
axs[i, j].set_xlabel(variable_names[i])
axs[i, j].set_ylabel('PDF')
axs[i, j].legend(fontsize=7)
# Scatter plots
elif i < j:
axs[i, j].grid()
# Gaussian samples first
axs[i, j].scatter(
samples[:, i],
samples[:, j],
alpha=0.30,
s=20,
label='Gaussian samples'
)
# Observations on top
axs[i, j].scatter(
data[:, i],
data[:, j],
alpha=0.9,
s=30,
edgecolor='black',
facecolor='none',
label='Observations'
)
axs[i, j].set_xlabel(variable_names[i])
axs[i, j].set_ylabel(variable_names[j])
# Put observations first in legend
handles, labels = axs[i, j].get_legend_handles_labels()
axs[i, j].legend(
[handles[1], handles[0]],
[labels[1], labels[0]],
fontsize=7,
markerscale=0.7
)
else:
axs[i, j].axis('off')
plt.show()
Task 3.2:
Suppose we are interested in soil with the following properties:
unit weight: \(\gamma < 16\ \mathrm{kN/m^3}\)
water content: \(w > 80\%\)
plasticity index: \(PI > 30\%\)
Using samples generated from the fitted multivariate Gaussian distribution, estimate the probability that all three conditions are satisfied at the same time.
Repeat the calculation using different numbers of samples,
\[ >N = 100,\quad 1{,}000,\quad 10{,}000,\quad 100{,}000,\quad 1{,}000{,}000. >\]How does the estimated probability behave as the number of samples increases? Comment on the stability of the estimate.
Hint: For each generated sample, check whether all three conditions are satisfied. The probability can be estimated as the proportion of generated samples that satisfy all three conditions.
# Investigate the effect of sample size
sample_sizes = [100, 1000, 10000, 100000, 1000000]
for N in sample_sizes:
# Draw samples from the fitted multivariate Gaussian distribution
samples = ### YOUR CODE HERE ###
# Check the three screening criteria simultaneously
condition = ### YOUR CODE HERE ###
# Estimate the probability as the proportion satisfying all three conditions
prob_meeting_condition = ### YOUR CODE HERE ###
print(f'N = {N:7d}: estimated probability = {prob_meeting_condition:.4f}')
Part 4: Conditioning the multivariate Gaussian distribution#
Task 4.1:
Suppose a laboratory measurement gives a water content of
\[ >w = 120\% >\]which is above the sample mean. What does this observation imply for the joint distribution of unit weight \(\gamma\) and plasticity index \(PI\)?
Compute the conditional distribution of \(\gamma\) and \(PI\) given \(w = 120\%\).
Compute the conditional probability
\[ >P\left( >\gamma < 16\ \mathrm{kN/m^3}, >\quad >PI > 30\% >\mid >w = 120\% >\right). >\]
# Conditional distribution of X = [γ, PI] given Y = w = 120%
a = 120.0
# Partition the mean vector for X = [gamma, PI] and Y = w
mu_X = np.array([
mu[0], # mean of gamma
mu[2] # mean of PI
])
mu_Y = mu[1] # mean of w
# Partition the covariance matrix
Sigma_XX = np.array([
[sigma[0, 0], sigma[0, 2]],
[sigma[2, 0], sigma[2, 2]]
])
Sigma_XY = np.array([
[sigma[0, 1]],
[sigma[2, 1]]
])
Sigma_YY = sigma[1, 1]
# ---------------------------------------------------------
# Conditional Gaussian equations
#
# mu_X|Y = mu_X + Sigma_XY * Sigma_YY^(-1) * (a - mu_Y)
#
# Sigma_X|Y = Sigma_XX
# - Sigma_XY * Sigma_YY^(-1) * Sigma_YX
# ---------------------------------------------------------
# Conditional mean vector
cond_mu_vector = ### YOUR CODE HERE ###
# Conditional covariance matrix
cond_sigma_matrix = ### YOUR CODE HERE ###
# Conditional standard deviations
cond_std = np.sqrt(np.diag(cond_sigma_matrix))
print('Conditional mean [γ, PI] =')
print(np.round(cond_mu_vector, 3))
print('\nConditional covariance matrix =')
print(np.round(cond_sigma_matrix, 3))
print('\nConditional standard deviations [γ, PI] =')
print(np.round(cond_std, 3))
# ---------------------------------------------------------
# Generate samples from the conditional Gaussian distribution
# ---------------------------------------------------------
N = 100000
cond_samples = ### YOUR CODE HERE ###
# ---------------------------------------------------------
# Estimate P(gamma < 16, PI > 30 | w = 120%)
# ---------------------------------------------------------
condition = ### YOUR CODE HERE ###
cond_prob = ### YOUR CODE HERE ###
print(
f'\nConditional probability '
f'P(γ < 16, PI > 30 | w = 120%) = {cond_prob:.4f}'
)
Task 4.2:
Using the results from the previous task, compare the conditional and unconditional distributions of \(\gamma\) and \(PI\).
Plot the conditional and unconditional marginal distributions of \(\gamma\) and \(PI\).
Plot the conditional and unconditional joint probability density contours of \(\gamma\) and \(PI\).
Compare the two cases and explain:
how the mean values change after conditioning on \(w = 120\%\);
how the uncertainty changes;
how these changes are related to the correlations of \(\gamma\) and \(PI\) with \(w\).
# Compare conditional and unconditional marginals of γ and PI
values_gamma = np.linspace(
max(0, mu[0] - 4*np.sqrt(sigma[0,0])),
mu[0] + 4*np.sqrt(sigma[0,0]),
200
)
uncond_γ_pdf = ### YOUR CODE HERE ###
cond_γ_pdf = ### YOUR CODE HERE ###
values_PI = np.linspace(
max(0, mu[2] - 4*np.sqrt(sigma[2,2])),
mu[2] + 4*np.sqrt(sigma[2,2]),
200
)
uncond_PI_pdf = ### YOUR CODE HERE ###
cond_PI_pdf = ### YOUR CODE HERE ###
# Plot marginal distributions
fig, ax = plt.subplots(1, 2, figsize=(9, 4), layout='constrained')
# γ
ax[0].plot(
values_gamma,
uncond_γ_pdf,
label='Unconditional distribution'
)
ax[0].plot(
values_gamma,
cond_γ_pdf,
'--',
label='Conditional distribution'
)
ax[0].grid()
ax[0].set_xlabel(r'$\gamma$ (kN/m$^3$)')
ax[0].set_ylabel('PDF')
ax[0].legend()
# PI
ax[1].plot(
values_PI,
uncond_PI_pdf,
label='Unconditional distribution'
)
ax[1].plot(
values_PI,
cond_PI_pdf,
'--',
label='Conditional distribution'
)
ax[1].grid()
ax[1].set_xlabel(r'$PI$ (%)')
ax[1].set_ylabel('PDF')
ax[1].legend()
plt.show()
# Compare joint conditional and unconditional distributions of γ and PI
# ---------------------------------------------------------
# 1. Create values for γ and PI
# ---------------------------------------------------------
values_gamma = np.linspace(9, 21, 200)
values_PI = np.linspace(0, 120, 200)
# Create a two-dimensional grid of γ and PI
G, I = np.meshgrid(values_gamma, values_PI)
# ---------------------------------------------------------
# 2. Define the unconditional distribution of [γ, PI]
# ---------------------------------------------------------
uncond_mu = [
mu[0], # mean of γ
mu[2] # mean of PI
]
uncond_cov = np.array([
[sigma[0, 0], sigma[0, 2]],
[sigma[2, 0], sigma[2, 2]]
])
# ---------------------------------------------------------
# 3. Calculate the joint PDFs
# ---------------------------------------------------------
# Create pairs of [γ, PI] values
X = ### YOUR CODE HERE ###
# Hint: use np.concatenate on G.T and I.T, then stack the two 1D arrays as columns.
# Unconditional joint PDF
Z_uncond = ### YOUR CODE HERE ###
# Conditional joint PDF given w = 120%
Z_cond = ### YOUR CODE HERE ###
# Convert the PDF values back to the two-dimensional grid
Z_uncond = ### YOUR CODE HERE ###
Z_cond = ### YOUR CODE HERE ###
# ---------------------------------------------------------
# 4. Plot the contours
# ---------------------------------------------------------
fig, ax = plt.subplots(figsize=(7, 5))
ax.contour(
G, I, Z_uncond,
levels=8,
linestyles='-'
)
ax.contour(
G, I, Z_cond,
levels=8,
linestyles='--'
)
ax.grid()
ax.set_xlabel(r'$\gamma$ (kN/m$^3$)')
ax.set_ylabel(r'$PI$ (%)')
ax.set_title(
r'Unconditional and conditional distributions ($w=120\%$)'
)
plt.show()
By Zheng Guan, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.