Prediction Interval for Neural Net in Python

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!

from Prediction Interval for Neural Net in Python