Analysis of Passing Maneuvers#
Part 1: Introduction and set up#
In this assignment, you will work with data on passing maneuvers from a fixed-based driving simulator experiment (you can read more about it here!) to assess the safety of these maneuvers. This type of analysis is key to better understand conflicts between drivers and better estimate possible crashes.
To do so, you will use data with three random variables:
The speed of the passing vehicle (\(S\), m/s)
Gap from the passing vehicle to the front of the passed vehicle at the end of the maneuver (\(d\), m)
The duration of the pass maneuver (\(P\), s)
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('http://files.mude.citg.tudelft.nl/'+fname, fname)
findfile('modif_data_passing.csv')
data = np.loadtxt('modif_data_passing.csv', delimiter = ',', skiprows=1)
data
Part 2: Data exploration#
Task 2.1:
Plot the data: have a look at the histograms of each variable as well as the scatter plots of each of them against the others. What can you conclude?
### YOUR CODE LINES HERE
Task 2.2:
Quantify the stregth of the correlation between each pair of variables using Pearson’s correlation coefficient. What can you conclude?
### YOUR CODE LINES HERE
Part 3: Multivariate Gaussian distribution#
Task 3.1:
Model the joint distribution of \(S\), \(d\) and \(P\) using a multivariate Gaussian distribution. Follow the steps:
Define the vector of means
Define the covariance matrix
Define the bivariate Gaussian distribution and draw 500 samples to compare with the observations.
Do you see differences between the observations and the samples from the bivariate Gaussian distribution? Compare both the univariate distributions and the bivariate scatter plots.
# Define the vector of means
mu1 = ### YOUR CODE HERE ###
mu2 = ### YOUR CODE HERE ###
mu3 = ### YOUR CODE HERE ###
mu = [mu1, mu2, mu3] #vector of means
# 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 a bivariate Gaussian distribution
samples = ### YOUR CODE HERE ###
# Plot against observations
### YOUR CODE LINES HERE
Task 3.2:
Now that you have a model, you want to get insights into how safe the circulation is on the road you are studying. In previous studies in a similar area, they have defined the following safety thresholds \(^*\):
High speeds lead to reduced driver’s reaction time and increased stopping distances. Therefore, \(S < 30 \ m/s\) is recommended.
Reduced distances between the passing and passed vehicle after the maneuver may indicate more risky maneuvers increasing the chance of collision. Thus, \(d>45 \ m\) is recommended.
Long passing times may lead to dangerous situations as it increases the likelihood of the passing maneuver ending in a no-passing zone due to traffic, road conditions etc. Hence, \(P<7 \ s\) is recommended.
Given the above safety thresholds, what is the probability of meeting them?
Hint: It can be computed analytically similarly to the 2-dimensional case. However, it is easier to use simulations to compute this probability.
\(^*\) These thresholds are not based in actual evidence and are just for academic purposes.
# Draw 10000 samples from the multivariate Gaussian distribution
samples = ### YOUR CODE LINES HERE
# Count number of points that meet the condition
condition = ### YOUR CODE LINES HERE
# Compute the probability of the condition
prob_meeting_condition = ### YOUR CODE LINES HERE
print(f"The probability of meeting the safety limits (S < 30, d > 45 and P < 7) is approximately {prob_meeting_condition:.4f}")
Part 4: Conditionalizing the multivariate Gaussian distribution#
Task 4.1:
You have now installed a sensor that measures the distance between the passed and passing vehicle at the end of the maneuver (\(d\)). For a passing vehicle, \(d=53 \ m\). What is the probability that it was a safe maneuver?
To answer the question, follow the steps:
Compute the distribution of \(S\) and \(P\) given \(d=53 \ m\).
Compute the probability of \(P(S \leq 30 \ m/s, P \leq 7 \ s|d=53 \ m)\).
cond_mu_vector = ### YOUR CODE LINES HERE
cond_sigma_matrix = ### YOUR CODE LINES HERE
Task 4.2:
Using the results of the previous task, compare the conditional and unconditional distributions of \(S\) and \(P\). To do so, follow the next steps:
Plot the conditional and unconditional univariate marginal distributions of \(S\) and \(P\).
Plot the contours of the PDF of the conditional and unconditional joint distributions of \(S\) and \(P\).
Explain the differences in terms of both the location and the width of the distribution.
# Compare the conditional and unconditional margins
# Evaluate both the conditional and unconditional univariate PDFs of S and P
### YOUR CODE LINES HERE
# Plot both PDFs
### YOUR CODE LINES HERE
# Compare the joint conditional and unconditional distribution
# Define the mesh of values where we want to evaluate the random variables
#already done in previous cell
# Define the grid
X1,X2 = ### YOUR CODE LINES HERE
X = ### YOUR CODE LINES HERE
# Evaluate the CDF
Z_uncond = ### YOUR CODE LINES HERE
Z_cond = ### YOUR CODE LINES HERE
# Create contours plot
fig = plt.figure(figsize=(6,5))
ax = fig.add_subplot(111)
contours = ax.contour(X1, X2, Z_uncond.reshape(X1.shape), 25, cmap='YlGnBu', vmin=0, vmax = 0.03, label = 'Unconditional joint PDF')
contours = ax.contour(X1, X2, Z_cond.reshape(X1.shape), 25, cmap='YlGnBu', vmin=0, vmax = 0.03, label = 'Conditional joint PDF')
ax.grid()
ax.set_xlim([0, 50])
ax.set_xlim([0, 50])
ax.set_ylim([0, 14])
ax.set_xlabel('${S} [m/s]$')
ax.set_ylabel('${P} [s]$')
fig.colorbar(contours, ax=ax, label='Probability Density')
By Patricia Mares Nasarre, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.