<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Bayesian Long Short-Term Memory Network Example
Licensed under the Apache License, Version 2.0.
This example implements a Bayesian version of LSTM (Hochreiter, Schmidhuber, 1997) using tf.keras and edward2.
It's based on the discussion in this issue.
Step1:  Data generation
We generate random, valid and invalid strings from the embedded reber grammar as inputs. The objective is a binary classification between valid and invalid strings.
Step2:  Hyperparameters
Step3:  Training a baseline LSTM
This is our baseline model. It's a simple RNN with a single layer of LSTM cells.
As it would be expected the model can quickly be fitted to differentiate the valid from the invalid strings.
Step4:  Training a Bayesian LSTM
In comparison we build a similar RNN we we replace the standard LSTM cells with the LSTMCellFlipout class from edward2.
Step5:  From the training curves and the test results it's clear that the model is not converging and we are not able to differentiate the valid from the invalid strings.
Step6:  Hyperparameter tuning
In this section we run the experiments over a grid of, for BNN important, parameters.
Step7:  Analysing the impact of hyperparameters
The following plot shows the test accuracy over the hyperparameters. It's clear that scaling the loss of the Bayesian layers in the model is necessary to reach convergence.
To plot those results NaN values in the clipvalues column are set to -1. In the corresponding runs gradients were not clipped.
What is the correct scaling factor for the regularizers?
BNN use the Kullback-Leibler Divergence (KL-divergence) as a regularization term to minimize the difference between the prior distribution $p$ and the approximated posterior distribution $q$ over the parameters of the network. This divergence term is written as $KL(q||p)$.  
The KL-divergence should be included once for an entire dataset. Due to the implementation  in tensorflow the regularization term is added to the value of the loss function after every batch though, leading to a larger KL-divergence term compared to the loss. To bring both terms on the same scale, we apply a scaling factor < 1 to the KL-divergence. The loss for each batch is divided by the batch size and therefore $KL(q||p)$ (which is in relation to the complete dataset) is divided by the number of training samples (in this case x_train.shape[0]).
Step8:  Choosing the correct batch size as well as possibly a value to clip gradients by during the training is more complicated.
From those experiments it looks like a smaller batch size will be benefitial, even though this is also dependent on the dataset. Running this experiment with less samples shows that a larger batch size can greatly increase accuracy. A possible explaination for this behaviour is that in cases where less data is available, larger batch sizes compensate for the additional variance that is caused by the Bayesian properties.
Using this dataset, clipping the gradients doesn't seem to have a huge affect on the outcome although it is often an important technique with Bayesian LSTMs.
Step9:  Training the best Bayesian Model
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import seaborn as sns
import pandas as pd
import tensorflow as tf
import edward2 as ed
import matplotlib.pyplot as plt
from tqdm import tqdm
from sklearn.model_selection import train_test_split, ParameterGrid
from tensorflow.keras.preprocessing import sequence
import embedded_reber_grammar as erg
def plot_hist(hist):
    plt.figure(figsize=(12, 4))
    plt.subplot(1, 2, 1)
    plt.plot(hist.history['val_accuracy'], label='val_accuracy')
    plt.plot(hist.history['accuracy'], label='train_accuracy')
    plt.legend()
    plt.subplot(1, 2, 2)
    plt.plot(hist.history['val_loss'], label='val_loss')
    plt.plot(hist.history['loss'], label='train_loss')
    plt.legend()
    plt.show()
Explanation: Bayesian Long Short-Term Memory Network Example
Licensed under the Apache License, Version 2.0.
This example implements a Bayesian version of LSTM (Hochreiter, Schmidhuber, 1997) using tf.keras and edward2.
It's based on the discussion in this issue.
End of explanation
x, y = [], []
n = 3000
for i in range(n):
    x.append(np.asarray(erg.encode_string(erg.generate_valid_string(erg.embedded_rg))))
    y.append(1)
for i in range(n):
    x.append(np.asarray(erg.encode_string(erg.generate_invalid_string(erg.embedded_rg))))
    y.append(0)
x = sequence.pad_sequences(x)
x_train, x_test, y_train, y_test = train_test_split(np.asarray(x), np.asarray(y))
print(f"Number of training samples: {x_train.shape[0]}")
print(f"Number of test samples: {x_test.shape[0]} \n")
sequence_length = x_train.shape[1]
num_chars = x_train.shape[2]
print(f"Length of sequences: {sequence_length}")
print(f"Number of characters: {num_chars}")
Explanation: Data generation
We generate random, valid and invalid strings from the embedded reber grammar as inputs. The objective is a binary classification between valid and invalid strings.
End of explanation
batch_size = 64
epochs = 40
Explanation: Hyperparameters
End of explanation
model = tf.keras.Sequential()
model.add(tf.keras.layers.InputLayer(input_shape=(sequence_length, num_chars)))
model.add(tf.keras.layers.RNN(
    tf.keras.layers.LSTMCell(128)
))
model.add(tf.keras.layers.Dense(32))
model.add(tf.keras.layers.Dense(1, activation='sigmoid'))
model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
model.summary()
hist = model.fit(x_train, y_train, validation_split=0.2, epochs=epochs, batch_size=batch_size, verbose=0)
test_results = model.evaluate(x_test, y_test)
print(f"Test loss: {test_results[0]}")
print(f"Test accuracy: {test_results[1]}")
plot_hist(hist)
Explanation: Training a baseline LSTM
This is our baseline model. It's a simple RNN with a single layer of LSTM cells.
As it would be expected the model can quickly be fitted to differentiate the valid from the invalid strings.
End of explanation
model = tf.keras.Sequential()
model.add(tf.keras.layers.InputLayer(input_shape=(sequence_length, num_chars)))
model.add(tf.keras.layers.RNN(
    ed.layers.LSTMCellFlipout(128)
    ))
model.add(tf.keras.layers.Dense(32))
model.add(tf.keras.layers.Dense(1, activation='sigmoid'))
model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
model.summary()
hist = model.fit(x_train, y_train, validation_split=0.2, epochs=epochs, batch_size=256, verbose=0)
Explanation: Training a Bayesian LSTM
In comparison we build a similar RNN we we replace the standard LSTM cells with the LSTMCellFlipout class from edward2.
End of explanation
test_results = model.evaluate(x_test, y_test)
print(f"Test loss: {test_results[0]}")
print(f"Test accuracy: {test_results[1]}")
plot_hist(hist)
Explanation: From the training curves and the test results it's clear that the model is not converging and we are not able to differentiate the valid from the invalid strings.
End of explanation
params = {
    'loss_scaling': [1., 1./x_train.shape[0]],
    'batch_size': [64, 128, 256],
    'clipvalue': [None, 0.1, 0.5],
}
param_grid = ParameterGrid(params)
results = pd.DataFrame(columns=list(params.keys())+['test_loss', 'test_accuracy'])
def training_run(param_set):
    sf = param_set['loss_scaling']
    bs = int(param_set['batch_size'])
    cv = param_set['clipvalue']
    model = tf.keras.Sequential()
    model.add(tf.keras.layers.InputLayer(input_shape=(sequence_length, num_chars)))
    model.add(tf.keras.layers.RNN(
        ed.layers.LSTMCellFlipout(
            128, 
            kernel_regularizer=ed.regularizers.NormalKLDivergence(scale_factor=sf),
            recurrent_regularizer=ed.regularizers.NormalKLDivergence(scale_factor=sf)
            ),
        ))
    model.add(tf.keras.layers.Dense(32))
    model.add(tf.keras.layers.Dense(1, activation='sigmoid'))
    optimizer = tf.keras.optimizers.Adam(clipvalue=cv)
    model.compile(loss='binary_crossentropy', optimizer=optimizer, metrics=['accuracy'])
    hist = model.fit(x_train, y_train, validation_split=0.2, epochs=epochs, batch_size=bs, verbose=0)
    return model, hist
for param_set in tqdm(param_grid):
    model, hist = training_run(param_set)
    test_results = np.mean(np.asarray([model.evaluate(x_test, y_test, verbose=0) for _ in range(10)]), axis=0)
    new_line = param_set
    new_line['test_loss'] = test_results[0]
    new_line['test_accuracy'] = test_results[1]
    results = pd.concat([results, pd.DataFrame(new_line, index=[0])], ignore_index=True, axis=0)
Explanation: Hyperparameter tuning
In this section we run the experiments over a grid of, for BNN important, parameters.
End of explanation
results_ = results.drop(columns=['test_loss']).fillna(-1)
sns.pairplot(results_, y_vars=['test_accuracy'], x_vars=['loss_scaling', 'batch_size', 'clipvalue'])
Explanation: Analysing the impact of hyperparameters
The following plot shows the test accuracy over the hyperparameters. It's clear that scaling the loss of the Bayesian layers in the model is necessary to reach convergence.
To plot those results NaN values in the clipvalues column are set to -1. In the corresponding runs gradients were not clipped.
What is the correct scaling factor for the regularizers?
BNN use the Kullback-Leibler Divergence (KL-divergence) as a regularization term to minimize the difference between the prior distribution $p$ and the approximated posterior distribution $q$ over the parameters of the network. This divergence term is written as $KL(q||p)$.  
The KL-divergence should be included once for an entire dataset. Due to the implementation  in tensorflow the regularization term is added to the value of the loss function after every batch though, leading to a larger KL-divergence term compared to the loss. To bring both terms on the same scale, we apply a scaling factor < 1 to the KL-divergence. The loss for each batch is divided by the batch size and therefore $KL(q||p)$ (which is in relation to the complete dataset) is divided by the number of training samples (in this case x_train.shape[0]).
End of explanation
results_ = results_[results_['loss_scaling'] != 1.]
results_lowclip = results_[results_['clipvalue'] == 0.1].drop(columns=['loss_scaling'])
results_highclip = results_[results_['clipvalue'] == 0.5].drop(columns=['loss_scaling'])
plt.scatter(results_lowclip['batch_size'], results_lowclip['test_accuracy'], c='b', label="clipvalue=0.1")
plt.scatter(results_highclip['batch_size'], results_highclip['test_accuracy'], c='r', label="clipvalue=0.5")
plt.xlabel("batch size")
plt.ylabel("accuracy")
plt.legend()
plt.show()
results_ = results_[results_['loss_scaling'] != 1.]
results_64 = results_[results_['batch_size'] == 64].drop(columns=['loss_scaling'])
results_128 = results_[results_['batch_size'] == 128].drop(columns=['loss_scaling'])
results_256 = results_[results_['batch_size'] == 256].drop(columns=['loss_scaling'])
plt.scatter(results_64['clipvalue'], results_64['test_accuracy'], c='b', label="batch_size=64")
plt.scatter(results_128['clipvalue'], results_128['test_accuracy'], c='r', label="batch_size=128")
plt.scatter(results_256['clipvalue'], results_256['test_accuracy'], c='g', label="batch_size=256")
plt.xlabel("clipvalue")
plt.ylabel("accuracy")
plt.legend()
plt.show()
Explanation: Choosing the correct batch size as well as possibly a value to clip gradients by during the training is more complicated.
From those experiments it looks like a smaller batch size will be benefitial, even though this is also dependent on the dataset. Running this experiment with less samples shows that a larger batch size can greatly increase accuracy. A possible explaination for this behaviour is that in cases where less data is available, larger batch sizes compensate for the additional variance that is caused by the Bayesian properties.
Using this dataset, clipping the gradients doesn't seem to have a huge affect on the outcome although it is often an important technique with Bayesian LSTMs.
End of explanation
best_params = results_.iloc[np.argmax(results_['test_accuracy'])].to_dict()
if best_params['clipvalue'] < 0:
    best_params['clipvalue'] = None
model, hist = training_run(best_params)
plot_hist(hist)
Explanation: Training the best Bayesian Model
End of explanation
<END_TASK>