{'fraction_non_alphanumeric': 0.0312231529209622, 'fraction_numerical': 0.006819158075601375, 'mean_word_length': 4.656643887623386, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Propagating Uncertainty through the tanh Function with Application to Reservoir Computing

Many neural networks use the tanh activation function, however when given a probability distribution as input, the problem of computing the output distribution in neural networks with tanh activation has not yet been addressed. One important example is the initialization of the echo state network in reservoir computing, where random initialization of the reservoir requires time to wash out the initial conditions, thereby wasting precious data and computational resources. Motivated by this problem, we propose a novel solution utilizing a moment based approach to propagate uncertainty through an Echo State Network to reduce the washout time. In this work, we contribute two new methods to propagate uncertainty through the tanh activation function and propose the Probabilistic Echo State Network (PESN), a method that is shown to have better average performance than deterministic Echo State Networks given the random initialization of reservoir states. Additionally we test single and multi-step uncertainty propagation of our method on two regression tasks and show that we are able to recover similar means and variances as computed by Monte-Carlo simulations.


Introduction
Neural networks can often have inputs that are uncertain. Input uncertainty can come from measurement error, adversarial noise [1], or even from output feedback. The majority of work on the subject of uncertainty in neural networks revolves around uncertainty in the model, not necessarily uncertainty in the input itself. Bayesian neural networks perform inference using a prior over the weights of the network, while dropout [2] samples iterations of the network with a probabilistic mask to build the posterior distribution of the model. In both cases the output uncertainty emerges as an explicit function of the model uncertainty, while the uncertainty of the input is classified as aleatoric and not explicity propagated through the model. For non-parametric probabilistic methods, such as Gaussian process regression [3], special care must be taken for multi-step prediction with uncertain inputs [4]. In this paper, we focus on addressing the problem of uncertainty propagation through the tanh function, which is a popular choice for activation function in neural networks, in particular the Echo State Network.
In light of the challenges in input uncertainty propagation and their role in recurrent neural networks, we aim to contribute the following: • A theoretical and numerical analysis for 3 methods of propagating input uncertainty through the tanh activation function, with an extension to other nonlinear activation functions.
• A new method, named the Probabilistic Echo State Network (PESN), which aims to reduce the time required to achieve the echo state property.
Related Works: There exists recent work on propagating input uncertainty through feed-forward bayesian neural networks 5. Here, the authors perform approximate inference to propagate uncertain inputs through feedforward neural networks for classification tasks. We differentiate our work in three fundamental ways: 1) we consider the tanh function which is not addressed in [5], 2) a novel contribution utilizing splines to propagate gaussian input uncertainty through continuous activations and 3) a general focus on improving the reservoir computing framework. To the best of the authors' knowledge we are the first to utilize a spline approximation to the integrand of an expectation of a gaussian in order to perform approximate inference. Similar ideas include [6], where the authors utilize b-splines to approximate the density in order to perform approximate inference. Work in improving reservoir computing has traditionally focused on the structure of the reservoir, such as optimizing hyperparameters [7] or finding the minimum reservoir size [8]. This work is one of the first to attempt to reduce the time required to converge to the echo state property.

Reservoir Computing
Reservoir computing (RC) is a paradigm for training recurrent neural networks (RNNs) [9]. It was introduced in early 2000's by Jaeger under the name 'Echo State Networks' (ESN) for time-series predictions [10], and by Maass under the name 'Liquid State Machines' [11] for modeling computation in biological networks of neurons. Despite being developed from very different communities, these two approaches are largely mathematically identical. In this work we focus on ESN since we aim to improve the efficiency of RNN learning for engineering applications.
ESN differs from other RNNs in terms of its training scheme. Generally, RC consists of two steps: 1) drive a network with sparse and fixed connections with an input and output sequence. 2) Train the output (or readout) layer so that the network output is similar to the teacher output. The readout layer is usually trained using regression techniques such as linear regression [9] and Gaussian process regression [12]. One problem is determining the initial state of the hidden layer (reservoir). According to the echo state property [13], the effect of initial conditions can be 'washed out' therefore the state can be initialized randomly. However, there are three main drawbacks: first, a significant amount of training data is wasted because the initial period of training run needs to be discarded. Second, at test time, an initial input sequence needs to be fed into the network before using it for prediction tasks. Third, the state forgetting property is usually not guaranteed so the performance of the ESN may still depend on the initial state.

Propagating Uncertainty through the tanh
In our work we analyze three distinct methods to propagate a gaussian input through the tanh activation function. The simplest and most well known method is simply Monte-Carlo (MC) sampling where we sample the input distribution then pass each sample through the activation to compute an estimate of the moments. While easy to implement and understand, Monte Carlo has an obvious drawback in terms of computational time. The second method utilizes a well known approximation to the tanh activation function in the form of the logistic cumulative distribution function. The connection between the logistic and gaussian distributions are utilized to approximate the mean of the activation output. The variance approximation fits the moments of the Gaussian pdf to the function (1 − tanh(x) 2 ) in order make the expectation tractable. The cross covariance terms are ignored in our derivation. While this method is certainly faster, the accuracy is limited by our approximations. The final method leverages spline approximations and analytical expressions for the expectation of polynomial functions. This method strikes a balance between computational complexity and accuracy by adjusting the width of the spline mesh. We compare absolute error of the moment approximations, computational complexity of the two analytical and spline methods, and provide error bounds for the spline approximation.

Analytical Approximation to Mean and Variance
First we relate the tanh function with the logistic cumulative distribution function CDF, and approximate the logistic distribution with an appropriate gaussian distribution. We assume that the dimensions of the logistically distributed random vector x are independent given the location parameter z. The expectation and variance of x given z are E(x|z) = z, Var(z|z) = π 2 12 I respectively. We can use moment matching to approximate x with a gaussian Here the logistic distribution is parameterized by the location λ and the scale ξ. The tanh function is a function of a logistic CDF with location λ = z and scale ξ = 0.5. Additionally let us assume that z ∼ N (µ, σ). The mean of the distribution, p(tanh(z)), is given by: Next, we ignore the cross covariance terms between the elementwise tanh, and directly compute the diagonals of the output covariance matrix. The authors in [14] investigated different approximations to the error function. One preliminary solution was the use of the tanh function to approximate the error function. We can take advantage of this approximation and the appropriate scaling factors to approximate the function (1 − tanh(x)) 2 . We use a Gaussian pdf where σ 2 = 2 π and y = 2 √ π x. Then we can again apply similar algebraic manipulations to compute the variance terms. The full derivation for the mean and variance expressions can be found in the supplementary material.

Spline Approximation to Mean and Variance
In this section, we present a formulation for computing the moments of the tanh transformation utilizing a spline approximation of the argument of the expectation. Again let us assume that z ∼ N (µ, σ), the expectation and variance of the elementwise tanh are given by the following equations: Let us assume that we are dealing with scalar inputs to the activation function, since the vector case is handled elementwise and we are ignoring the off-diagonal elements of the input variance. Instead of direct numeric integration, we propose taking advantage of three facts.
• The tails of the Gaussian pdf approach zero, thus when multiplied against a function with constant tails, such as the tanh, the integrand also goes to zero. Thus we can approximate the integral over the real line using an integral over some compact subset centered around the mean of the input. • The tanh function is continuous so we can uniformly approximate it with polynomials on the compact subset of integration. • We can derive analytic forms for the definite integral of a product of polynomials and exponentials.
As a result, we can derive expressions for the mean and variance in terms of these analytic integrals. Let f (z) be a scalar valued function that is continuous on the real interval [a, b]. Additionally, let P (z) be a piecewise continuous cubic polynomial that interpolates f with N nodes, where each node is the beginning of intervals [z j , z j+1 ], where z 0 = a and z N = b.
Within the inner sum, the four terms can be computed analytically utilizing expressions for the integral of the product between polynomials and exponentials. We provide the full expressions for polynomials up to 3rd order in the supplementary material. The nodes and coefficients can be computed prior to evaluation of the network, and the integrals are computed once at every layer. This method is in fact general for any continuous function f (z) whose tails are bounded, or polynomial outside of a compact set. Additionally, we can move forward to approximate higher moments of the output distribution and get expressions for the skew and kurtosis, at the cost of having spline approximations for f (z) 3 and f (z) 4 . We provide expressions for higher moments and for alternative activation functions in the supplement. The mean and variance are given below, and Algorithm 1 describes how to use this technique to compute the moments of the activation.

Computational Complexity Analysis
Next we look at the computational complexity of both the analytical and spline method in terms of the number of hidden states in the network. Let D be the number of hidden states in the network and let N be the number of mesh points for the spline approximation. Since the spline method makes frequent use of the error function, to approximate the operation count, we use the error function approximation introduced in [15]. We can see from Table 1

Error Bounds for Mean and Variance Estimation
In this section we provide an upper bound on the mean and variance estimates as a function of the input mean, input variance, the number of nodes in the mesh, and the location of the mesh in space. Let g(z) ∈ C 4 , and define τ as the maximum interval length in the mesh spanning some compact set z ∈ [a, b]. Then for a cubic spline approximation P (z), we can place an upper bound on the maximum error over the interval [a, b] [17].
Algorithm 1 (µ, σ) approximation of tanh(z) with spline method Input: a, b, n points 1: z ← evenly spaced n points numbers over the interval[a, b] 2: c 1 , c 2 ← cubic spline coefficients of Eqs. (53), (62) Let us use g(z) = tanh(z). Given the fact that z ∼ N (µ, σ), we can compute an error bound over the mean of the output distribution by applying the triangle inequality. where Let us additionally compute a bound for the variance. Let µ true denote the true value of the mean, µ spline be the approximated value of the mean. Additionally, letg(z) be our polynomial approximation of tanh(z) 2 . Then we go to the expression for variance Equation (78): 1 is the spline error bound where the unknown function is tanh(z) 2 . Additionally, we know that |µ spline − µ true | < µ (a, b, µ, σ, τ ) and that since we are dealing with the tanh activation function, the activation mean is upper bounded by 1. Utilizing this method allows us to design the number of uniformly spaced mesh points and the size of the interval according to a desired accuracy. For reference, we utilize the interval [-10, 10] for the mean of the tanh activation, and use 101 mesh points, resulting in an error bound of 4.21321e−5 for the mean of an input distribution with mean 3 and variance 0.2. Again, the full derivation for both error bounds can be found in the supplement.

Probabilistic Echo State Networks
The proposed algorithm, Probabilistic Echo State Networks (PESN) derives its internal equations from the deterministic echo state network [13]. Let z k ∈ R Nx+Nu be input to the network, y k ∈ R Nx be the output of the network, h k ∈ R N h is the hidden state at time k. The deterministic echo state network is characterized by 5 hyperparameters: reservoir size N h , leak rate L, noise magnitude M e , sparsity fraction s, and spectral radius r.
Here dw is Gaussian noise with variance 1. The sparsity fraction refers to the fraction on nonzero elements in the matrix W , while the spectral radius parameter is the spectral radius of W . For the echo state network, we randomly generate the internal weight matrices W in , W fb , and W , and we only train the readout weight matrix W out . In our case, training is accomplished using batch least squares. At test time, we recursively update the readout matrix using recursive least squares regression. [18] For probabilistic echo state networks, assume that the input to the networks are Gaussian. We use moment matching to propagate the Gaussian uncertainty through the nonlinearity in the network, thus we retrieve a Gaussian output. Our notation is as follows: a Gaussian random vector z has a mean µ z and covariance matrix Σ z . We assume that all Gaussian vectors in our echo state network are additionally jointly Gaussian, so that linear combinations of these Gaussian vectors produce another Gaussian vector with appropriate mean vector and covariance matrix. Define the input as z(k), the previous target as y(k − 1) and the hidden state as h(k − 1).
Using our element-wise approximation, we drop the off-diagonal terms in a and compute the approximate mean µ tanh(a) and variance Σ tanh(a) using Equations (77) and (78). We also ignore the cross correlation terms, which will cause underestimation of uncertainty.
Here N h , L, M e , s, r are the five hyperparameters for the echo state network, z and y are input and output data, generated or collected from experiments, n train and n test are the lengths of the train and test data, λ is a forgetting factor, P is δI where I is the identity matrix. Algorithm 2 is used for training. Algorithm 3 is used for single and multi-step prediction for the probabilistic echo state network.

Numerical Analysis of Moment Approximation
First we compare the absolute error between the mean and variance estimated from Monte Carlo against the mean and variance from spline method and analytic method. We evaluate our moment propagation on a range input means and select variance values to assess accuracy. At each point, we have an input Gaussian distribution determined by a mean and variance. We sample this input distribution and compute an estimated mean and variance through Monte Carlo, then compare our two methods. Figure 4 shows that the spline method results in a more accurate approximation, and as the input saturates (mean approaches 5 and -5), the error quickly drops. In general we see the analytic approximation over estimates the variance of the tanh output. Full contour plots of the spline and analytical comparison for a grid of both means and variance are available in the supplementary materials. We see in Figure 3 how computational time effects absolute error of moments. The analytic method is constant, but has comparatively high error compared to the others. The spline method has higher initial computational cost, but converges in error quickly.

Effect of Probabilistic Hidden State on Washout Period
In Table 2 we numerically analyze the difference in multistep prediction error between the probabilistic echo state and the deterministic echo state algorithms as a function of the number of washout points. Each of the 50 trials were run with varying washout lengths (in timesteps), where we ran 50 samples of the deterministic echo state algorithm and compared the average absolute error between the ground truth and predicted state. The table shows the the probabilistic echo state network outperforms the monte-carlo deterministic echo state in absolute error for the majority of washout lengths. When the Monte Carlo deterministic ESN outperforms the PESN, it does so by a fairly small margin. This implies that the PESN reduces the time required to attain the "echo state property" that is characteristic for echo state networks.  In the "mean" column, the minimum value between probabilistic and monte-carlo is bolded. In the "Washout Length" column, the minimum error length is bolded.  In an attempt to quantify the level of synchronization in the hidden states of both the ESN and the PESN, we can create a probability distribution of the values of the reservoir values during washout. Synchronization would imply that there is less uncertainty in the values of the hidden state, thus we can quantify this via the Shannon Entropy of this probability distribution. Table 3 lists the Shannon entropy values for hidden state trajectories. Intuitively, the entropy is small for every washout length, which implies there is less randomness to "wash out".

Model Learning
We test the PESN on the task of learning dynamics of cart pole and Gazebo ARDrone. We compare the PESN against Monte Carlo estimates obtained with the deterministic ESN. The PESN successfully propagates uncertainty for single step predictions and for up to 20 timesteps for multi step prediction. The results, found in the supplementary material Section 3.3, demonstrate the ability of our method to capture input uncertainty as it passes through the recurrent network.

Conclusions
In this work, we investigate two different methods of propagating uncertainty through the tanh function. We demonstrate that the analytical method is the fastest way to propagate uncertainty through the tanh function, however it tends to over estimate uncertainty. The spline method is more general and has tuneable accuracy, however the computational cost is higher. Utilizing these new developments, we propose the probabilistic echo state network (PESN), which attempts to solve a fundamental problem of reservoir computing, which is the time required to fulfill the asymptotic behavior of the echo state property. We are able demonstrate, for multi step regression tasks, the PESN has better performance than a deterministic ESN (as measured by absolute error). The better performance is accompanied by a lower shannon entropy in the distribution of hidden state value distribution, which we believe indicates faster convergence to the echo state property. We additionally test the PESN method on ARDrone data collected from the Gazebo simulator, to demonstrate that the method can scale to higher dimensional systems and still propagate uncertainty effectively. However, these methods are not without their drawbacks. For the spline method, input, output scaling and spectral radius of the PESN must be tuned fit inside interval [a, b] so that we satisfy our moment error bounds. The multi-step regression is heavily dependent on incremental updating of the output weights in order to maintain low absolute error. Both the spline and analytic methods are tractable techniques to propagate uncertainty through the tanh function, and the resulting Probabilistic Echo State network is able to improve the convergence to the echo state property.

Supplementary Material
This is the supplementary material for the paper on Propagating Uncertainty through the Tanh function for Reservoir Computing. The supplement is organized as follows: Section 6.1 includes derivations for both the analytic and spline methods, derivations for the error bound for mean and variance of the spline tanh approximation. Section 6.2 repeats the equations for the Probabilistic Echo State Network (PESN). Section 6.3 contains additional numerical examples. These include propagation of uncertainty through a feedforward neural network for regression, dynamical systems regression, and activation function comparisons.
6.1 Propagating Uncertainty through Tanh

Analytical Approximation of Mean and Variance
First we present some useful identities and approximations: The cumulative distribution function (CDF) of a logistic distribution over a random variable x is: Here the logistic distribution is parameterized by the location λ and the scale σ. The tanh function can be written as follows and relation can be derived: Thus the tanh function is a function of a logistic CDF with location λ = z (the input to the tanh) and scale σ = 0.5. Let us assume that the dimensions of the logistically distributed random vector x are independent given the location parameter z. From the definition of the logistic distribution, the mean and variance of the conditional distribution p(x|z, 0.5) are: We also make use of the formula for the product of two Gaussian distributions( [19]). where Let us assume that z ∼ N (µ, σ). We want to compute the moments of the distribution: p(tanh(z)). We will make use of moment matching to approximate a logistic distribution with a Gaussian, and vice versa.

E(tanh(z)) = tanh(z)p(z)dz (43)
In Equation (50), we use the product rule of Gaussians, and specify: µ a = x, µ b = µ, Σ a = I π 2 12 , and Σ b = σ. Then we are left with the constant coefficient in front of a Gaussian pdf of the input z. As a result, the constant coefficient is not a function of the input, thus the integral over this Gaussian pdf integrates to 1.
Thus, given the input mean and diagonal covariance function, we can compute the posterior mean of an elementwise tanh operation.
Variance: The author in [14] investigated different approximations to the error function. One preliminary solution was the use of the tanh function to approximate the error function. Using this fact, we see that with the appropriate scaling factors, the Gaussian pdf can be an approximation to the function (1 − tanh(x)) 2 . Again, we are ignoring cross correlation terms, thus the output variance will be diagonal and we can compute each term separately.

Spline Approximation of Mean and Variance
In this section, we present a formulation for computing the moments of the tanh transformation utilizing a spline approximation of the argument of the expectation. Again let us assume that z ∼ N (µ, σ), the expectation and variance of the elementwise tanh are given by the following equations: In both cases the integrands are not easily integrable, and depend on the mean and variance of the input. Let us assume that we are dealing with scalar inputs to the activation function, since the vector case is handled elementwise and we are ignoring the off-diagonal elements of the input variance. Instead of direct numeric integration, we propose taking advantage of three facts. First, we know that the tails of the Gaussian pdf approach zero, thus when multiplied against a function with constant tails, the integrand also goes to zero. The tanh, sigmoid both have constant tails. For the ReLu and swish activation functions, the negative tail is 0, and the positive tail is x, which has an analytical form when integrated against a Gaussian. Thus we can approximate the integral over the real line using an integral over some compact subset centered around the mean of the input. Second, most activation functions are continuous, and on the compact subset of integration, we can uniformly approximate these functions with polynomials, in fact we will utilize piecewise cubic polynomials for this approximation. Finally, we can derive analytic forms to the integrals of polynomials multiplied against exponentials. The results are expressions for both the output mean and variance that are functions of the exponential function, error function, and input mean and variance. Let f (z) be a scalar valued function that is continuous on the real interval [a, b]. Additionally, let P (z) be a piecewise continuous cubic polynomial that interpolates f with N nodes, where each node is the beginning of intervals [z j , z j+1 ], where z 0 = a and z N = b.
Within the inner sum, the four terms can be computed analytically, however let us first make a change of variables to simplify the expressions. Let x = (z−µ) √ 2σ , then dz = √ 2σdx. Additionally we should transform our bounds: α = 2σ . Next we can compute the integrals: Thus using Equation (68), and Equations (70) -(76), we can analytically compute the expectation of any smooth function f (z) whose tails are constant, or polynomial outside of a compact set. The nodes and coefficients can be computed prior to evaluation of the network, and the integrals are computed once at every layer. The mean and variance the functions to approximate are f (z) and f (z) 2 respectively. Both functions have a distinct set of coefficients, but we use the same set of equally spaced nodes for each of computation. Additionally, we can move forward to approximate higher moments of the output distribution and get expressions for the skew and kurtosis, at the cost of having spline approximations for f (z) 3 and f (z) 4 . We can express the mean µ, variance σ, skewness γ and kurtosis κ.
Error Bounds for Spline Mean and Variance Estimation: In this section we provide an upper bound on the mean and variance estimates as function of the input mean, input variance, the number of nodes in the mesh, and the location of the mesh in space. Let g(z) ∈ C 4 ( some examples of common activation functions with fourth derivatives are tanh, linear, sigmoid, and swish), and define τ as the maximum interval length in the mesh spanning some compact set z ∈ [a, b]. Then for a cubic spline approximation P (x), we can place an upper bound on the maximum error over the interval [a, b] [17].
Let us use g(z) = tanh(z) as an example to compute an expression for the error for an approximation of E tanh(z).g Now we approximate expectation of the function g(z), given the fact that z ∼ N (µ, σ), by splitting the integral into pieces where the function is approximately constant outside of the interval [a, b], and the integral over [a, b] itself. Then we can bound the approximations of g(x) by their maximum values, pass the constants outside the integral, and evaluate integrals over a Gaussian distribution. Let us define the following: Now we can compute an error bound over the mean of the output distribution by applying the triangle inequality, assumptions about the tail behavior of tanh, our approximation over the interval [a, b].
. We can generalize this bound to any function in the expectation, we simply have to compute the appropriate polynomial coefficients and have an assumption about the tail behavior of the function. For example, let us additionally compute a bound for the variance. Let µ true denote the true value of the mean, µ spline be the approximated value of the mean. Additionally, letg(z) be our polynomial approximation of tanh(z) 2 . Then we go to the expression for variance Equation (78): In this case, we apply Equation (91) with the coefficients of the integrand function g(z) = tanh(z) 2 , thus we can come up with a bound for the following integral in Equation (95): Additionally, we know that |µ spline − µ true | < µ (a, b, µ, σ, τ ) and that since we are dealing with the tanh activation function, the activation mean is upper bounded by 1. Thus: Utilizing this method allows us to design the number of uniformly spaced mesh points and the size of the interval according to a desired accuracy. For reference, we utilize the interval [-10, 10] for the mean of the tanh activation, and use 101 mesh points, resulting in an error bound of 4.21321e−5 for the mean of an input distribution with mean 3 and variance 0.2.

Probabilistic Echo State Networks
The proposed algorithm, Probabilistic Echo State Networks (PESN) derives its internal equations from the deterministic echo state network [13]. Let z k ∈ R Nx+Nu be input to the network, y k ∈ R Nx be the output of the network, h k ∈ R N h is the hidden state at time k. The deterministic echo state network is characterized by 5 hyperparameters: reservoir size N h , leak rate L, noise magnitude M e , sparsity fraction s, and spectral radius r.
(97) Here dw is Gaussian noise with variance 1. The sparsity fraction refers to the fraction on nonzero elements in the matrix W , while the spectral radius parameter is the spectral radius of W . For the echo state network, we randomly generate the internal weight matrices W in , W fb , and W , and we only train the readout weight matrix W out . In our case, training is accomplished using batch least squares. At test time, we recursively update the readout matrix using recursive least squares regression. [18] We now present expressions for probabilistic echo state networks. Here we assume that in that the input to the networks are Gaussian, we use moment matching to propagate the Gaussian uncertainty through the nonlinearity in the network, and we are able to retrieve a Gaussian output. Our notation is as follows: a Gaussian random vector z has a mean µ z and covariance matrix Σ z . We assume that all Gaussian vectors in our echo state network are additionally jointly Gaussian, so that linear combinations of these Gaussian vectors produce another Gaussian vector with appropriate mean vector and covariance vector. That is if a and b are Gaussian random vectors that are also jointly Gaussian, then: (100) Analogously to Equation (96) if the input is given by z(k), the previous target by y(k − 1) and the hidden state by h(k − 1), the argument to the activation function is given by Using our element-wise approximation, we drop the off-diagonal terms in a and compute the approximate mean µ tanh(a) and variance Σ tanh(a) using Equations (77) and (78).
Here N h , L, M e , s, r are the five hyperparameters for the echo state network, z and y are input and output data, generated or collected from experiments, n train and n test are the lengths of the train and test data, λ is a forgetting factor, P is δI where I is the identity matrix, W in , W out , W fb , and W are the weights of the network, h is the hidden state.

Activation Function Moment Comparison
First we compare the absolute error between the mean and variance estimated from Monte Carlo against the mean and variance from spline method and analytic method. We evaluate our moment propagation on a grid of input means and variances to assess accuracy on a range. At each gridpoint, we have an input Gaussian distribution determined by a mean and variance. We sample this input distribution and compute an estimated mean and variance through Monte Carlo, then compare against our two methods. Figure 4 shows that the spline method results in a more accurate approximation, and as the input saturates (mean approaches 5 and -5), the error quickly drops logarithmically. In general we see the analytic approximation over estimates the variance of the tanh output.
Additionally, we can compute the skew and kurtosis any continuous activation function. The expressions for these moments are given in Equations (77) -(80). We present the absolute error in moments for mean, variance, skew, and kurtosis for 3 activation functions utilizing the spline method: tanh, sigmoid, and swish (where β = 1 for the swish activation function) in Figure 5.

Multiple Layer Uncertainty Propagation
In this section, we have constructed two feedforward neural networks with the tanh activation function to demonstrate the performance of the moment approximations. In this case, we randomly initialized all weight matrices. In both cases, the network input is 1024 units and the network output is a single unit. Network 1 has a width of 5 units for 5 hidden layers and network 2 has a width of 50 units and 10 hidden layers. The uncertain input is Gaussian distributed with mean -0.5, and variance 0.01. We see in Table 4a that for network 1, the output is close to a delta function, as a result of the saturation of neurons. The analytical method over estimates the variance in this case, while the spline method is able to recover the moments well. For network 2, shown in Table 4b, the distribution goes from a delta function to Gaussian, and both the analytical method and the spline method can recover the output distribution. We see in the table, the absolute error of approximation is smaller for the spline case. Our error criteria are defined for a given approximation X, layer i, and number of hidden units N hi . µ M C,i,j is the Monte Carlo approximation mean at layer i, hidden unit j. The Monte Carlo We have large errors in skew and kurtosis wherever the variance is less than 10e − 6, since these moments are standardized using the variance.
approximation was computed with 50000 samples.   In the shallow network, the output cdf is close to a delta function. In the deep and wide network, the output cdf is close to a Gaussian distribution.

Model learning
In the field of robotics and autonomy, recently there is an increasing interest in using machine learning to get the accurate model approximation. Model accuracy is important for many applications including planning and control. However, analytical models cannot catch the complexity of a system very well, so we approximate the system model by learning from data. In this section, we learn unknown system dynamics from data using probabilistic echo state networks(ESN) and evaluate the learned model. We suppose our systems' dynamics are unknown and represent the dynamics as the following differential equation   ,where x ∈ R Nx is the state and u ∈ R Nu is the control. Given the state-control pair and the state transition dx, we are able to infer the dynamics as well as the state transition. With the goal of inferring dx given (x, u), we compare our probabilistic ESN using spline approximation of tanh to Monte Carlo(MC) sampling of the deterministic ESN using tanh under uncertain inputs. To make the inputs uncertain for the deterministic ESN, we added Gaussian noise to the state x when we do the MC sampling. For fair comparison, the two networks shared same weights W in , W out , W f b , and W . Data for the Cart Pole system was generated by equations of motion derived using first principles and the ARDrone data was collected using a Gazebo simulator [20]. As the kinematics of each system is straightforward and easy to derive, we directly use the kinematics equations to predict next position states and we only learn the dynamics of the systems with ESNs.
Cart Pole: We consider the Cart Pole having a state vector x = [x, θ,ẋ, ω] T , where x is cart position, θ is pole angle,ẋ is cart velocity and ω is pole angular velocity. The control input to the system u is the force applied to the cart. Through several times of simulation with different controls, we collected various data. We use [x T , u] T as an input to the network. To learn the dynamics of the system, we set the target output as the two velocity terms [ẋ, ω] T .
ARDrone: In the simulated ARDrone, we use 12 state of the system. x = [x, y, z,ẋ,ẏ,ż, φ, θ, ψ,φ,θ,ψ] T , where (x, y, z) are positions, (φ, θ, ψ) are angles, (ẋ,ẏ,ż) are linear velocities and (φ,θ,ψ) are roll rate, pitch rate and yaw rate each. The control input to the ARDrone system is u = [v x,cmd , v y,cmd , v z,cmd , v ψ,cmd ] T . Here (v x,cmd , v y,cmd , v z,cmd ) are linear x, y, z velocity commands and v ψ,cmd is the angular z velocity command. Running the simulated ARDrone with various control commands, we collected 80,000+ data. Using [x T , u T ] T as an input to the network, we set the output target as velocity terms [ẋ,ẋ,ẏ,φ,θ,ψ] T to learn the dynamics of the drone.      variance of the MC-sampled deterministic ESN with uncertain inputs. However, in the multi step prediction, the predicted mean of probabilistic ESN diverges from the MC-sampled ESN after 4-5 time steps, but the variance of the distribution grows similarly.