import tensorflow as tf

    for h, n_out in enumerate(layers_width[1:]):
        # Hidden layer
        if(h < len(layers_width)-2):
            # Perform affine mapping at each layer of the neural network
            Z = tf.layers.dense(Z, n_basis//2)
            # Define variational parameters
            alpha_mean = tf.get_variable('alpha_mean_layer'+str(h),
                shape=[1, 1, n_basis, n_out],
                initializer=tf.random_normal_initializer())
            alpha_logstd = tf.get_variable('alpha_logstd_layer'+str(h),
                shape=[1, 1, n_basis, n_out],
                initializer=tf.random_normal_initializer())
            alpha_std = tf.exp(alpha_logstd)
            # Compute epsilon from {n_samples} standard Gaussian
            # epsilon = tf.random_normal([n_samples, 1, n_out*2, n_out])
            epsilon = tf.random_uniform([n_samples, 1, n_basis, n_out])
            hyp_params = tf.get_variable('hyp_params_layer'+str(h),
                shape=[2],
                initializer=tf.random_normal_initializer())
            l1, l2 = tf.nn.sigmoid(hyp_params[0]), tf.exp(hyp_params[1])
            epsilon = tf.sinh(epsilon*l2)/tf.cosh(epsilon*l2)**l1/l2
            # Compute A_{h+1}
