Modality refers to a certain type of information and/or the representation format in which information is stored.

According to wikipedia, « Multimodal learning attempts to model the combination of different modalities of data, often rising in real-world applications. An example of joint data is combining text (typically represented as discrete word count vectors) with imaging data consisting of pixel intensities and annotation tags. As these modalities have fundamentally different statistical properties, combining them is non-trivial, which is why specialized modelling strategies and algorithms are required. »

Can we simultaneously train structured (tabular) and unstructured (text, image) data in the same neural network model while optimizing the same output target ?

We’ll perform this experiment using a very simple example by implementing a « concat » combination method that will concatenate text features, tabular features, and image features all at once before our final regressor layer(s).

Specifically, we will build a multimodal neural network that contains 3 input heads :

• Image : Inception V3 input layers (Transfer Learning from Imagenet)
• Text : Bidirectionnal LSTM with Attention layers pretrained on GlOve embeddings
• Structured : Multi-layers perceptron

And 1 output head :

• House Price

Of course, more advanced methods and strategies exist to solve multimodal problems. For example, the great HuggingFace offers a multimodal-toolkit that you can experiment with.

Table of contents

• Dataset
• Tabular data engineering
• Text cleaning
• Image cleaning — Fine-tune VGG16 to predict good/bad images
• Architecture#1 — Structured Data fine-tuning with MLP (Python code)
• Architecture#2 — Text Data fine-tuning using GlOve BiLSTM with Attention (Python code)
• Architecture#3 — Image Data fine-tuning with Inception V3 (Python code)
• Final Architecture — Load each previously fine-tuned weights and create the final hybrid neural networks training architecture by using the concat method. (Python Code)
• Results and Conclusion
• Final words

Dataset

First, let’s take a quick look at the dataset we’re going to use. I searched a very long time to find a dataset containing images, text and tabular data and finally found one! The link to download the dataset is here: https://www.kaggle.com/datasets/ericpierce/austinhousingprices (You must have a Kaggle account to download it).

As it is written on the Kaggle page:

“Context — Austin Housing

This data was originally sourced in support of my capstone project at Northwestern University. The Austin Housing market is one of the hottest markets in 2021, and these listings show how that market has changed over the past couple of years.”

The data contains a csv file with all the tabular features, a text description feature and a column that specify the image filename associated with the dataset row. The csv dataset comes with a huge zipped file containing the house images. Unfortunatly, there is only one image available for each house but we will work with that.

The goal here is to predict house price using all the data we have.

Tabular data engineering

Although it is said that the data has been cleaned, there are still many data quality issues left. For example, there are impossible feature values like houses with 0 bathrooms or 0 bedrooms or negative square foot of living area.

Here is a table that resume some of the major data transformations and cleaning steps performed:

An important step is to clean up our target value, house prices, using inflation rates (concept drift problem). Past selling prices do not match actual selling prices due to inflation; the definition of 1$ unit is not the same over time. For this purpose, I used a table of inflation rates to correct the price amount.

Once cleaned, our dataset contains 56 features and 10537 rows.

Here is a geospatial plot (plotly) of all Austin Houses we have:

We can see that there are some zipcode clusters where the price is higher (red) than other regions (blue). I’ve created a new feature that represent the zipcode clusters called regionCluster using the kmeans algorithm:

Text cleaning

To clean up house text descriptions, we apply normalization, remove stop words and find the root of each word using lemmatization.

Here is an example of the result:

Image cleaning — Fine-tune VGG16 to predict good/bad images

This one is a bit tricky. After randomly exploring some images, there seem to be very different / heterogenous types of house images. Here are a few:

As you can see, there are several heterogeneous types of house images. Another challenge is that we only have one image available per house, so the image must be useful.

Since there are about 10,000 images, we can’t look at all of them to decide which data to keep and which data to delete. To help us, we will fine-tune a VGG16 to predict good/bad images for us.

First, we need to determine what is a good image for us. I decided to keep only images where you can see the complete architecture of the house from the outside:

Bad images would be any interior image or any image where we cannot see the full architecture of the house from the outside.

As we are using transfer learning, fewer samples from each class should do the job.

We need to create two directories, one for the good images and one for the bad ones, then manually dispatch some of them in the right place. After an hour, I sorted around 800 images for the train set (400 good and 400 bad) and around 250 images for the test set.

We are now ready to fine-tune our VGG16 model!

Python Code:

Here is the final confusion matrix from the test data sample:

We can now predict every image we have using our “VGG16 cleaner” and remove all rows of our dataset where the available image has been classified as “bad”.

For the next steps, the strategy we will use is to individually train each neural network, save the weights and then load the trained weights into our final hybrid model and retrain with all the parts merged together.

Architecture#1 — Structured Data fine-tuning with MLP:

Python code to create the neural network architecture, fine-tune, and save the weights:

Multi-layer perceptron (MLP) architecture:

Here are the most important shapley values for our tabular features (Using the DeepLift algorithm (DeepExplainer) from the awesome SHAP package) :

This barplot shows us that the 3 most important features to predict price of our structured data neural network are: livingAreaSqFt, zipcode and regionCluster.

numPriceChanges, yearBuilt and lotSizeSqFt follow next.

If we look at some of these features:

The biggest livingAreaSqFt, the biggest shapley value, the biggest price (seems logical!)

The more recent yearBuilt, the biggest shapley value, the biggest price (seems logical!)

Architecture#2 — Text Data fine-tuning GlOve BiLSTM with Attention:

Python code to create the neural network architecture, fine-tune, and save the weights:

BiLSTM Attention with GloVe embeddings architecture:

Again using shapley values, here are some useful ngrams of the house text description to predict price:

The “large private backyard”, “open live area” and “wood laminate floor” grams appear to have the most positive influence on price according to the shapley values.

In contrast, “bright open floor”, “beautiful wood floor”, “yes sprinkler sy” and “vinyl plank floor” seem to have the most negative influence on the price.

Architecture#3 — Image Data fine-tuning with Inception V3:

We use the pre-trained InceptionV3 model on generic image (imagenet weights) and fine-tune only the head layer using our houses prices.

Python code to create the neural network architecture, fine-tune, and save the weights:

The final architecture (InceptionV3 is presented here as a single layer but represents a more complex underlying architecture):

Using the gradcam method, we can try to understand what our CNN is looking at:

The neural network seems to focus mostly on the house part of the images, which is fine, but there is no evidence which parts of the houses are more important than others… maybe our model is not trained enough !

Grad-CAM, Gradient-weighted Class Activation Mapping, uses the feature maps / class-specific gradient information flowing into the final convolutional layer of a CNN. It produces a visual explanation heat-maps / localization map that will help us understand how the model make decisions / the important regions in the image.

Python code for the gradcam:

Final Architecture — Load fine-tuned weights and create final hybrid architecture for training

Here is the python code to load all our fine-tuned neural network architecture and then create the final architecture:

The final architecture:

We can see the concatenate layer in the center that concatenate our 3 previous neural network architectures; image, text and structured data respectively.

Results and Conclusions

If we take a quick look at our first results:

The worst model of all is the dummy mean house price model where the price of each house = average price of the whole training set.

The second worst model of all is the Image Data NN Architecture alone where there seems to be very slight house price predictive power, meaning that house images doesn’t correlate very much with the house prices.

Does it make sense ?

The price of a house is not limited to the appearance of the exterior. The local real estate market seems a lot more important. For example, if we take the same house and simply move it to another regionCluster, the price could be very different. Therefore, it should not be surprising that our Image Data NN Architecture alone trained on image didn’t perform as well as the Structured Data NN Architecture trained on the attributes containing much more context about the local real estate market. Since our result is bad, we could use some strategies to improve our image data neural network and dataset:

• Having a dataset where photos are more consistent and taken from the same angles could be helpful.
• It is difficult to represent a full house using a single image. Having a dataset with multiple photos for each house (interior, exterior, kitchen, bathroom, etc.) could have helped our model a lot. Multiple images can be done by creating a montage that combines all images into a single image and then pass the montage through the CNN input.
• Direct regression task on images maybe not the best choice. Using bounding boxes to detect objects in images like different types of garages, windows, doors, patio, swimming pool, etc. could be an interesting alternative and surely reduced some unwanted pixel noise in images but requires much more work and time to prepare datasets.
• Maybe more epochs could have helped… choosing a different model than InceptionV3… hyper-parameter tuning… etc.

Our Text Data NN Architecture alone performed better than the dummy mean house price model. There seems to be predictive information / correlations on house prices when using the text description.

However, adding the textual description with our structured data (Hybrid [Structured+Text] NN architecture) seems to have only added noises because the performance is almost the same as structured data alone. One explanation could be that the house textual description contains redundant information from our structured features.

Third place is our Final Hybrid neural network architecture where images and text seem to have added noise to the structured features.

The 2 winners (almost tied) are the Structured Data Neural Network with a MAE of 79 105$ and the Hybrid (Structured + Text) Neural Network with a MAE of 77 911$ (Sounds high but our Austin houses present in the dataset are expensives with an average price of 499 635$)

Final words

In a context where we have better informative images and better informative texts correlated to our target, it would be very interesting to experiment again with this kind of hybrid architecture and see how it goes !

Also, keep in mind that in practice, it is difficult to combine different noise levels and conflicts between modalities. Moreover, modalities have different quantitative influence on the prediction output.

The most common method is to combine high-level embeddings from the different inputs by concatenating them, just like we did. The problem with this approach is that it would give an equal importance to all the sub-networks / modalities which is highly unlikely in real-life situations. For example, using a weighted combination of networks could be an alternative.

Table of contents :

• 1- Confusion Matrix (CM)
• 2- Multiclass Confusion Matrix (MCM)
• 3- Regression Confusion Matrix
• 4- RCM Accuracy AUC Score
• 5- RCM F1 AUC Score
• 6- Conclusion

1-Confusion Matrix (CM)

Confusion matrix is a popular tool for summarizing the performance of a classification algorithm (A model which is used to predict a discrete variable) by giving a better idea of what is correctly predicted and what types of errors your algorithm is making. The confusion matrix is a N x N matrix, where N is the number of classes or outputs. For 2 classes, we obtain a 2 x 2 confusion matrix.

The name « confusion » stems from the fact that it makes it easy to see whether the system is confusing / mislabeling two classes. It have two dimensions “actual” vs “predicted” and for 2 X 2 matrix (binary problems), it produce 4 count metrics :

• True Positive (TP) : predicted(1), actual(1)
• False Positive (FP) : predicted(1), actual(0)
• True Negative (TN) : predicted(0), actual(0)
• False Negative (FN) : predicted(0), actual(1)

Those 4 count metrics are also very usefull for calculating a bunch of classification model performance metrics like precision, recall, f1 score on go on.

2-Multiclass Confusion Matrix (MCM)

The binary confusion matrix can be extended for multiclass classification problems. For example, for 3 classes we obtain a 3 X 3 confusion matrix.

How to calculate TP, FN, FP, TN? Pretty much the same logic here but we have to calculate the sum of values for each class. For example, using the 3 X 3 confusion matrix above we obtain:

Once we have our TP, FN, FP, TN for each class there are two main approaches in defining performance metrics called “micro” and “macro”. The difference between macro and micro averaging is that macro weighs each class equally whereas micro weighs each sample equally. The macro and micro will have the same score if we have an equal number of samples for each class.

3-Regression Confusion Matrix

Continuous data is measured while discrete data is counted… Is it possible to extend the use of the confusion matrix from a discrete variable problem (classification) to a continuous variable problem (regression)?

To make it possible we have to convert our continous variable to a discrete form first!

How can we do that ? Simply by discretizing the continuous values into regions (bins). Discretization is the process of transforming continuous variables into a discrete form. We do this by creating a set of continuous intervals:

Let’s do it!

Suppose we have this simple regression problem X vs Y :

1. The first step is to fit our regression model :

Code example :

2.Second, we discretize the target Y in N bins. We will start with 4 bins. We can use pandas « qcut » function to do that:

Code example :

Result :

We can see that we have almost the same sample size for each discretized bin label. Here is the new dataframe :

3. Next, we have to merge the Y True continous intervals with the Y Pred results, like this :

Result :

We are now ready to calculate our multiclass Confusion Matrix for Regression ! (Don’t forget to use the test set or validation set as we don’t want to evaluate the performance on the fitted data used for training)

LINEAR REGRESSION — CONFUSION MATRIX (4 BINS)

With this simple 4 X 4 confusion Matrix, we can see that our model is not breathtaking. A perfect model would result in a perfectly centered dark diagonal which is far from the case here !:

Perfect Model example :

With the Linear Regression metric, we have more like a reverse « L » shape than a diagonal. We have a roc auc of 0.5250 and an accuracy of 0.2877. If we look at the regression metrics we have a R2 of 0.0695 and a RMSE of 129.303.

This is nice because we can now use pretty much all the classification metrics to benchmark the performance of our regression model.

Now if we do the same with a default CatBoostRegressor :

CATBOOST REGRESSION — CONFUSION MATRIX (4 BINS)

Seems a lot better ! We begin to see a centered diagonal. We now have a roc auc of 0.64873 and an accuracy of 0.47308. If we now look at the regression metrics we have a R2 of 0.6323 and a RMSE of 81.2865.

What happen if we split our continuous target into more bins ?

In this context, more bins means a more accurate confusion matrix of our regression model.

More bins means also smaller sample by class which result in a less credible evaluation. In general, the more bins we add, the more our model should tend to lose performance when using the confusion matrix, since we are looking at the data much more precisely.

Let’s now try with 10 bins and the Linear Regression model :

LINEAR REGRESSION — CONFUSION MATRIX (10 BINS)

Wow, now we can see more precisely how the model fits the continuous data and where the errors are. Our Linear regression model is not capable of predicting the first bins (< 279) and the last bins (> 447) as there are 0 counts in the predicted columns. As we expected, the accuracy have now dropped from 0.2877 to 0.09055. Note that the RMSE and R2 remains the same as those metrics are calculated with the overall results.

If we do the same with the CatBoost model :

CATBOOST REGRESSION — CONFUSION MATRIX (10 BINS)

Again we see a more centered diagonal compared to the linear regression. Accuracy decreased from 0.47308 to 0.14438, which is a smaller drop (~67%) than with linear regression (~80%) and still has better overall accuracy.

The nice thing about using the confusion matrix for regression problems is that it gives you a powerful visual tool to interpret the overall fit of the model and help spot problems in specific regions of the continuous data, even for multivariate problems where it is not possible to plot X vs Y (you can also see my previous article visualization trick for multivariate regression problems for more tools).

We can add a last model, a polynomial regression and compare the 3 models for 4, 10 and 24 bins:

POLYNOMIAL REGRESSION

CONFUSION MATRICES BENCHMARK

4-RCM Accuracy AUC Score

Can we calculate a classification metric that takes into account multiple bins to get a better performance evaluation and benchmark models together ?

We will experiment someting here… Let’s create a new regression-classification metrics that we can call RCM Accuracy AUC Score (RegressionConfusionMatrix Accuracy AreaUnderCurve Score)

First, let’s calculate accuracy (can be any classification metrics) using multiple bins, for each model.

Then we can do a simple interpolation and draw a curve where Y = Performance Score and X = Number of bins, like this :

Now we can calculate the area under the curve, for each model using the sklearn “auc” function that calculates area under curve using the trapezoidal rule :

Based on this new RCM Accuracy AUC Score, Catboost seems to have the bigger area!

Is it the same story with the F1 Score ?

5-RCM F1 AUC Score

Based on this new RCM F1 AUC Score, again, Catboost seems to have the bigger area under the curve!

6-Conclusion

Final model benchmark:

With the traditional regression model performance metrics, the polynomial regression model seems a little bit better than the CatBoost model in all metrics (RMSE, MAE, R2). However, with the classification metrics applied on discretized continuous range, the Catboost seems better than the polynomial regression… Why ? … This is an interesting experiment that need more digging!

To conclude, as we experimented, confusion matrix can also be used as a powerfull tool for regression model visualization, analysis and benchmark.

Even when the regression problem is multivariate, the regression confusion matrix can offer an interesting visualization tool to help understand how the model fits our target and where the error lies.

For more tools to help multivariate regression visualization you can also try this trick : Visualization trick for multivariate regression problems

There are billions of deep learning forecasting tutorials out there (exagerating a bit). So what’s special about this one?

• So far I have seen very few tutorials explaining in detail how to work with 3D tensors. I will explain it with some simple visuals to understand the steps of transformation from 2d matrix to 3d tensor.
• I’ve seen many deep learning tutorials showing good optimistic results in the learning process, but what about going beyond the test set, being completely blindsided in the future? (Sometimes this step is missing but can demonstrate big overfitting issues from the data preparation, especially when creating the tensor slices).
• How about extrapolating long-term multisteps into the futur ? 1 year into the futur ?
• We will also experiment a colossal collection of advanced neural networks architectures with our 3d tensors : Bidirectional LSTM, Seq2Seq, Seq2Seq with Attention, Transformers, etc.
• We will experiment time series vector representation embeddings with the awesome Time2Vec
• Not all time series are equal; some are deterministic, other more complex with multiple seasonal structures and cycles and some are completely stochastic. We will try different ones!
• We will implement a Multivariate-Multisteps-Multioutput deep learning forecasting model.
• At the end, I’ll give you access to the jupyter notebooks with all the codes in a way that it is easy for you to re-use it with your own time series.

Warning*: Like the title suggest, this is a more hands-on article, soo don’t expect to deep dive into all the concepts and theory (but I’m a nice guy soo I will give you some links on the way if your interested in learning more).

Also, in order not to write a long endless novel, I did not include any time series analysis (TSA) parts which is not a good practice at all ! We should always start with a good TSA to understand our time series before embarking on our forecasting journey.

Table of contents

• 1. Time serie#1 — Monthly easy deterministic
• 2. Feature Engineering
• 3. Converting 2d matrice to 3d tensor
• 4. RdR Score Performance Metric
• 5. Bidirectional LSTM (BiLSTM)
• 6. Seq2Seq
• 7. CNN-BiLSTM
• 8. TCN-BiLSTM
• 9. MDN-BiLSTM
• 10. Attention-BiLSTM
• 11. Seq2Seq with Attention
• 12. MultiHeadAttention-BiLSTM
• 13. Time2Vec-BiLSTM
• 14. Time2Vec with Transformer
• 15. Average Ensemble Model
• 16. Final Benchmark
• 17. Time serie#2 — Daily multi-seasonalities
• 18. Time serie#3 — Daily stochastic (Amazon Stock)
• 19. Multivariate-Multisteps-Multioutput
• 20. Link to all the CODE
• 21. Conclusion

1. Time serie#1 — Monthly easy deterministic

Here is the first time series we will use :

This is a monthly time serie. As you can see just by looking at the time serie, there is a strong deterministic trend with a strong deterministic monthly seasonality. If the trend and the seasonality remains the same over time, it should not be a big deal to find a model that performs very well on this.

Normally, as this is a very predictable and univariate time series, I would approach this problem using a simple SARIMA, HoltWinters, FbProphet or Thymeboost model, but the goal here is to experiment deep learning models with tensors !

For this first experiment, we will try to forecast 24 months into the future !

2. Feature Engineering

This diagram illustrates the feature engineering pipeline we’ll be experimenting with (the green cells are the steps used for this particular time series but the full pipeline is implemented in the notebook):

The parameters:

The corresponding pipeline:

As you can see in this schema, by splitting at the start of the feature engineering process, we ensure that we do not add any overfitting signals from the test sequence in our train sequence.

We want our deep learning model to better understand the historical patterns (autoregressive correlations) soo we will add some lagged values as features (X).

We also need to compute the leaded values that will represent the future timesteps targets our model will train on (Y).

Lagged and leaded values:

Seasonal Means:

Since there seems to be strong monthly seasonality in our time series, we will add monthly and quarterly averages as features (X).

The following function take the periodicity of the time series as input (ex: monthly = 12, daily = 365, etc.) and compute automaticly the proper seasonal averages.

We could have applied differencing to our Y time series and un-differencing after the forecast to treat the trending problem (non-stationarity of mean).

We could also have applied a log transformation on the Y time series to treat the non-stationarity of variance (seems heteroskedastic) and un-log it after the forecast.

The code is all here in the jupyter notebook but I did not included thoses processing steps as the results were better without them for this particular time series. You could definitively try this at home if you want!

Final feature engineering pipeline:

The final preprocessing step is to convert the 2d matrice into 3d Tensor. We have to be very careful here to not add any Y signals in the same X row as it would result in a highly overfitted model, very good test score and very bad / strange results in the forecasting step. A single X tensor slice should contains time series information at time T0, T-1, T-2, etc. and a single y tensor slice should contain information at time T+1, T+2, T+3, etc.

3. Converting 2d matrice to 3d tensor

How can we transform a 2d matrice into a 3d tensor ? We already have the rows (time series past observation in time), and the features (lags, means, exogenous, etc.), we now need to add the sequence length dimension.

There are other ways we can do this 3d transformation. In fact, you can try multiple 3d configurations that can lead you to very different results ! Here, I’ll show you one 3d configuration that gave me good forecasting results with the studied time series.

We will define a function that will create time slices (sequence length) from our 2d matrice where the length of the sequence match the steps we want to predict into the future :

The following function will require 4 inputs;

• 2d matrice
• column name(s) list of the target(s) we want to forecast (can manage multioutput)
• date column name
• number of steps we want to forecast.

There will also be 4 outputs;

• X_train 3dtensor
• y_train 3d tensor
• X_pred 3d tensor
• List of the final train features (metadata).

The 2d matrice to 3d tensor function:

Calling the function:

This call is for the test-validation sets. We just have to input the full time series to obtain the full training set and the forecasting set. Easy as this!

What exactly happen in this small function ? Simply, at each step, a new tensor slice is append to a list and then the list is converted to a tensor. Let’s analyze how those tensor slices are created, step by step with some simple visuals!

For example, if we want to forecast a 2 inputs, 1 output time series with 2 steps into the future, here what the function do with our 2d matrice input:

The output tensors shapes are X_train (9, 2, 2), y_Train (9, 2, 1) and X_Forecast (1,2,2).

At each iteration, a new slice was created (orange) :

The last slice:

As you can see, X and y are are offset from each other (we want our model to learn to extrapolate into the future).

Now, what if we want to forecast a 2 inputs, 1 output time series, but this time 3 steps into the future ? Well… not very different;

The shapes of our output tensors are now; X_Train (7, 3, 2), y_Train (7, 3, 1) and X_Forecast (1,3,2).

The last slice:

Ok… and maybe a tricky one! What if I have 2 features and 2 outputs to forecast 3 steps into the future ? To make this working, we can use the actual transformation function to get our 3d tensors but we will have to change our neural network model architecture a bit:

Observe that now our Y tensor is (7, 3, 2) versus a single output. This require a new neural network architecture to make this work… if you read the article until the end I will give you an example ! 😉

The shapes of our tensors are now; X_train (7, 3, 2), y_train (7, 3, 2) and X_Forecast (1,3,2).

The last slice:

Now we are ready to go to the model part!

4. RdR Score Performance Metric

To benchmark the performance of the model forecasts, we will use my custom RdR score which is deeply explained here.

In short, we can interpret the RdR score as the percentage of difference (In error and shape similarity) between your model and a simple random walk model, based on RMSE score and DTW (Dynamic Time Warping) score.

If the percentage is negative, your model is [X]% worst than randomness.

If the percentage is positive, your model is [X]% better than randomness.

In another words, the [X]% will change depending on the RMSE errors and the DTW distance around 0, which is the bound of naïve randomness.

Why RMSE ? I think that penalizing the average error with large errors reflect more the reality and the stability of the model but we could have replaced the RMSE with SMAPE or any regression metric we want.

5. Bidirectional LSTM (BiLSTM)

LSTM preserves information from inputs that has already passed through it using the hidden state. Using bidirectional LSTM will run the inputs in two ways, one from past to future and one from future to past. If you want to know more about LSTM, go here

keras implementation:

Mean RdR Score (On 12 random seeds): 0.554537

Forecast Result:

6. Seq2Seq

Seq2Seq models are LSTM with encoder-decoder architecture. If you want to know more about Seq2Seq, go here:

keras implementation:

The TimeDistributed layer is a keras wrapper that allows you to apply a same Dense (fully-connected) operation to every temporal slice, one time step at a time, of an input 3D Tensor. Example, if TimeDistributed receives data of shape (None, 100, 32, 256) then the wrapped layer (e.g. Dense) will be called for every slice of shape (None, 32, 256). Here none corresponds to samples/batch size.

The RepeatVector is used to repeat the input for set number, n of times. Here, it produce N copies of Encoder’s output. For example, if RepeatVector with argument 3 is applied to layer having input shape (batch_size, 16), then the output shape of the layer will be (batch_size, 16, 3). The RepeatVector layer acts as a bridge between the encoder and decoder modules.

Mean RdR Score (On 12 random seeds): 0.482833

Forecast Result:

The encoder-decoder did is job! The time series seems denoised.

7. CNN-BiLSTM

A CNN-BiLSTM is Convolutional layers followed by BiLSTM layers with a Dense layer on the output. The Convolutional layers are used for feature extraction and the BiLSTM layers for learning the extracted features across time steps.

keras implementation:

Mean RdR Score (On 12 random seeds): 0.585092

Forecast Result:

8. TCN-BiLSTM

In a TCN model (Temporal Convolutional Network), the convolutions are causal, meaning that there is no information leakage from future to past (Causal here simply means a filter at time step t can only see inputs that are no later than t). The network use dilated convolutions (a.k.a. atrous convolutions), a technique that expands the kernel (input) by inserting holes / spaces / defined gaps between the values. TCN can also be improved by adding residual blocks between each layer to speed up convergence and enable deeper models (At the core of residual blocks, there is a direct connection which skips some layers in between. This connection is called a « skip connection »). For this special neural network architecture we will use this library « pip install keras-tcn ». If you want more information about the TCN models, go here

keras implementation:

Mean RdR Score (On 12 random seeds): 0.483058

Forecast Result:

9. MDN-BiLSTM

A MDN-BLSTM is BLSTM layers followed by a Mixture Density Network layer that learns the parameters (mean, variance and probability) of the distributions mixture in the data. For this special neural network architecture we will use this library « pip install keras-mdn-layer ». If you want to learn more about the MDN neural network, you can read an article I wrote here :

keras implementation:

Mean RdR Score (On 12 random seeds): 0.441485

Forecast Result:

10. Attention-BiLSTM

Adding attention layer can help our BiLSTM in memorizing large sequences of data by enhancing some parts of the input data while diminishing other parts (focus). The final model should devote more focus to a smaller but more important part of the data. For this special neural network architecture we will use this library « pip install keras-self-attention » but keras also offer an Attention layer (Luong-style attention). If you want to learn more about attention, you can go here and here

keras implementation:

Mean RdR Score (On 12 random seeds): 0.536199

Forecast Result:

11. Seq2Seq with Attention

Encoder-decoder architecture where the attention alignement is applied between the decoder LSTM state of the seq2seq model and the prediction (Luong Attention Model). Is you want to learn more about seq2seq with attention go here

keras implementation:

Mean RdR Score (On 12 random seeds): 0.604704

Forecast Result:

12. MultiHeadAttention-BiLSTM

MultiHead attention runs the attention mechanism several times in parallel for consistent performance improvements over conventional attention. The independent attention outputs are then concatenated and linearly transformed into the expected dimension. For this special neural network architecture we will use this library « pip install keras-multi-head » but keras also offer a MultiHeadAttention layer as described in the paper “Attention is all you Need” (Vaswani et al., 2017). If you want to learn more about multihead attention mechanism, go here

keras implementation:

Mean RdR Score (On 12 random seeds): 0.527637

Forecast Result:

13. Time2Vec BiLSTM

Time series embeddings are a representation of time data in the form of vector embeddings. The official paper is here. According to the paper the Time2Vec representations have 3 main properties:

• 1-Capturing both periodic and non-periodic (linear) patterns with the help of periodic activation functions with learnable parameters
• 2-Invariant to time rescaling
• 3-Being simple enough so it can be combined with many models

Time2Vec Layer:

keras implementation (model):

Mean RdR Score (On 12 random seeds): 0.642626

Forecast Result:

Impressive results!

14. Time2Vec with Transformer

And finaly, the last model! Here, we use the awesome Time2Vec to create time data embeddings and we add a transformer block to it. In the Transformer, the Attention module repeats its computations multiple times in parallel. Deeper explanation available here

Transformer:

keras implementation (Model):

Mean RdR Score (On 12 random seeds): 0.532287

Forecast Result:

15. Average Ensemble Model

If we take all the results of the preceding models and make an average, we’ve got our average ensemble model!

Mean RdR Score (On 12 random seeds): 0.611903

Forecast Result:

15. Final Benchmark

Since neural networks are extremely stochastic by nature, I executed all the models with different seeds 12 times to get the final results.

Mean RdR Scores (after 12 runs):

We can see that, in this case, all the model are superior than a naïve Random Walk model.

The Time2Vec-BiLSTM is, by far, the great Winner here! What about the stability ?:

After runing 12 times the same notebook with different random seed, we can see that, given the actual hyperparameters, some model are much more stable than others (yes, 12 samples is few for establishing a representative distribution but still, it’s better than nothing considering the computing time needed to get more samples!).

Standard Deviations (all metrics):

The Attention-BiLSTM and the TCN-BiLSTM seems to be the less stable models (wider distributions).

The Average Ensemble and the Time2Vec-BiLSTM seems to be both, the more stable models (narrow distribution) and the more performant models (the highest mean RdR Score), using this experiment configuration, of course!

In forecasting, doing an Average Ensemble model is always a good idea as it can ensure higher performance stability over time than using a single model.

To improve the stability on other models we could try to tune some hyper-parameters like the batch size, the number of hidden neurons, weights decay (regularization), etc.

16-Time serie#2 –Daily multi-seasonalities

Here is the second time series we will use :

This is a daily time serie. Just by looking at the time series no obvious trend can be seen. However, the seasonal structure is much more complex than the previous one as it seems to have multiple seasonal components: Weekly, Monthly, Day of month, Week of year, etc.

For this second experiment, we will try to forecast 365 days into the future:

Because Time2Vec-BiLSTM was the best model from the last experiment, we will try it on our new time series (code remains the same).

The new RdR Score is: 0.389154

We can see that the RdR score is still a lot better than a naïve random walk model but lower than the previous time series score. In fact, this time series is more complex and should be more difficult to predict which is the case here !

Forecast Result:

If we zoom (Blue = Validation forecast, Orange = Ground truth):

Not bad at all for a first try on a 365 days forecast into the future !

17-Time serie#3 –Daily stochastic (Amazon Stock)

Here is the last time series we will use :

This is a daily time series, again. This time, no obvious seasonal structure can be seen. The trend seems to go up, exponentially-alike then stabilize and seems to go down. Very difficult to predict what will happen next just by looking at it ! We will try to forecast 120 days into the future.

If we start with a simple BiLSTM model:

The new RdR Score is: -0.676

Forecast Result:

If we zoom:

• Blue: Ground Truth
• Orange: Validation forecast (test)
• Red: Validation Naïve Random Walk
• Green: Forecast

What happened here ? Our model is worst than a simple naïve random walk model !

Well, stocks are very stochastic in nature, so there are rarely good historical recurring patterns that our model can learn from. By definition, a stochastic signal is unpredictable. In this case, we can try to improve the results by adding some exogenous features (ex: other stocks, tweets, etc.) and hope that there is some temporal cross correlations between the exogenous and endogenous time series that our model will catch during the training process.

This little diagram can help to understand the different types of time series:

The same experiment seems to be better with the Time2Vec BiLSTM (RdR Score: 0.288) … however great care needs to be taken here due to the unpredictable stochastic nature of the time series itself (A good one-time validation test does not necessarily imply that the model will always going to perform well in the future):

Adding exogenous features:

We add 3 exogenous features and set “exog = True”:

• MSFT_CLOSE (Microsoft stocks)
• covid (Google Trends: covid search term)
• amazon (Google Trends: amazon search term)

BiLSTM:

The new RdR Score is: -0.034. Small performance increase, but still worst than Naïve Random Walk:

Time2Vec BiLSTM:

The new RdR Score is 0.603 vs 0.288 before. Considerable performance increase here.

Again, as the stochastic nature of the time series should trigger us an alert, we should dig deeper into the temporal cross correlations analysis between the added exogenous features and the amazon stocks as both model got an increased performance with it. We could also test with different random seeds and different time sequences to measure the stability of the performance over time.

18-Multivariate-Multisteps-Multioutput

To finish, we will use the same time series from the first experiment but we will extract the Trend, Seasonal and Residuals components and we will simultaneously forecast the 3 resulting time series using a Multivariate Multisteps Multioutput neural network model:

We now have 3 time series we want to forecast at the same time using the same model:

• EASY — Trend (Output #1)

• EASY — Seasonal (Output #2)

• EASY — Residuals (Output #3)

*WARNING*: There is something very important we need to take into consideration here. If you look at the 3 time series, the Seasonal seems stationary but the Trend have a mean non-stationary issue and the Residuals have a variance non-stationary issue (heteroskedastic variance).

If we want our model to perform on both 3 time series at the same time, we can’t just apply the same transformations (ex: differencing or log) on all the time series, we need to apply specific transformation for specific time series depending on the specific needs.

This is why, if you look at the new parameter of the configuration dictionary, stabilize_mean and stabilize_variance now takes list of columns as input instead of a binary indicator TRUE / FALSE. This means that the transformations will be applied only on the specified times series in the lists.

According to the parameters:

• Trend time series will receive a differencing treatment
• Residuals time series a log transformation
• Seasonal time series will not receive any treatment

For the 2d matrice conversion to 3d tensor, we can re-use our precious function as it take a list of y targets colnames as input. The transformation logic is the same!

• The X_Train 3d tensor shape is : (139, 24, 45)
• The y_Train 3d tensor shape is : (139, 24, 3)

If we continue with our favorite Time2Vec-BiLSTM model, here is 2 possible ways of constructing our neural network architecture to be able to manage multiple outputs:

1- Using the TimeDistributed layer at the end, like this:

2- Using the Reshape layer at the end, like this:

If we take the TimeDistributed architecture, here is the results:

• RdR Score for Trend: 0.84317
• RdR Score for Seasonal: 0.75001
• RdR Score for Residuals: 0.30810

Trend Forecast Result:

Seasonal Forecast Result:

Here, there seems to have some problems. We should go back to hyper-parameters and feature engineering tuning steps ! This time series should be very easy to understand by our model as it’s the same repeating pattern over time.

Residual Forecast Result:

Here, normally, the residual component should be pure random white noise. We may think that the initial decomposition didn’t do its job because there is some structure left in the remaining residues.

For fun, if we combine the 3 forecasts together (additive) we obtain:

That’s a very interesting experiment !

19-Link to all the code

In my public github, here, you have all the code needed to calculate easily my custom RdR Score as well as the custom “differencing” and “Log” preprocessing classes I made for time series.

I also added 4 jupyter notebooks:

• Advanced DL Forecasting — EASY.ipynb
• Advanced DL Forecasting — MEDIUM.ipynb
• Advanced DL Forecasting — HARD.ipynb
• Advanced DL Forecasting — EASY Mulivariate-Multisteps-Multioutput.ipynb

These are templates that you can easily adapt with your own time series:

Each notebook has a table of contents where you can quickly access the parts of interest:

Each model is split into 4 parts:

• Model Architecture
• Model Training
• Model Validation
• Model Forecasting

Enjoy!

20-Conclusions

In this tutorial, we experimented various deep learning architectures with a set of diverse time series using 3d tensors. Here are some on the conclusions we can retain from this experiments:

• Preparing the 3d tensor from the 2 matrice is a crucial step in the deep learning forecasting journey.
• We have to be very careful not to add any future “y signals” in the “X features” when preparing the tensors. You can start by spliting your time series before doing any preprocessing step.
• There is no “one size fit all” model or feature engineering for time series problems.
• Multivariate Multisteps Multioutput models are powerful but tricky and not always a good idea if the mixed time series to be produced are too heterogeneous (better to treat them individually with different parameters or preprocessing).
• Time2Vec is a great vector embeddings tool to try for improving model performance.
• A more complex model does not always mean a better model.
• Using an average ensemble of multiple heterogenous forecasting models can help with the stability of the results over time.
• To perform well on an autoregressive (univariate) time series forecasting problem, the time series itself must have a minimum of historical deterministic structural patterns. Generally speaking, the model will perform poorly on univariate stochastic time series.
• We can try to improve the results of stochastic time series by adding some exogenous features (temporal cross correlations)
• A good one-time validation test does not necessarily imply that the model will always perform well in the future. We can try different random seeds and different historical time sequences to measure the stability of the performance over time. We can also dig into temporal correlations to find results explanations.
• As we extrapolate over time, the results should always be used with extreme caution. As the model learns from past structure and past correlations, it cannot predict future changes it has never seen before such as data drift, context drift or new disruptive events.

In this text, I will propose you an experimental technic to evaluate the performance of time series forecasting models but before, we will quickly surf on popular time series scoring technics:

1. MAE, RMSE and AIC
2. Mean Forecast Accuracy
3. Warning: The time series model EVALUATION TRAP!
4. RdR Score Benchmark

This new RdR Score technique will give several benefits such as being able to:

• Compare models together and select the best one
• Facilitate explanation to manager or business team
• Help decide whether the forecasting model should be use or not
• How much a forecasting model is good, alone or compared to other models
• Use the shape similarity of the forecast as an important evaluation criterion
• Use randomness as an important evaluation criterion; is the forecasting model better than a naïve random decision? How much better?

The proposed RdR metric use:

• R: Naïve Random Walk
• d: Dynamic Time Warping
• R: Root Mean Squared Error

*Warning: This is very experimental and not issued from any research paper. I named this technic RdR score only to give this experimentation a name and to facilitate the understanding of what the metric actually do. Use it at your own risk!

The code of this experimental RdR score technic is available on my github.

To experiment this metric, we will work on three different datasets and four models to resolve multistep forecasting problems:

• SARIMA (Box-Jenkins method)
• Holt-Winters (Triple Exponential Smoothing)
• LightGBM (Gradient Boosting) (as a Multivariate Multi-Target Regressor)
• Seq2Seq (Deep Learning) (as a Multivariate Multi-Target Regressor)

I will assume that you know those models. If not, those models are a mix of popular old school econometric forecasting, machine learning and deep learning. There is plenty information on them if you search on Google.

1MAE, RMSE and AIC

Currently, the most popular metrics for evaluating time series forecasting models are MAE, RMSE and AIC.

If you want to read more about those metrics:

· MAE vs RMSE: https://medium.com/human-in-a-machine-world/mae-and-rmse-which-metric-is-better-e60ac3bde13d

· AIC: https://en.wikipedia.org/wiki/Akaike_information_criterion, https://stats.stackexchange.com/questions/24116/one-sentence-explanation-of-the-aic-for-non-technical-types, https://towardsdatascience.com/the-akaike-information-criterion-c20c8fd832f2

To briefly summarize, both MAE and RMSE measures the magnitude of errors in a set of predictions. The major difference between MAE and RMSE is the impact of the large errors. For example, if some prediction data points are large outliers errors when compared to the ground truth, those large errors will be diluted in the mean of MAE while RMSE score will be higher because of the square operation.

AIC measure the loss of information and penalize the complexity of a model. It is the negative log likelihood penalized for a number of parameters. The main idea behind AIC is that model with less number of parameters is better. AIC lets you test how well your model fits the dataset without overfitting it.

With MAE and RMSE, the perfect score is 0 (The goal is to have the lowest score possible). Both values range from 0 to ∞, and depend on the scale of the target we want to forecast.

MAE is the easiest to interpret. For example if MAE is 450$, we could say that our forecasting model have an average error of 450$ per forecast (+ or -). Very easy to explain to a manager or a business team. Is it good? Should we use it? Well, it depends on the use case context, the distribution of the errors (skewed or not, outliers or not) and many other things.

RMSE interpretation is less intuitive. If RMSE is 450$, we could say that our forecasting model have an “penalized” average error of 450$ per forecast, which is somewhat not very intuitive.

For AIC, there is no perfect score, the lower is the better. Soo it is not possible to evaluate the AIC score alone, it is only use to compare models together. For example, if used alone, an AIC of -950 does not give any clue about the model performance and impossible to explain to a manager.

2Mean Forecast Accuracy

The Mean Forecast Accuracy is also an interesting metric. This metric is very intuitive and easy to explain to a manager (Our model have an average of 66% forecasting accuracy, which also means that our models have an average forecasting error of 34%). It gives a good approximation of how well a forecasting model performs. For example, take these forecasting results:

Well… Nothing is perfect!

We can see that this metric have a major flaw. The mean has one main disadvantage; it is particularly susceptible to the influence of outliers. When the forecast results are really bad (when the error alone is higher than the expected ground truth result), the percentage can goes really low (In this example: (1 — (225/25)) gives Minus 800%) which will have a big negative impact on the global Mean Forecast Accuracy.

· A solution to this problem could be to clip the minimum percentage value to 0%, to reduce the impact of isolated / outlier forecast results.

· We could also have use the median instead of the mean

For example, if we used median instead of zero clipping, it gives us a Median Forecast Accuracy of 85% (the outlier have been ignored), compared to 66% with Mean Forecast Accuracy with zero clipping.

In general, when your error distribution is skewed, you should use Median instead of Mean. In some case, the Mean Forecast Accuracy could also mean totally nothing. If you remember your statistics; the coefficient of variation (CV) represents the ratio of the standard deviation to the mean (Coefficient of Variation = (Standard Deviation / Mean) * 100). A big CV value means big variability, which also means greater level of dispersion around the mean. For example, we could consider anything above a CV of 0.7 as highly variable and not really forecastable. In addition, it can show that your forecasting model prediction’s ability is very unstable! (The mean have a poor meaning in this case) (source: https://blog.arkieva.com/do-you-use-coefficient-of-variation-to-determine-forecastability/)

3Warning: The time series model EVALUATION TRAP!

Okay, before we deep dive into the experimental RdR technic, I would like to talk about one important mistake I have seen so many times on the internet: Please, do not calculate performance metrics on the fitted data / training data!

In traditional statistics, time series forecasting models were often evaluated on the “fit” (“find the best fit”) predictions results of the model, which I think, makes absolutely no sense at all! In machine learning process, are we evaluating / choosing our models on the fit results? NEVER (I hope so!). Unless you want to know if your model is overfitting or underfitting (you can compare validation and training set results)

I think this way of thinking comes from old fashion curve fitting technics where the goal is to fit a curve as best as possible. When you add the time dimension, the problem is that when we try to extrapolate this overfitted curve through time, the results are rarely good:

Another big problem with that way of thinking is, when you use advanced technics like gradient boosting ensemble or deep learning, the fit will usually be very good or perfect and hide overfitting problems that you will only discover in the extrapolation process through time.

I have read several experiments where a deep learning model is perfectly fit on a time series with the conclusion: “Wow, deep learning with time series is so revolutionary and incredible!” (And of course, you don’t see the extrapolation process through time that should come next and that probably didn’t work very well…). Fitting is not a problem anymore (You just have to set a big number of estimators to an ensemble model or a big number of epochs to a neural network model), overfitting is!

4RdR Score Benchmark

The proposed RdR score technic will mainly answer three questions:

How can we take into account the shape similarity of a time series? :

• Answer: Dynamic Time Warping

How can we know if we should use our forecasting model or not? :

• Answer: Is it better or worst than a Naïve Random Walk?

How can we take into account the errors? :

• Answer: Root Mean Squared Error (RMSE)

Now, suppose we have those two following models, with exactly the same RMSE:

Which model would you choose?

One big flaw about using MAE, RMSE or AIC; those metrics does not take into account an important criterion of the forecast: THE SHAPE SIMILARITY!

Why use Dynamic Time Warping (DTW) as a similarity metric? :

• Euclidean distance between time series: Bad choice because there is distortion in the time axis

• DTW: Find the optimal (minimum distance) warping path between two time series, by « synchronizing » / « align » the different signals on the time axis

• The DTW distance is the square root cost of optimal warping path

• The lowest the distance on the warping path, the more similar are the time series.

If you want to know more about Dynamic Time Warping (DTW): https://www.slideshare.net/DavideNardone/accelerating-dynamic-time-warping-subsequence-search-with-gpu

If a manager ask you, “should we use our forecasting model to predict the future? Is it good or not?” An interesting answer could be: Well, based on the performance evaluation metrics, our model is 65% better than if we randomly make a decision. This is where random walk shows off.

We can interpret the RdR score as the percentage of difference (In error and shape similarity) between your model and a simple random walk model, based on RMSE score and DTW score. If the percentage is negative, your model is [X]% worst than randomness. If the percentage is positive, your model is [X]% better than randomness. In another words, the [X]% will change depending on the RMSE errors and the DTW distance around 0, which is the bound of naïve randomness. Why RMSE instead of MAE? I think that penalizing the average error with large errors reflect more the reality and the stability of the model.

-∞% means that sometimes, you may have a very big negative percentage score like -98695%; This is bad news and it means that you should definitely not use the model!

There could also be a strange situation where the perfect score equals the random walk model (A straight line). Well, in his case, you do not really need a model!

Enough talking, let’s try it and see what happens!

Experiment Dataset#1: Easy (Good autocorrelation and seasonal structure, deterministic system):

The multistep unseen test data we will use to validate the performance of our models:

The naïve random walk model (RdR Score = 0):

The forecasts from a random walk model (without drift) are equal to the last observation; Scenario where future movements are unpredictable and are equally likely to be up or down (Stochastic). It gives a straight line because if we do infinite random simulations from the last data point, as the chances are equally likely to be up or down, the mean will be equal to the last observation. Like this:

Can we beat this?

Here, the multistep forecasts of 12 periods in the future (unseen test data) of our four models:

Let’s calculate the RdR score of each model:

With this graph, we can see that all the performance are very near. The best model is Seq2Seq.

If we want to have more detailed information on the RdR score, we can plot the RMSE vs DTW like this:

The y axis is for the penalized errors while the x axis is for the shape similarity of the time series (between prediction and unseen test dataset). In this graph, we can see that Seq2Seq was the best model, in both error and shape and that it is halfway between the random walk score and the perfect score (which represent the 50,32% RdR score). We can see that LightGBM have slightly more errors than the econometric models (RMSE) but have a slightly better shape similarity (DTW distance). All models were a lot better than the Random Walk model performance.

Interpretation of the best model (Seq2Seq):

If we zoom the prediction of the best model (Seq2Seq):

Experiment Dataset#2: Medium (Autocorrelation and seasonal structure, with a lot of noise):

The multistep unseen test data we will use to validate the performance of our models:

The naïve random walk model (RdR Score = 0):

Can we beat this? :

Here, the multistep forecasts of 12 periods in the future (unseen test data) of our four models:

Let’s calculate the RdR score of each model:

We can see that the SEQ2SEQ model was the best. SARIMA performed the worst but all models are still better than the random walk. As we expected, we can also see that globally, this round was more difficult than the previous one; the RdR scores are lower.

If we want to have more detailed information on the RdR score, we can plot the RMSE vs DTW like this:

The y axis is for the penalized errors while the x axis is for the shape similarity of the time series. In this graph, we can see that Seq2Seq was the best model, in both error and shape. SARIMA was almost equal RMSE than random walk but the shape of SARIMA time series was better (DTW distance).

Interpretation of the best model (Seq2Seq):

If we zoom the prediction of the best model (Seq2Seq):

Experiment Dataset#3: Hard (Poor autocorrelation, poor seasonal structure, stock price alike):

The multistep unseen test data we will use to validate the performance of our models:

The naïve random walk model (RdR Score = 0):

Can we beat this?

Here, the multistep forecasts of 12 periods in the future (unseen test data) of our four models:

Let’s calculate the RdR score of each model:

Just by looking at this bar chart, we should not use both Holt-Winters and LightGBM as they are worst (in error and shape similarity) than a simple naïve random walk model.

If we want to have more detailed information on the RdR score, we can plot the RMSE vs DTW distance like this:

The y axis is for the penalized errors while the x axis is for the shape similarity of the time series. In this, we can see that Seq2Seq was the best model while Holt-Winters and SARIMA were very near the random walk performance. LightGBM was out, the worst model in this case.

If we zoom the prediction of the best model (Seq2Seq):

Interpretation of the best model (Seq2Seq):

With RdR score, we can also try to get an overall score on the three datasets together:

Seq2Seq win this round with mean RdR score of 38,08%, Holt-Winters in second position, SARIMA then LightGBM.

Of course, we could have tune the hyper-parameters, make more data transformations, tune neural network architecture, add exogenous data, etc. but the goal here was not to put the focus on the models.

As usual, the python code of this experimental “RdR score” is avalaible on my github as jupyter notebook experiment: github.

I hope that this experiment have been useful to you!

Dave Cote, Data Scientist, M.Sc. BI

In a univariate regression problem, it is easy to see the model decision function / how your model is generalizing the problem by doing a simple X vs Y plot:

With this plot, we can easily see if the model found a good generalization solution, if it overfit, underfit, etc. It is also a nice tool to compare different models together:

When you add 1 feature to the regression problem, we can still plot the model decision function but it will be in 3D, and the generalization solution is not a line anymore but become a 3D surface:

It is more difficult to evaluate the model decision function in 3D but still, it is possible !

Now, what happen if we add other features to the regression model, how can we visualize the model regression function in 4D, 5D, 6D… nD ??? Sometimes, it is difficult to evaluate a model performance only by using the RMSE, MAE or R2 scores because it does not give any clues about how the model is reacting when the data changes, so wee need a visualization to help us out!

Off course, the visualization trick that I will present you is not perfect, nothing is! Still, it helps to understand how the model reacts when the target change and it also help to validate the overall model generalization performance. Something nice about this trick: It is always the same, no matter how many features you have in your model!

First, when you got your train fit and test predictions ready, create 2 new datasets with 2 columns each:

Dataset_1: fit_results_df:

1. Y_TRUE: Y train set original values

2. Y_FIT: Y train fit values

Dataset_2: test_pred_results_df:

1. Y_TRUE: Y test set original values

2. Y_PRED: Y test prediction values

Then, for both datasets, simply sort all the rows by Y_TRUE values in an ascending order and finally, reset the indexes.

If we plot the Y_TRUE values of both datasets, it will show:

This is all the ordered y values of the train and test set in ascending order.

Then, we plot the corresponding Y_FIT and Y_PRED values of both datasets:

This model is a simple Decision Tree with a max depth limited to three. With this visualization trick, we can see that:

· The “simplified” 2d model decision function is “step-alike” (Give a clue that this is a Decision Tree!)

· The model is probably underfitting (Could have more steps!). George St-Pierre would probably say: “I am not impressed by your performance”

· The performance seems very similar between the train fit and the predictions test set

· The model seems to follow the target variations; when the y values go up, the prediction and fit goes up too most of the time ! (Sometimes, it happens that the model is only good with high or low target values and do not follow the target variations correctly… the visualization help us to see that very easily!)

· There are patterns, bunches of fit/predictions that are out! (Maybe outliers? Need some investigations…):

Let us improve the model a bit and see what will happen! :

This new model is a Random Forest with a max depth limited to 5 and 750 estimators.

We can easily say that this model is better than the old one, without checking at the Scores! (Although, you should check the scores!). We can also see that there are still bunches of fit/predictions data points that are out in the right-bottom corner of the plot, those data points need an investigation!

If we want, we can also compare multiple models at a time (Let say we are only interested in the predictions of the test set):

If you prefer, you can also use a line instead of data points to represent the “2D simplified” model decision function on Y Target:

You can plot the train and the test sets at the same time if you want!

Example: Random Forest on boston dataset, test and train datasets

You can also add a prediction interval, like the inductive conformal prediction, around Y prediction target if needed:

Example: Deep Neural Network with an Inductive Conformal Prediction

Conclusion:

We saw that this tool:

• Is not perfect but can help to validate and visualize a “2D simplified” multivariate regression model decision function when you have more than two features in your problem
• Can be used to compare different models together
• Can help to identify patterns, for example: If the model follows correctly the target variations or if there are data points clusters that did not fit/predict correctly / need investigations (Maybe some feature engineering / cleaning steps to help the model with those data points!)
• Can be use to compare Train vs Test performance
• Can be customized depending on the preferences or needs of the data sicence practitioner (Scatter vs line plot, inductive conformal prediction, train vs test, etc.)
• Off course, this tool is not sufficient; you should continue to use regression scores like RMSE to validate your models!

I hope that this shared visualization trick could be useful to you!

The simplified codes of this experimentation is available on my github (Jupyter notebook experimentation on the boston dataset): github

Dave Cote, Data Scientist, M.Sc. BI

In this article, I will first explain briefly what a MDN is and then give you the python code to make your own MDN model with only a few lines of code.

Before we start…

If you prefer to follow the jupyter notebook: https://github.com/CoteDave/blog/blob/master/Made%20easy/MDN%20regression/Made%20easy%20-%20MDN%20regression.ipynb

If you want to play around and experiment yourself with my experimental custom MDN python class model : https://github.com/CoteDave/blog/blob/master/Made%20easy/MDN%20regression/mdn_model.py

First… Mixture Density Network regression … huh What !??? Let’s start slow and explain separately and quickly all the main terms here:

1. Regression

2. Density

3. Mixture Density

4. Network

5. Regression Network

6. Density Network

7. Mixture Density Network regression

8. Practice fun — Univariate dataset with custom MDN class

9. Practice fun — Multivariate dataset with custom MDN class

10. Conclusion

1REGRESSION

First, the term “regression”. If you have click on this article, I assume you know what regression is all about, but if not, here a really quick and simple explanation from Jason Brownlee of Machine Learning Mastery (Great guy by the way!):

“Regression predictive modeling is the task of approximating a mapping function (f) from input variables (X) to a continuous output variable (y) […] A regression problem requires the prediction of a quantity. A problem with multiple input variables is often called a multivariate regression problem […] For example, a house may be predicted to sell for a specific dollar value, perhaps in the range of $100,000 to $200,000 […] ” (Source: https://machinelearningmastery.com/classification-versus-regression-in-machine-learning/#:~:text=A%20regression%20problem%20requires%20the,called%20a%20multivariate%20regression%20problem.)

Here another visual explanation that differentiate a classification problem from a regression problem:

A Last example for the luck…

2DENSITY

Okay, now what the term “density” is all about? Here a quick “cheesy” example:

Suppose you are delivering pizza for pizza-hut. You decide to capture the time (in minutes) of each delivery you just made. After 1000 deliveries, you decide to visualize your data to see how well you performed. Here is the result :

Here is the “density” of the distribution of your pizza delivery times data. Soo, in average, it took you 30 minutes per delivery (The peak in the graph). It also says that 95% of the time (2 standard deviations (SD)), your delivery took between 20 and 40 minutes to accomplish. Here, the density kind of represents the “frequency” of the time results. The difference between “frequency” and “density” is that:

· Frequency: If you draw a histogram under this curve and counts all the bins, it will sum up to any whole number (depending on the total number of observations captured in your dataset).

· Density: If you draw a histogram under this curve and counts all the bins it will always sum up to 1. We can also call this curve the probability density function (pdf).

In statistical terms, this is a beautiful normal / Gaussian distribution. This normal distribution have two parameters:

· The mean

· The standard deviation (According to Wikipedia: “Standard deviation is a number used to tell how measurements for a group are spread out from the average (mean), or expected value. A low standard deviation means that most of the numbers are close to the average. A high standard deviation means that the numbers are more spread out.”)

The change of mean and standard deviation will affect the shape of the distribution. For Example:

Of course, there exists a variety of different distributions types with different kind of parameters. For example:

3MIXTURE DENSITY

Okay now let’s have a look at those 3 distributions:

If we take this bimodal distribution (Also called a general distribution):

Mixture Density Network use the assumption that any general distribution like this bimodal distribution can be broken down into a mixture of normal distributions (The mixture can also be customize with other type of distributions like Laplace for example):

Soo if you add up those 2 normal distributions (red and blue), where each distribution have is own mean and standard deviation, it will gives back the bimodal general distribution (purple).

4NETWORK

Network for Artificial neural network. Soo yes, a Mixture Density Network is a type of artificial neural network. Here a classical example of a neural network:

You have the input layer (yellow), the hidden layer (green) and the output layer (red). (I assume here that everyone know what an artificial neural network is, if not, go on google to find out!).

5REGRESSION NETWORK

A regression network is an artificial neural network where the goal is to learn to output a continuous value given some input features. In the example above, given the age, sex, education and other features, the artificial neural network learns to predict the salary.

6DENSITY NETWORK

A density network is an artificial neural network where the goal is not to simply learn to output a single continuous value, but to learn to output the distribution parameters (here, mean and standard deviation), given some input features. In the example above, given the age, sex, education and other features, the artificial neural network learns to predict the mean and the standard deviation of the expected salary distribution. Predicting the distribution versus a single value have advantages like being able to gives uncertainty boundaries with the prediction. This is a “Bayesian” approach of resolving a regression problem. Here is a great illustration of predicting the distribution of each expected continuous values:

Another great illustration that shows us the distribution of the expected values, for each predicted instance:

7MIXTURE DENSITY NETWORK

Finally! A mixture density network is an artificial neural network where the goal is to learn to output all the parameters (here, the mean, standard deviation and Pi) of all the distribitions mixed in the general distribution, given the specific input features. The new parameter “Pi” is the mixture parameter that gives the weights/probability of a given distribution in the final mixture.

The final result:

8PRACTICE FUN — UNIVARIATE DATASET WITH CUSTOM MDN CLASS

Okay, enough of this theoretical stuff, let’s do it!

Here we have this famous “half-moon” dataset:

Just by looking at the data, we can see that there is two overlapping clusters.

If we plot the density distribution of the target value (y):

We have a nice “ghost-like” multimodal distribution (general distribution). If we try a standard Linear Regression on this dataset to predict y with X:

Did not work very well!

Let’s now try a nonlinear model (Radial Basis Function Kernel Ridge regression):

Better, but not quite yet!

The main reason why both model did not make it is because, if you look on the X-axis, multiple different y values exists for the same X value… more specifically, there seems to exists more than one possible y distribution for the same X. The regression model just tried to find the optimal function that minimize the error and did not took in consideration a mixture of densities, which is the case! Those X in the middle do not have a unique Y solution, they have two possible solutions, one high and one low!

Let’s now try a MDN model to see what this beast can do! For that, I implemented for you a quick nice and easy to use “fit-predict”, “sklearn alike” customized python MDN class. Here is the link to my python code if you want to use it yourself (Be aware: this MDN class is experimental and has not been extensively tested): https://github.com/CoteDave/blog/blob/master/Made%20easy/MDN%20regression/mdn_model.py

To be able to use this class, you will need to have sklearn, tensorflow probability, Tensorflow < 2 (If you want to make this MDN model class compatible with TF2, feel free to contribute!), umap and hdbscan (for custom visualization class function).

· n_mixtures: number of distribution mixtures to use by the MDN. If set to -1, it will “auto” find the optimal number of mixtures using a gaussian mixtures model (GMM) and a HDBSCAN model on X and y.

· dist: Distribution type to use in the mixture. For now, there are two choices available; “normal” or “laplace”. (Based on some experimentations, Laplace distribution gave me better results than normal distribution).

· input_neurons: Number of neurons to use in the input layer of the MDN

· hidden_neurons: Architecture of the hidden layers. Of the MDN. List of neurons per hidden layer. This parameter give you the ability to choose the number of hidden layers and the number of neurons per hidden layer.

· gmm_boost: Boolean. If set to True, will add cluster features to the dataset.

· optimizer: Optimization algorithm to use.

· learning_rate: Learning rate of the optimization algorithm

· early_stopping: To avoid overfitting when training. This trigger will decide when to stop training when no change in metric over a given number of epochs.

· tf_mixture_family: Boolean. If set to True, will use the tf_mixture family (Recommended): The Mixture object implements batched mixture distributions. The mixture model is defined by a Categorical distribution (the mixture) and a python list of Distribution objects.

· input_activation: Activation function of the input layer

· hidden_activation: Activation function of the hidden layer

Now that our MDN model is fitted on the data, let’s samples from the mixture density distribution and plot the probability density function:

Wow! Our MDN model fitted pretty well the true general distribution!

We can broke down the final mixture distribution into each distribution to see how it looks:

Let’s sample again some Y data using our learned mixture distribution and this time, plot the generated samples against the true ones:

Again, pretty close from reality!

Given X, we can also generate multiple batch of sample to produce statistics like quantile, mean, etc.:

Now, we can plot the mean of each of our learned distributions, with their respective mixture weights (pi):

Wow, MDN just nailed it, again!

As we have the mean and the standard deviation from each distribution, we can also plot the uncertainties with the full tail; let’s say we plot the mean with a confidence interval of 95%:

We can also mix the distribution together and, when we have multiple y distributions for the same X, we choose the most probable mixture using the highest Pi parameter value:

· Y_preds = for each X, choose the Y mean of the distribution with the max probability/weight (Pi parameter)

If we add the 95% confidence interval:

Not ideal here, because there are clearly two distinct clusters overlapping in our data with near equivalent density. The errors would be higher than a standard regression model. This also means that maybe there is an important feature missing in our dataset that could help to avoid the cluster overlapping in a higher dimensionality.

We can also choose to mix the distribution using the Pi parameters and the means of all the distribution together:

· Y_preds = (mean_1 * Pi1) + (mean_2 * Pi2)

Which gives:

If we add the 95 confidence interval:

This option gives us pretty much the same result as a non-linear regression model, mixing everything to minimize the distance between the dots and the function. This gives a single function that pass through the two clusters.

In this very particular case, my favorite option would be to assume that in some region of the data, X have more than one Y while in other regions; only one of the mixture is used. That give use something like:

For example, when X = 0, there is two possible distinctive Y solutions from each mixture. When X = -1.5, there is a unique Y solution from mixture 1. Depending on the use case or business context, an action or decision can be triggered when more than one solution exists for the same X.

With this option, the rows are duplicated when there is an overlapping distribution (If both mixture probabilities are >= given probability threshold). This give something like:

With the 95% confidence interval:

Dataset rows passed from 2500 to 4063 and the final predictions dataset looks like:

In this data table, we see that, for example, when X = -0.276839, Y can be 1.43926 (with 0.351525 of probabilities of mixture_0) but can also be -0.840593 (with 0.648475 of probabilities of mixture_1).

The instances with multiple distributions gives also an important information that something is going on in the data and that maybe it needs more analysis. Maybe it is some data quality problems or maybe it indicates that an important feature is missing in the dataset!

Another example: “traffic scene prediction is a good candidate for where Mixture Density Networks can be used. In traffic scene prediction, we need a distribution over behaviors an agent can exhibit — for example, an agent could turn left, turn right, or go straight. Thus, Mixture Density Networks can be used to represent “behaviors” in each of the mixture it learns, where the behavior consists of a probability and trajectory ((x, y) coordinates up to some time horizon in the future).(source: http://blog.adeel.io/tag/mixture-density-networks/)“

9PRACTICE FUN — MULTIVARIATE DATASET WITH CUSTOM MDN CLASS

Okay finally, does MDN perform well with multivariate regression problems? Let’s find out !

First of all, let’s have a look at the dataset we will use:

· age: age of primary beneficiary

· sex: insurance contractor gender, female, male

· bmi: Body mass index, providing an understanding of body, weights that are relatively high or low relative to height, objective index of body weight (kg / m ^ 2) using the ratio of height to weight, ideally 18.5 to 24.9

· children: Number of children covered by health insurance / Number of dependents

· smoker: Smoking

· region: the beneficiary’s residential area in the US, northeast, southeast, southwest, northwest.

· charges: Individual medical costs billed by health insurance. This is the target that we want to forecast

I took this data dictionary from https://www.kaggle.com/mirichoi0218/insurance. If you go to this link, you can also have an interactive visualization of the dataset.

If you want to do some expriments, the dataset is also available here: https://github.com/stedy/Machine-Learning-with-R-datasets

The problem statement is: Can you accurately predict insurance costs (charges)?

To see a quick exploratory data exploration on this dataset, you can read one of my previous article here: https://www.linkedin.com/pulse/demonstrating-power-feature-engineering-part-ii-how-i-dave-c%C3%B4t%C3%A9/

Now, let’s import the dataset:

…and do some basic data preparation!:

The data is ready to feed the MDN model!

Let’s now predict the test dataset using the “best mixture probabilities (Pi parameter) strategy” and plot the result (y_pred vs y_test):

For this visualization technic, you can read one of my previous article: https://www.linkedin.com/pulse/visualization-trick-multivariate-regression-problems-dave-c%C3%B4t%C3%A9/

Wow! This is completely insane! With a R2 of 89.09 and a MAE of 882.54, The MDN just destroyed it!

Let’s plot the fitted vs true distribution for fun:

Almost exactly the same!

If we broke down the mixture model:

There is a total of six different distributions in the general mixture distribution.

We can generate multivariate samples from the fitted mixture model (a PCA is applied to visualize the result in 2D):

The generated samples are very close to the real ones!

We can predict from each distribution if we want:

We can check if there are overlapping cluster in the dataset:

y_pred_overlaps have the same number of rows than the X_test dataset, so there is no overlapping! (We can adjust the pi parameter threshold for the sensibility)

10CONCLUSION

· MDN did a great job in the univariate regression dataset where two clusters overlapped on each other and where X could have more than one Y output, compared to linear or nonlinear classical ML models.

· MDN did also a great job in the multivariate regression problem and can compete with popular model like XGBoost

· MDN is a great and unique tool to have in your ML toolbox and can resolve particular problems that other models can’t do (capable of learning from data obtained from a mixed distribution)

· As MDN learn the distribution, you can also compute uncertainty with your predictions or generate new samples from the learned distribution

I hope this article helped you in your Machine Learning journey!