# Paris Saclay Center for Data Science¶

# RAMP on predicting the number of air passengers

Balázs Kégl (LAL/CNRS), Alex Gramfort (Inria), Djalel Benbouzid (UPMC), Mehdi Cherti (LAL/CNRS)

## Introduction¶

The data set was donated to us by an unnamed company handling flight ticket reservations. The data is thin, it contains

- the date of departure
- the departure airport
- the arrival airport
- the mean and standard deviation of the number of weeks of the reservations made before the departure date
- a field called `log_PAX` which is related to the number of passengers (the actual number were changed for privacy reasons)

The goal is to predict the `log_PAX` column. The prediction quality is measured by RMSE.

The data is obviously limited, but since data and location informations are available, it can be joined to external data sets. The challenge in this RAMP is to find good data that can be correlated to flight traffic.

In [1]:

    %matplotlib inline
    import os
    import importlib
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    pd.set_option('display.max_columns', None)

## Load the dataset using pandas¶

The training and testing data are located in the folder `data`. They are compressed `csv` file (i.e. `csv.bz2`). We can load the dataset using pandas.

In [2]:

    data = pd.read_csv(
        os.path.join('data', 'train.csv.bz2')
    )

In [3]:

    data.info()

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 8902 entries, 0 to 8901
    Data columns (total 6 columns):
    DateOfDeparture     8902 non-null object
    Departure           8902 non-null object
    Arrival             8902 non-null object
    WeeksToDeparture    8902 non-null float64
    log_PAX             8902 non-null float64
    std_wtd             8902 non-null float64
    dtypes: float64(3), object(3)
    memory usage: 417.4+ KB

So as stated earlier, the column `log_PAX` is the target for our regression problem. The other columns are the features which will be used for our prediction problem. If we focus on the data type of the column, we can see that `DateOfDeparture`, `Departure`, and `Arrival` are of `object` dtype, meaning they are strings.

In [4]:

    data[['DateOfDeparture', 'Departure', 'Arrival']].head()

Out[4]:

| DateOfDeparture | Departure | Arrival
--- | --- | --- | ---
0 | 2012-06-19 | ORD | DFW
1 | 2012-09-10 | LAS | DEN
2 | 2012-10-05 | DEN | LAX
3 | 2011-10-09 | ATL | ORD
4 | 2012-02-21 | DEN | SFO

While it makes `Departure` and `Arrival` are the code of the airport, we see that the `DateOfDeparture` should be a date instead of string. We can use pandas to convert this data.

In [5]:

    data.loc[:, 'DateOfDeparture'] = pd.to_datetime(data.loc[:, 'DateOfDeparture'])

In [6]:

    data.info()

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 8902 entries, 0 to 8901
    Data columns (total 6 columns):
    DateOfDeparture     8902 non-null datetime64[ns]
    Departure           8902 non-null object
    Arrival             8902 non-null object
    WeeksToDeparture    8902 non-null float64
    log_PAX             8902 non-null float64
    std_wtd             8902 non-null float64
    dtypes: datetime64[ns](1), float64(3), object(2)
    memory usage: 417.4+ KB

When you will create a submission, `ramp-workflow` will load the data for you and split into a data matrix `X` and a target vector `y`. It will also take care about splitting the data into a training and testing set. These utilities are available in the module `problem.py` which we will load.

In [7]:

    import problem

The function `get_train_data()` loads the training data and returns a pandas dataframe `X` and a numpy vector `y`.

In [8]:

    X, y = problem.get_train_data()

In [9]:

    type(X)

Out[9]:

    pandas.core.frame.DataFrame

In [10]:

    type(y)

Out[10]:

    numpy.ndarray

We can check the information of the data `X`

In [11]:

    X.info()

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 8902 entries, 0 to 8901
    Data columns (total 5 columns):
    DateOfDeparture     8902 non-null object
    Departure           8902 non-null object
    Arrival             8902 non-null object
    WeeksToDeparture    8902 non-null float64
    std_wtd             8902 non-null float64
    dtypes: float64(2), object(3)
    memory usage: 347.9+ KB

Thus, this is important to see that `ramp-workflow` does not convert the `DateOfDeparture` column into a `datetime` format. Thus, keep in mind that you might need to make a conversion when prototyping your machine learning pipeline later on. Let's check some statistics regarding our dataset.

In [12]:

    print(min(X['DateOfDeparture']))
    print(max(X['DateOfDeparture']))

    2011-09-01
    2013-03-05

In [13]:

    X['Departure'].unique()

Out[13]:

    array(['ORD', 'LAS', 'DEN', 'ATL', 'SFO', 'EWR', 'IAH', 'LAX', 'DFW',
           'SEA', 'JFK', 'PHL', 'MIA', 'DTW', 'BOS', 'MSP', 'CLT', 'MCO',
           'PHX', 'LGA'], dtype=object)

In [14]:

    _ = plt.hist(y, bins=50)

In [15]:

    _ = X.hist('std_wtd', bins=50)

In [16]:

    _ = X.hist('WeeksToDeparture', bins=50)

In [17]:

    X.describe()

Out[17]:

| WeeksToDeparture | std_wtd
--- | --- | ---
count | 8902.000000 | 8902.000000
mean | 11.446469 | 8.617773
std | 2.787140 | 2.139604
min | 2.625000 | 2.160247
25% | 9.523810 | 7.089538
50% | 11.300000 | 8.571116
75% | 13.240000 | 10.140521
max | 23.163265 | 15.862216

In [18]:

    X.shape

Out[18]:

    (8902, 5)

In [19]:

    print(y.mean())
    print(y.std())

    10.99904767212102
    0.9938894125318564

## Preprocessing dates¶

Getting dates into numerical columns is a common operation when time series data is analyzed with non-parametric predictors. The code below makes the following transformations:

- numerical columns for year (2011-2012), month (1-12), day of the month (1-31), day of the week (0-6), and week of the year (1-52)
- number of days since 1970-01-01

In [20]:

    # Make a copy of the original data to avoid writing on the original data
    X_encoded = X.copy()

    # following http://stackoverflow.com/questions/16453644/regression-with-date-variable-using-scikit-learn
    X_encoded['DateOfDeparture'] = pd.to_datetime(X_encoded['DateOfDeparture'])
    X_encoded['year'] = X_encoded['DateOfDeparture'].dt.year
    X_encoded['month'] = X_encoded['DateOfDeparture'].dt.month
    X_encoded['day'] = X_encoded['DateOfDeparture'].dt.day
    X_encoded['weekday'] = X_encoded['DateOfDeparture'].dt.weekday
    X_encoded['week'] = X_encoded['DateOfDeparture'].dt.week
    X_encoded['n_days'] = X_encoded['DateOfDeparture'].apply(lambda date: (date - pd.to_datetime("1970-01-01")).days)

In [21]:

    X_encoded.tail(5)

Out[21]:

| DateOfDeparture | Departure | Arrival | WeeksToDeparture | std_wtd | year | month | day | weekday | week | n_days
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
8897 | 2011-10-02 | DTW | ATL | 9.263158 | 7.316967 | 2011 | 10 | 2 | 6 | 39 | 15249
8898 | 2012-09-25 | DFW | ORD | 12.772727 | 10.641034 | 2012 | 9 | 25 | 1 | 39 | 15608
8899 | 2012-01-19 | SFO | LAS | 11.047619 | 7.908705 | 2012 | 1 | 19 | 3 | 3 | 15358
8900 | 2013-02-03 | ORD | PHL | 6.076923 | 4.030334 | 2013 | 2 | 3 | 6 | 5 | 15739
8901 | 2011-11-26 | DTW | ATL | 9.526316 | 6.167733 | 2011 | 11 | 26 | 5 | 47 | 15304

We will perform all preprocessing steps within a scikit-learn pipeline which chains together tranformation and estimator steps. This offers offers convenience and safety (help avoid leaking statistics from your test data into the trained model in cross-validation) and the whole pipeline can be evaluated with `cross_val_score`.

To perform the above encoding within a scikit-learn pipeline we will a function and using `FunctionTransformer` to make it compatible with scikit-learn API.

In [22]:

    from sklearn.preprocessing import FunctionTransformer

    def _encode_dates(X):
        # With pandas < 1.0, we wil get a SettingWithCopyWarning
        # In our case, we will avoid this warning by triggering a copy
        # More information can be found at:
        # https://github.com/scikit-learn/scikit-learn/issues/16191
        X_encoded = X.copy()

        # Make sure that DateOfDeparture is of datetime format
        X_encoded.loc[:, 'DateOfDeparture'] = pd.to_datetime(X_encoded['DateOfDeparture'])
        # Encode the DateOfDeparture
        X_encoded.loc[:, 'year'] = X_encoded['DateOfDeparture'].dt.year
        X_encoded.loc[:, 'month'] = X_encoded['DateOfDeparture'].dt.month
        X_encoded.loc[:, 'day'] = X_encoded['DateOfDeparture'].dt.day
        X_encoded.loc[:, 'weekday'] = X_encoded['DateOfDeparture'].dt.weekday
        X_encoded.loc[:, 'week'] = X_encoded['DateOfDeparture'].dt.week
        X_encoded.loc[:, 'n_days'] = X_encoded['DateOfDeparture'].apply(
            lambda date: (date - pd.to_datetime("1970-01-01")).days
        )
        # Once we did the encoding, we will not need DateOfDeparture
        return X_encoded.drop(columns=["DateOfDeparture"])

    date_encoder = FunctionTransformer(_encode_dates)

In [23]:

    date_encoder.fit_transform(X).head()

Out[23]:

| Departure | Arrival | WeeksToDeparture | std_wtd | year | month | day | weekday | week | n_days
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | ORD | DFW | 12.875000 | 9.812647 | 2012 | 6 | 19 | 1 | 25 | 15510
1 | LAS | DEN | 14.285714 | 9.466734 | 2012 | 9 | 10 | 0 | 37 | 15593
2 | DEN | LAX | 10.863636 | 9.035883 | 2012 | 10 | 5 | 4 | 40 | 15618
3 | ATL | ORD | 11.480000 | 7.990202 | 2011 | 10 | 9 | 6 | 40 | 15256
4 | DEN | SFO | 11.450000 | 9.517159 | 2012 | 2 | 21 | 1 | 8 | 15391

## Random Forests¶

Tree-based algorithms requires less complex preprocessing than linear-models. We will first present a machine-learning pipeline where we will use a random forest. In this pipeline, we will need to:

- encode the date to numerical values (as presented in the section above);
- oridinal encode the other categorical values to get numerical number;
- keep numerical features as they are.

Thus, we want to perform three different processes on different columns of the original data `X`. In scikit-learn, we can use `make_column_transformer` to perform such processing.

In [24]:

    from sklearn.preprocessing import OrdinalEncoder
    from sklearn.compose import make_column_transformer

    date_encoder = FunctionTransformer(_encode_dates)
    date_cols = ["DateOfDeparture"]

    categorical_encoder = OrdinalEncoder()
    categorical_cols = ["Arrival", "Departure"]

    preprocessor = make_column_transformer(
        (date_encoder, date_cols),
        (categorical_encoder, categorical_cols),
        remainder='passthrough',  # passthrough numerical columns as they are
    )

We can combine our preprocessor with an estimator (`RandomForestRegressor` in this case), allowing us to make predictions.

In [25]:

    from sklearn.ensemble import RandomForestRegressor
    from sklearn.pipeline import make_pipeline

    n_estimators = 10
    max_depth = 10
    max_features = 10

    regressor = RandomForestRegressor(
        n_estimators=n_estimators, max_depth=max_depth, max_features=max_features
    )

    pipeline = make_pipeline(preprocessor, regressor)

We can cross-validate our `pipeline` using `cross_val_score`. Below we will have specified `cv=5` meaning KFold cross-valdiation splitting will be used, with 8 folds. The mean squared error regression loss is calculated for each split. The output `score` will be an array of 5 scores from each KFold. The score mean and standard deviation of the 5 scores is printed at the end.

In [26]:

    from sklearn.model_selection import cross_val_score

    scores = cross_val_score(
        pipeline, X, y, cv=5, scoring='neg_mean_squared_error'
    )
    rmse_scores = np.sqrt(-scores)

    print(
        f"RMSE: {np.mean(rmse_scores):.4f} +/- {np.std(rmse_scores):.4f}"
    )

    RMSE: 0.6277 +/- 0.0186

## Linear regressor¶

When dealing with a linear model, we need to one-hot encode categorical variables instead of ordinal encoding and standardize numerical variables. Thus we will:

- encode the date;
- then, one-hot encode all categorical columns, including the encoded date as well;
- standardize the numerical columns.

In [27]:

    from sklearn.preprocessing import OneHotEncoder
    from sklearn.preprocessing import StandardScaler

    date_encoder = FunctionTransformer(_encode_dates)
    date_cols = ["DateOfDeparture"]

    categorical_encoder = OneHotEncoder(handle_unknown="ignore")
    categorical_cols = [
        "Arrival", "Departure", "year", "month", "day",
        "weekday", "week", "n_days"
    ]

    numerical_scaler = StandardScaler()
    numerical_cols = ["WeeksToDeparture", "std_wtd"]

    preprocessor = make_column_transformer(
        (categorical_encoder, categorical_cols),
        (numerical_scaler, numerical_cols)
    )

We can now combine our `preprocessor` with the `LinearRegression` estimator in a `Pipeline`:

In [28]:

    from sklearn.linear_model import LinearRegression

    regressor = LinearRegression()

    pipeline = make_pipeline(date_encoder, preprocessor, regressor)

And we can evaluate our linear-model pipeline:

In [29]:

    scores = cross_val_score(
        pipeline, X, y, cv=5, scoring='neg_mean_squared_error'
    )
    rmse_scores = np.sqrt(-scores)

    print(
        f"RMSE: {np.mean(rmse_scores):.4f} +/- {np.std(rmse_scores):.4f}"
    )

    RMSE: 0.6117 +/- 0.0149

# Merging external data¶

The objective in this RAMP data challenge is to find good data that can be correlated to flight traffic. We will use some weather data (saved in `submissions/starting_kit`) to provide an example of how to merge external data in a scikit-learn pipeline.

Your external data will need to be included in your submissions folder - see RAMP submissions for more details.

First we will define a function that merges the external data to our feature data.

In [30]:

    # when submitting a kit, the `__file__` variable will corresponds to the
    # path to `estimator.py`. However, this variable is not defined in the
    # notebook and thus we must define the `__file__` variable to imitate
    # how a submission `.py` would work.
    __file__ = os.path.join('submissions', 'starting_kit', 'estimator.py')
    filepath = os.path.join(os.path.dirname(__file__), 'external_data.csv')
    filepath

Out[30]:

    'submissions/starting_kit/external_data.csv'

In [31]:

    pd.read_csv(filepath).head()

Out[31]:

| Date | AirPort | Max TemperatureC | Mean TemperatureC | Min TemperatureC | Dew PointC | MeanDew PointC | Min DewpointC | Max Humidity | Mean Humidity | Min Humidity | Max Sea Level PressurehPa | Mean Sea Level PressurehPa | Min Sea Level PressurehPa | Max VisibilityKm | Mean VisibilityKm | Min VisibilitykM | Max Wind SpeedKm/h | Mean Wind SpeedKm/h | Max Gust SpeedKm/h | Precipitationmm | CloudCover | Events | WindDirDegrees
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 2011-09-01 | ATL | 35 | 29 | 24 | 21 | 18 | 14 | 79 | 56 | 32 | 1022 | 1019 | 1017 | 16 | 16 | 11 | 19 | 6 | 26.0 | 0.00 | 3 | NaN | 129
1 | 2011-09-02 | ATL | 36 | 29 | 22 | 17 | 15 | 14 | 61 | 46 | 30 | 1019 | 1016 | 1014 | 16 | 16 | 16 | 24 | 7 | 34.0 | 0.00 | 2 | NaN | 185
2 | 2011-09-03 | ATL | 35 | 29 | 23 | 17 | 16 | 14 | 64 | 47 | 30 | 1015 | 1013 | 1011 | 16 | 16 | 16 | 19 | 7 | 26.0 | 0.00 | 4 | NaN | 147
3 | 2011-09-04 | ATL | 27 | 24 | 22 | 22 | 19 | 16 | 93 | 72 | 51 | 1014 | 1012 | 1011 | 16 | 14 | 4 | 21 | 9 | 26.0 | 6.10 | 6 | Rain | 139
4 | 2011-09-05 | ATL | 26 | 24 | 22 | 23 | 22 | 20 | 94 | 91 | 87 | 1010 | 1005 | 999 | 16 | 13 | 3 | 32 | 16 | 45.0 | 16.00 | 8 | Rain-Thunderstorm | 149

In [32]:

    def _merge_external_data(X):
        filepath = os.path.join(
            os.path.dirname(__file__), 'external_data.csv'
        )

        X = X.copy()  # to avoid raising SettingOnCopyWarning
        # Make sure that DateOfDeparture is of dtype datetime
        X.loc[:, "DateOfDeparture"] = pd.to_datetime(X['DateOfDeparture'])
        # Parse date to also be of dtype datetime
        data_weather = pd.read_csv(filepath, parse_dates=["Date"])

        X_weather = data_weather[['Date', 'AirPort', 'Max TemperatureC']]
        X_weather = X_weather.rename(
            columns={'Date': 'DateOfDeparture', 'AirPort': 'Arrival'})
        X_merged = pd.merge(
            X, X_weather, how='left', on=['DateOfDeparture', 'Arrival'], sort=False
        )
        return X_merged

    data_merger = FunctionTransformer(_merge_external_data)

Double check that our function works:

In [33]:

    data_merger.fit_transform(X).head()

Out[33]:

| DateOfDeparture | Departure | Arrival | WeeksToDeparture | std_wtd | Max TemperatureC
--- | --- | --- | --- | --- | --- | ---
0 | 2012-06-19 | ORD | DFW | 12.875000 | 9.812647 | 34
1 | 2012-09-10 | LAS | DEN | 14.285714 | 9.466734 | 33
2 | 2012-10-05 | DEN | LAX | 10.863636 | 9.035883 | 22
3 | 2011-10-09 | ATL | ORD | 11.480000 | 7.990202 | 27
4 | 2012-02-21 | DEN | SFO | 11.450000 | 9.517159 | 16

Use `FunctionTransformer` to make our function compatible with scikit-learn API:

We can now assemble our pipeline using the same `data_merger` and `preprocessor` as above:

In [34]:

    from sklearn.preprocessing import OrdinalEncoder
    from sklearn.compose import make_column_transformer

    date_encoder = FunctionTransformer(_encode_dates)
    date_cols = ["DateOfDeparture"]

    categorical_encoder = OrdinalEncoder()
    categorical_cols = ["Arrival", "Departure"]

    preprocessor = make_column_transformer(
        (date_encoder, date_cols),
        (categorical_encoder, categorical_cols),
        remainder='passthrough',  # passthrough numerical columns as they are
    )

In [35]:

    n_estimators = 10
    max_depth = 10
    max_features = 10

    regressor = RandomForestRegressor(
        n_estimators=n_estimators, max_depth=max_depth, max_features=max_features
    )

    pipeline = make_pipeline(data_merger, preprocessor, regressor)

In [36]:

    scores = cross_val_score(
        pipeline, X, y, cv=5, scoring='neg_mean_squared_error'
    )
    rmse_scores = np.sqrt(-scores)

    print(
        f"RMSE: {np.mean(rmse_scores):.4f} +/- {np.std(rmse_scores):.4f}"
    )

    RMSE: 0.6330 +/- 0.0108

## Feature importances¶

We can check the feature importances using the function `sklearn.inspection.permutation_importances`. Since the first step of our pipeline adds the new external feature `Max TemperatureC`, we want to apply this transformation after adding `Max TemperatureC`, to check the importances of all features. Indeed, we can perform `sklearn.inspection.permutation_importances` at any stage of the pipeline, as we will see later on.

The code below:

- performs `transform` on the first step of the pipeline (`pipeline[0]`) producing the transformed train (`X_train_augmented`) and test (`X_test_augmented`) data
- the transformed data is used to fit the pipeline from the second step onwards

Note that pipelines can be sliced. `pipeline[0]` obtains the first step (tuple) of the pipeline. You can further slice to obtain either the transformer/estimator (first item in each tuple) or column list (second item within each tuple) inside each tuple. For example `pipeline[0][0]` obtains the transformer of the first step of the pipeline (first item of the first tuple).

In [37]:

    from sklearn.model_selection import train_test_split

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, random_state=42
    )

    merger = pipeline[0]
    X_train_augmented = merger.transform(X_train)
    X_test_augmented = merger.transform(X_test)

    predictor = pipeline[1:]
    predictor.fit(X_train_augmented, y_train).score(X_test_augmented, y_test)

Out[37]:

    0.5939918926708109

With the fitted pipeline, we can now use `permutation_importance` to calculate feature importances:

In [38]:

    from sklearn.inspection import permutation_importance

    feature_importances = permutation_importance(
        predictor, X_train_augmented, y_train, n_repeats=10
    )

Here, we plot the permutation importance using the training set. The higher the value, more important the feature is.

In [39]:

    sorted_idx = feature_importances.importances_mean.argsort()

    fig, ax = plt.subplots()
    ax.boxplot(feature_importances.importances[sorted_idx].T,
               vert=False, labels=X_train_augmented.columns[sorted_idx])
    ax.set_title("Permutation Importances (train set)")
    fig.tight_layout()
    plt.show()

We can replicate the same processing on the test set and see if we can observe the same trend.

In [40]:

    from sklearn.inspection import permutation_importance

    feature_importances = permutation_importance(
        predictor, X_test_augmented, y_test, n_repeats=10
    )

In [41]:

    sorted_idx = feature_importances.importances_mean.argsort()

    fig, ax = plt.subplots()
    ax.boxplot(feature_importances.importances[sorted_idx].T,
               vert=False, labels=X_test_augmented.columns[sorted_idx])
    ax.set_title("Permutation Importances (test set)")
    fig.tight_layout()
    plt.show()

With the current version of scikit-learn, it is not handy but still possible to check the feature importances at the latest stage of the pipeline (once all features have been preprocessed).

The difficult part is to get the name of the features.

In [42]:

    preprocessor = pipeline[:-1]
    predictor = pipeline[-1]

    X_train_augmented = preprocessor.transform(X_train)
    X_test_augmented = preprocessor.transform(X_test)

Let's find out the feature names (in the future, scikit-learn will provide a `get_feature_names` function to handle this case).

In [43]:

    date_cols_name = (date_encoder.transform(X_train[date_cols])
                                  .columns.tolist())
    categorical_cols_name = categorical_cols
    numerical_cols_name = (pipeline[0].transform(X_train)
                                      .columns[pipeline[1].transformers_[-1][-1]]
                                      .tolist())

In [44]:

    feature_names = np.array(
        date_cols_name + categorical_cols_name + numerical_cols_name
    )
    feature_names

Out[44]:

    array(['year', 'month', 'day', 'weekday', 'week', 'n_days', 'Arrival',
           'Departure', 'WeeksToDeparture', 'std_wtd', 'Max TemperatureC'],
          dtype='<U16')

We can repeat the previous processing at this finer grain, where the transformed date columns are included.

In [45]:

    from sklearn.inspection import permutation_importance

    feature_importances = permutation_importance(
        predictor, X_train_augmented, y_train, n_repeats=10
    )

Here, we plot the permutation importance using the training set. Basically, higher the value, more important is the feature.

In [46]:

    sorted_idx = feature_importances.importances_mean.argsort()

    fig, ax = plt.subplots()
    ax.boxplot(feature_importances.importances[sorted_idx].T,
               vert=False, labels=feature_names[sorted_idx])
    ax.set_title("Permutation Importances (train set)")
    fig.tight_layout()
    plt.show()

We can replicate the same processing on the test set and see if we can observe the same trend.

In [47]:

    from sklearn.inspection import permutation_importance

    feature_importances = permutation_importance(
        predictor, X_test_augmented, y_test, n_repeats=10
    )

In [48]:

    sorted_idx = feature_importances.importances_mean.argsort()

    fig, ax = plt.subplots()
    ax.boxplot(feature_importances.importances[sorted_idx].T,
               vert=False, labels=feature_names[sorted_idx])
    ax.set_title("Permutation Importances (test set)")
    fig.tight_layout()
    plt.show()

## Submission¶

To submit your code, you can refer to the online documentation.