Prediction Interval for Neural Net in Python

Question

I'm currently using keras to create a neural net in python. I have a basic model and the code looks like this:

    from keras.layers import Dense
    from keras.models import Sequential

    model = Sequential()
    model.add(Dense(23, input_dim=23, kernel_initializer='normal', activation='relu'))
    model.add(Dense(500, kernel_initializer='normal', activation='relu'))
    model.add(Dense(1, kernel_initializer='normal', activation="relu"))
    model.compile(loss='mean_squared_error', optimizer='adam')

It works well and gives me good predictions for my use case. However, I would like to be able to use a Variational Gaussian Process layer to give me an estimate for the prediction interval as well. I'm new to this type of layer and am struggling a bit to implement it. The tensorflow documentation on it can be found here:

https://www.tensorflow.org/probability/api_docs/python/tfp/layers/VariationalGaussianProcess

However, I'm not seeing that same layer in the keras library. For further reference, I'm trying to do something similar to what was done in this article:

https://blog.tensorflow.org/2019/03/regression-with-probabilistic-layers-in.html

There seems to be a bit more complexity when you have 23 inputs vs one that I'm not understanding. I'm also open to other methods to achieving the target objective. Any examples on how to do this or insights on other approaches would be greatly appreciated!

Answer 1

tensorflow_probability is a separate library but suitable to use with Keras and TensorFlow. You can add those custom layers in your code and change it to a probabilistic model. If your goal is just to get a prediction interval it would be simpler to use the DistributionLambda layer. So your code would be as follows:

from keras.layers import Dense
from keras.models import Sequential
from sklearn.datasets import make_regression
import tensorflow_probability as tfp
import tensorflow as tf 
tfd = tfp.distributions

# Sample data 
X, y = make_regression(n_samples=100, n_features=23, noise=4.0, bias=15)

# loss function Negative log likelyhood
negloglik = lambda y, p_y: -p_y.log_prob(y)

# Model
model = Sequential()
model.add(Dense(23, input_dim=23, kernel_initializer='normal', activation='relu'))
model.add(Dense(500, kernel_initializer='normal', activation='relu'))
model.add(Dense(2))
model.add(tfp.layers.DistributionLambda(
      lambda t: tfd.Normal(loc=t[..., :1],
                           scale=1e-3 + tf.math.softplus(0.05 * t[..., 1:]))))

model.compile(loss=negloglik, optimizer='adam')

model.fit(X,y, epochs=250, verbose=None)

After training your model, you can get your prediction distribution with the following lines:

yhat = model(X) # make predictions
means = yhat.mean() # prediction means 
stds = yhat.stddev() # prediction standard deviation