© Revista de Matemática: Teoría y Aplicaciones 2022 29(2) : 289–312 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v29i2. 48885

deep gaussian processes and infinite neural networks for the analysis of eeg signals in alzheimer’s diseases

procesos gausianos profundos y redes neuronales infinitas para el análisis de

señales eeg en la enfermedad de alzheimer

Krishna Román^{∗}  Andy Cumbicus^{†}  Saba Infante^{‡}

Rigoberto Fonseca-Delgado^{§}

Received: 1/Nov/2021; Revised: 23/Jun/2022; Accepted: 28/Jun/2022

∗Yachay Tech University, School of Mathematical and Computational Sciences, Urcuquí, Ecuador. E-mail: krishna.roman@yachaytech.edu.ec

†Misma dirección que/Same address as:K.Román.E-mail:andy.cumbicus@yachaytech.edu.ec ‡Misma dirección que/Same address as: K. Román. E-mail: sinfante@yachaytech.edu.ec §Misma dirección que/Same address as: K. Román. E-mail: rfonseca@yachaytech.edu.ec

290 k. román — a. cumbicus — s. infante — r. fonseca-delgado

Abstract

Deep neural network models (DGPs) can be represented hierarchically by a sequential composition of layers. When the prior distribution over the weights and biases are independently identically distributed, there is an equivalence with Gaussian processes (GP) in the limit of an infinite network width. DGPs are non-parametric statistical models used to characterize patterns of complex non-linear systems due to their flexibility, greater generalization capacity, and a natural way of making inferences about the parameters and states of the system. This article proposes a hierarchical Bayesian structure to model the weights and biases of a deep neural network. We deduce a general formula to calculate the integrals of Gaussian processes with non-linear transfer densities and obtain a kernel to estimate the covariance functions. In the methodology, we conduct an empirical study analyzing an electroencephalogram (EEG) database for diagnosing Alzheimer’s disease. Additionally, the DGPs models are estimated and compared with the NN models for 5, 10, 50, 100, 500, and 1000 neurons in the hidden layer, considering two transfer functions: Rectified Linear Unit (ReLU) and hyperbolic Tangent (Tanh). The results show good performance in the classification of the signals. Finally, we use the mean square error as a goodness of fit measure to validate the proposed models, obtaining low estimation errors.

Keywords: deep Gaussian process; Alzheimer disease; electroencephalogram.

Resumen

Los modelos de redes neuronales profundos (DGPs) se pueden representar jerárquicamente mediante una composición secuencial de capas. Cuando la distribución prior sobre los pesos y sesgos son independientes idénticamente distribuidos, existe una equivalencia con los procesos Gaussiano (GP), en el límite de una anchura de red infinita. Los DGPs son modelos estadísticos no paramétricos y se utilizan para caracterizar los patrones de sistema no lineales complejos, por su flexibilidad, mayor capacidad de generalización, y porque proporcionan una forma natural para hacer inferencia sobre los parámetros y estados del sistema. En este artículo se propone una estructura Bayesiana jerárquica para modelar los pesos y sesgos de la red neuronal profunda, se deduce una formula general para calcular las integrales de procesos Gaussianos con funciones de transferencias no lineles, y se obtiene un núcleo para estimar las funciones de covarianzas. Para ilustrar la metodología se realiza un estudio empírico analizando una base de datos de electroencefalogramas (EEG) para el diagnóstico de la enfermedad de Alzheimer. Adicionalmente, se estiman los modelos DGPs, y se comparan con los modelos de NN para 5, 10, 50, 100, 500 y 1000 neuronas en la capa oculta, considerando dos funciones de transferencia: Unidad Lineal Rectificada (ReLU) y tangenge hiperbólica (Tanh). Los resultados demuestran buen desempeño en la clasificación de las señales. Finalmente, utilizó como medida de bondad de ajuste el error cuadrático medio para validar los modelos propuestos, obteniéndose errores de estimación bajos.

Palabras clave: procesos gausianos profundos; enfermedad de Alzheimer; electroencefalogramas.

Mathematics Subject Classification: 60G15, 60H35.

1 Introduction

Deep machine learning systems are computational algorithms that provide powerful modern tools that allow the use of mathematical models with complex structures by including multiple intermediate layers at different levels combined with transfer functions with non-linear systems.

Neural networks are machine learning models that have received a lot of attention in recent years due to their success in many real-world applications: they have been used very frequently in filtering content on social networks, in natural language processing, pattern recognition in Big data, tracking objects in a sequence of images or videos, face and voice recognition, human mobility, mass information dissemination in networks, electronic commerce and classification of relevant information, among many other applications, LeCun et al. (2015) [13], Goodfellow et al. (2016) [8], Mosavi et al. (2020) [20], and Luca et al. (2020) [15].

In this work, a combined technique of Gaussian processes (GP) with neural network models is used, following the approach outlined by: MacKay (1992) [16], MacKay (1995) [?], Neal (1996) [21], and Williams (1996) [27]. GPs are used to model functional data because they are flexible, robust to outliers, and estimate calibrated uncertainty. Deep Gaussian Processes (DGP) are a generalization of a multilayer neural network viewed as a GP in the limit width. Salimbeni and Deisenroth (2017) [?] used an inference algorithm variational doubly stochastic algorithm that does not force independence between layers. Also, they demonstrated that a DGP model could be used effectively for many data points. They provide strong empirical evidence of the inference scheme for DGPs and show that they work well in practice in both classification and regression. Zhao et al. (2021) [30] used a state-space model for the deep Gaussian process (DGP) regression. They constructed a DGP by placing a GP prior transformed into the length and magnitude scales at each hierarchy level. They used a posterior maximum estimation procedure based on filtering algorithms and demonstrated the performance using non-stationary synthetic and gravitational wave signals.

Wilson and Izmailov (2020) [28] showed that deep ensembles effectively estimate the marginal distributions. They also investigated priors on the considered functions by defining vague priors on the neural network weights. In addition, they demonstrate properties, generalize the models from a probabilistic perspective, and obtain results equivalent to those reproduced using a GP. Lee et al. (2018) [14] demonstrated an equivalence between infinitely wide deep networks and GPs, developed computationally efficient methods to compute the covariance function of GPs and connected GPs with the theory of signal propagation in random neural networks. Investigations in the same direction highlight the work of Schoenholz et al. (2017) [26], Matthews et al. (2018) [18], Novak et al. (2019) [22], Garriga-Alonso, et al. (2019) [?], Agrawal et al. (2020) [2], Damianou and Lawrence (2013) [6], Khan et al. (2019) [11], Hazan and Jaakkola, (2015) [9]. We can mention related works to these topics: Infante et al. (2008) [10], Cedeño et al. (2021.a) [3] and Cedeño et al. (2021.b) [4].

The contribution in this work is based on the proposal of a hierarchical Bayesian structure to model the weights and biases of a neural network with infinite width. We established a general formula to calculate the integrals of Gaussian processes for sigmoid functions with non-linear structures, and we obtained a kernel to update the covariance functions. Finally, we made an application on a set of electroencephalogram signals from patients with Alzheimer’s disease.

The rest of the article is as follows: the Deep Neural Networks are defined in Section 2, the Gaussian processes are described in Section 3, the connection between the Deep Neural Networks and Gaussian Processes is made in Section 4, and the results are shown in the Section 5, and finally some discussions and conclusions are established in the Section 6.

2 Deep neural networks

Mathematically a neural network can be defined as a directed graph with vertices representing neurons and edges expressing connections. Each neuron’s input is a weighted sum of the output of all previous layer’s neurons connected to the input. There are many variants of neural networks that differ by their architecture. The simplest of these forms is the forward neural network, also known as the feedforward neural network.

Deep neural networks compose computations performed by many layers. Denoting the output of hidden layers by h_{i}^{(l)}, where l ∈ {1, …, L}, i ∈ {1, …, N_{l}} denote the indices of the neuron within the layer that receives the information from the neurons of the previous layer h_{j}^{(l−1)}, j ∈ {1, …, N_{l−1}}.

The output j represented by the j−th neuron in the output layer, is connected to the input vector x via a biased weighted sum and an non-linear activation function ϕ. The j−th component of the network output, h_{j}^{(l)}, is computed as:

h_{j}^{(l)} (x) = ϕ (z_{i}^{(l−1)}(x))  (1)

where:

z_{i}^{(l)} (x) = b_{i}^{(l)} + ∑_{j=1}^{N_{l−1}} W_{ij}^{(l)}h_{j}^{(l−1)}(x).  (2)

For convenience, the parameters of the neural network are combined into a vector of parameters θ = (b^{(1)}, W^{(1)}, …, b^{(L)}, W^{(L)}) and input data h^{(0)} = x = (1, x_{1}, …, x_{N_{l}}). Deep neural networks compose computations performed by many layers. The calculation for a network with L hidden layers is:

ŷ = f [h^{(L)} (z^{(L)} (h^{(L−1)} (… z^{(2)} (h^{(2)} (z^{(1)} (h^{(1)} (z^{(0)})))))))] . (3)

When L is large, it is called a deep neural network, and each pre-activation function z^{(L)}(x) is typically a linear operation with the matrix W^{(L)} and bias b^{(L)}, which can be combined with the parameters θ. An architecture of a feedforward neural network with L hidden layers is shown in Figure 1.

Figure 1: An architecture of a feedforward neural network with l hidden layers, (Witten et al., 2016) [29].

3 Gaussian processes

Assume we have access to a training dataset of n input-output observations

𝔻 = {(x_{i}, y_{i}) : i = 1, …, n}  (4)

y_{i} is assumed to be a noisy realisation of an underlying latent function f = f(x), that is,

y_{i} = f(x_{i}) + ϵ_{i}, ϵ_{i} ∼ N(0, σ_{ϵ}^{2})  (5)

where x_{i} ∈ ℝ^{q} and y_{i} ∈ ℝ.

The interest here is in estimating the function f, which generally is nonlinear. The GP provides a natural way to make inferences about these functions (Rasmussen and Williams, (2006), [23]). The GP prior over functions can be seen as an extension of the multivariate Gaussian distribution. By definition, a stochastic process is a set of random variables {f(x) : x ∈ X }, indexed by a set, X. A GP is a stochastic process such that for any finite set of function evaluations, f(x) = (f(x_{1}), …, f (x_{n}))^{T}, where f is multivariate Gaussian distributed. We say that for any finite set of elements drawn from X, f is a GP described by a mean, m(.), and covariance function, K(., .), which we write as:

f(x) ∼ GP (m(x), K (x, x′))  (6)

where

m(x) = 𝔼 [f(x)] , K (x_{i}, x_{j}) = 𝔼 {[f(x_{i}) − m(x_{i})] [f(x_{j}) − m(x_{j})]} (7)
m : X → ℝ, and K(., .) : X × X → ℝ.  (8)

Let X_{new} be a matrix with on each row a new input point x_{i}^{new}, i = 1, …, n. To sample a function, we first compute the covariances between all inputs in X_{new} = (x_{1}^{new}, …, x_{n}^{new}) and collect these in an n × n matrix:

K (X_{new}, X_{new}) = ( \begin{matrix}k (x_{1}^{new}, x_{1}^{new}) & … & k (x_{1}^{new}, x_{n}^{new}) \\ k (x_{2}^{new}, x_{1}^{new}) & … & k (x_{2}^{new}, (x_{n}^{new}) \\ ⋮ & ⋮ & ⋮ \\ k (x_{n}^{new}, x_{1}^{new}) & … & k (x_{n}^{new}, x_{n}^{new}) \end{matrix} ) .

Let

f_{new} = (f (x_{1}^{new}) , …, f (x_{n}^{new}))^{T} .

The joint distributionp(f, f_{new}) is given by:

( \begin{matrix}f \\ f_{new} \end{matrix} ) ∼ N [( \begin{matrix}m(X) \\ m(X_{new}) \end{matrix} ) , ( \begin{matrix}K (X, X) & K (X, X_{new}) \\ K (X_{new}, X) & K (X_{new}, X_{new}) \end{matrix} )]

where K(X, X), represents the kernel evaluated at X, K(X_{new}, X_{new}) is the covariance matrix between the new points, K (X, X_{new}) is the covariance matrix between the observed points and the new values, and K(X_{new}, X) is the covariance matrix between the new and the observed points. The covariance function must be positive definite, that is

∑_{j} ∑_{i} v_{i}K(x_{i}, x_{j})v_{j} ≥ 0, for all v_{i}, x_{i}.  (9)

To complete the prior specification, we need to specify the mean and covariance functions. Any positive definite covariance function can be chosen; for K(., .) are the following, see Table 1.

Table 1: Covariance functions.

Kernel  K (x, x′)
Exponential  σ^{2} exp (−\frac{|x−x′|}{l})
Squared Exponential  σ^{2} exp (−\frac{1}{2} \frac{|x−x′|^{2}}{l^{2}} )
Matérn  σ^{2} \frac{2^{1−ν}}{Γ(ν)} (\frac{|x−x′|}{l})^{ν} κ_{ν} (\frac{|x−x′|}{l})
Brownian Motion  min (x, x′)

The distribution predicted by the GP can be determined by the conditional rules of multivariate Gaussian distribution:

f_{new}|f, X, y ∼ N (𝔼 (f_{new}) , ℂov (f_{new}))  (10)

where

𝔼 (f_{new}) = m (X_{new}) + K (X_{new}, X) K^{−1} (X, X) [f − m (X)]  (11)

and

ℂov (f_{new}) = K (X_{new}, X_{new}) − K (X_{new}, X) K^{−1} (X, X) K (X, X_{new}) .

For a Gaussian likelihood:

y_{i} = f(x_{i}) + ϵ_{i} , ϵ_{i} ∼ N (0, σ_{ϵ}^{2}) , y|f ∼ N (f(x_{i}), σ_{ϵ}^{2}I) (12)

where I is the identity matrix. The noise can be included in the covariance function, as follows:

K (f(x_{i}), f(x_{j})) = K (x_{i}, x_{j}) + δ_{ij}σ_{ϵ}^{2}  (13)

where δ_{ij} is the Kronecker delta, and σ_{ϵ}^{2} is the noise variance between layers.

The uncertainty is now present in the observations, and the joint distribution over the unknown data and the known data is augmented in the covariance equation by

( \begin{matrix}f \\ f_{new} \end{matrix} ) ∼N (( \begin{matrix}m (X) \\ m(X_{new}) \end{matrix} ); ( \begin{matrix}K (X, X) + σ_{ϵ}^{2}I & K (X, X_{new}) \\ K (X_{new}, X) & K (X_{new}, X_{new}) \end{matrix} )) .

The marginal distribution is given by

f_{new}|X_{new}, X, f ∼ N (𝔼 (y_{new}) , ℂov (y_{new}))  (14)

where

𝔼 (y_{new}) = m (X_{new}) + K (X_{new}, X) (K (X, X) + σ_{ϵ}^{2}I)^{−1} [y − m (X)]

and

ℂov (y_{new}) = K (X_{new}, X_{new}) − K (X_{new}, X) (K (X, X) + σ_{ϵ}^{2}I)^{−1} K (X, X_{new}) .

In this case the integrals required to infer a posterior, p (f_{new}|X_{new}, X, f), are tractable.

4 Deep neural networks and Gaussian processes

In this article, we consider fully-connected ANNs with layers numbered from l = 0 (input) to l = L − 1 (output), each containing, N_{0}, …, N_{L−1} neurons, and with a Lipschitz, twice differentiable nonlinearity activation function ϕ : ℝ → ℝ, with bounded second derivative. For each x ∈ ℝ^{d_{in}} denote the input to the network (x = (x_{1}, …, x_{d_{in}})), and z^{(l)} ∈ ℝ^{d_{out}} denote its output. We use z_{i}^{(l)}(x), h_{i}^{(l)}(x) to represent the pre- and post-activation functions at layer l with input x, also, let h_{i}^{(0)} = x. The parameters consist of the connection matrices W_{ij}^{(l)} ∈ ℝ^{N_{l}×N_{l+1}} and bias vectors b_{i}^{(l)} ∈ ℝ^{N_{l+1}} for l = 0, …, L − 1; which are

independent and randomly selected, with zero mean and variances \frac{σ_{w}^{2}}{N_{l}}, and σ_{b}^{2}, respectively.

Now we are going to establish the relationship between a single-hidden layer neural networks, and Gaussian processes. Suppose that z_{j}^{(l)}(x) is a Gaussian process with mean and covariance functions μ^{(l)}(x), K^{(l)}(x, x′), respectively, that is,

z_{j}^{(l)}(x) = (\begin{matrix}z_{j}^{(l)}(x) \\ z_{j}^{(l)}(x′) \end{matrix}) ∼ GP (μ^{(l)}(x), K^{(l)} (x, x′))

where

μ^{(l)}(x) = 𝔼 \{z_{j}^{(l)}(x)\} = (\begin{matrix}μ(x) \\ μ(x′) \end{matrix})

and

K =K^{(l)} (x, x′) = ℂov \{z_{j}^{(l)}(x), z_{j}^{(l)}(x′)\}=(\begin{matrix}K^{(l−1)}(x, x) & K^{(l−1)}(x, x′) \\ K^{(l−1)}(x′, x) & K^{(l−1)}(x′, x′) \end{matrix}) .

The i−th component of the network output, z_{i}^{(1)}, is computed as:

z_{i}^{(1)}(x) = b_{i}^{(1)}+∑_{j=1}^{N_{1}} W_{ij}^{(1)}h_{j}^{(1)}(x), h_{j}^{(1)}(x) = ϕ (b_{j}^{(0)} + ∑_{k=1}^{d_{in}} W_{ij}^{(0)}x_{k}) . (15)

Note that there is a dependency on the input data vector x, and also since the weight and bias are considered i.i.d., the post-activations h_{j}^{(1)} and h_{k}^{(1)} are in-

dependent for j ≠ k. Also, z_{i}^{(1)}(x) is a sum of i.i.d terms, it follows from the Central Limit Theorem that in the limit of infinite width N_{1} → ∞, z_{i}^{(1)}(x) will be Gaussian distributed, Lee, et. al (2018) [14].

Suppose that z_{j}^{(l)}(x) is a Gaussian processes, i.i.d for every j, and that

h_{j}^{(l)}(x), also are independent and identically distributed. Then after l − 1 steps, the recurrence relation for a feedforward network is defined as

z_{i}^{(l)}(x) = b_{i}^{(l)} + ∑_{j=1}^{N_{l}} W_{ij}^{(l)}h_{j}^{(l)}(x), h_{j}^{(l)}(x) = ϕ (z_{j}^{(l−1)}(x)) . (16)

Prior on weights:

b_{j}^{(l)}|σ_{b}^{2} ∼ N (0, σ_{b}^{2}) , W_{ij}^{(l)}|Σ_{W} ∼ N (0, Σ_{W}) .  (17)

Prior on hyperparameters:

σ_{b}^{2}|α, β ∼ IG (α, β) Σ_{W} |ν, R ∼ IG (ν, R) .  (18)

298 k. román — a. cumbicus — s. infante — r. fonseca-delgado Note, z_{i}^{(l)}(x) is a sum of i.i.d. random terms so that, as N_{l} → ∞, any finite collection {z_{i}^{(l)}(x), l = 1, …, L} will have joint multivariate Gaussian distribution, i.e,

z_{i}^{(l)}(x) ∼ GP (μ^{(l)}(x), K^{(l)} (x, x′)), z_{i}^{(l)}(x) = (z_{i}^{(l)}(x_{1}), …, z_{i}^{(l)}(x_{n})) .(19) A general equation is now established to approximate the covariance for a bivariate Gaussian process:

K^{(l)} (x, x′) = 𝔼 {[z_{i}^{(l)}(x) − 𝔼 (z_{i}^{(l)}(x))] [z_{i}^{(l)}(x′) − 𝔼 (z_{i}^{(l)}(x′))]}

= 𝔼 {[b_{i}^{(l)} + ∑_{j=1}^{N_{l}} W_{ij}^{(l)}h_{j}^{(l)}(x)] [b_{i}^{(l)} + ∑_{j=1}^{N_{l}} W_{ij}^{(l)}h_{j}^{(l)}(x′)]}

= σ_{b}^{2} + Σ_{W} ∑_{j=1}^{N_{l}} 𝔼 [h_{j}^{(l)}(x)h_{j}^{(l)}(x′)]

= σ_{b}^{2} + Σ_{W} ∑_{j=1}^{N_{l}} 𝔼 [ϕ (z_{j}^{(l−1)}(x)) ϕ (z_{j}^{(l−1)}(x′))]  (20)

where the calculation of the expectation is a two dimensions Gaussian integral:

𝔼 (ϕ (z_{j}^{(l−1)}(x)) ϕ (z_{j}^{(l−1)}(x′))) = ∫_{−∞}^{∞} ∫_{−∞}^{∞} ϕ (z_{j}^{(l−1)}(x)) ϕ (z_{j}^{(l−1)}(x′))
× p (z_{j}^{(l−1)}(x), z_{j}^{(l−1)}(x′))
× dz_{j}^{(l−1)}(x)dz_{j}^{(l−1)}(x′).  (21)

Since:

(\begin{matrix}z_{j}^{(l)}(x) \\ z_{j}^{(l)} (x′) \end{matrix}) ∼ GP (μ^{(l)}(x), K^{l} (x, x′))  (22)

then

∫_{−∞}^{∞} ∫_{−∞}^{∞} ϕ (z_{j}^{(l−1)}(x)) ϕ (z_{j}^{(l−1)}(x′)) \frac{1}{2π|K|^{\frac{1}{2}}}

× exp \{−\frac{1}{2}(z_{j}(x)^{(l−1)} − μ(x), z_{j}^{(l−1)}(x′) − μ(x′)) K^{−1} (\begin{matrix}z_{j}(x)^{(l−1)} − μ(x) \\ z_{j}^{(l−1)}(x′) − μ(x′) \end{matrix})\}

× dz_{j}^{(l−1)}(x)dz_{j}^{(l−1)}(x′).

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(2): 289–312, Jul–Dec 2022

deep gaussian processes and infinite neural networks... 299 Let

K^{−1} = (\sqrt{K}^{T} )^{−1} (\sqrt{K})^{−1} .  (23)

Consider the following transformation:

(\begin{matrix}ξ_{i_{1}} \\ ξ_{i_{2}} \end{matrix}) = (\sqrt{K})^{−1} (\begin{matrix}z_{j}(x) − μ(x) \\ z_{j}(x′) − μ(x′) \end{matrix}) ⇒ (\sqrt{K}) (\begin{matrix}ξ_{i_{1}} \\ ξ_{i_{2}} \end{matrix}) = (\begin{matrix}z_{j}(x) − μ(x) \\ z_{j}(x′) − μ(x′) \end{matrix}) .

To simplify the notation, the following variable change is made:

z_{j}(x) = z_{j}(x)^{(l−1)}, z_{j}(x′) = z_{j}(x′)^{(l−1)}.  (24)
\sqrt{K}ξ_{i_{1}} = z_{j}(x) − μ(x) ⇒ z_{j}(x) = \sqrt{K}ξ_{i_{1}} + μ(x)

⟺ ξ_{i_{1}} = \frac{z_{j}(x) − μ(x)}{\sqrt{K}} ∼ N(0, 1)

and

\sqrt{K}ξ_{i_{2}} = z_{j}(x′) − μ(x′) ⇒ z_{j}(x′) = \sqrt{K}v + μ(x′)

⟺ ξ_{i_{2}} = \frac{z_{j}(x′) − μ(x′)}{\sqrt{K}} ∼ N(0, 1).

The Jacobian of the transformation is:

J = (\begin{matrix}\frac{∂z_{j}(x)}{∂ξ_{i_{1}}} & \frac{∂z_{j}(x)}{∂ξ_{i_{2}}} \\ \frac{∂z_{j}(x′)}{∂ξ_{i_{2}}} & \frac{∂z_{j}(x′)}{∂ξ_{i_{2}}} \end{matrix} ) = (\begin{matrix}\sqrt{K} & 0 \\ 0 & \sqrt{K} \end{matrix}) ⇒ |J| = |K|^{\frac{1}{2}}  (25)

then

\frac{1}{2π|K|^{\frac{1}{2}}} ∫_{−∞}^{∞} ∫_{−∞}^{∞} ϕ (\sqrt{K}ξ_{i_{1}} + μ(x)) ϕ (\sqrt{K}ξ_{i_{2}} + μ(x′))

× \frac{1}{2π} exp (−\frac{1}{2} (ξ_{i_{1}}^{2} + ξ_{i_{2}}^{2})) |J|dξ_{i_{1}}dξ_{i_{2}}

= \frac{1}{\sqrt{2π}} \frac{1}{\sqrt{2π}} ∫_{−∞}^{∞} ∫_{−∞}^{∞} ϕ (\sqrt{K}ξ_{i_{1}} + μ(x)) ϕ (\sqrt{K}ξ_{i_{2}} + μ(x′))

× exp (−\frac{1}{2}ξ_{i_{1}}^{2}) exp (−\frac{1}{2}ξ_{i_{2}}^{2}) dξ_{i_{1}}dξ_{i_{2}}

= \frac{1}{\sqrt{2π}} ∫_{−∞}^{∞} ϕ (\sqrt{K}ξ_{i_{1}} + μ(x)) exp (−\frac{1}{2}ξ_{i_{1}}^{2}) dξ_{i_{1}}

× \frac{1}{\sqrt{2π}} ∫_{−∞}^{∞} ϕ (\sqrt{K}ξ_{i_{2}} + μ(x′)) exp (−\frac{1}{2}ξ_{i_{2}}^{2}) dξ_{i_{2}}.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(2): 289–312, Jul–Dec 2022

The approximation of the last integral depends on the choice of the sigmoid function ϕ(.).

In the case of deep neural networks, the transfer function ϕ(z) is a bounded function where all moments are bounded. Then we can apply Central Limit Theorem to show that the stochastic process is a Gaussian process, Williams (1996) [27] and Lee et. al. (2018) [14].

Cho and Saul (2009) [5] developed a new family of covariance functions which allows computing the correlation between two vectors x, z ∈ ℝ^{N_{l}}. They define the n−th order arc-cosine kernel function via the integral representation:

K^{(l)}(x, z)=2 ∫ \frac{1}{(2π)^{\frac{N_{l}}{2}}} exp (−\frac{∥w∥^{2}}{2})Θ(w .x)Θ (w .z) (w .x)^{l} (w .z)^{l} dw, (26)

where Θ (t) = \frac{1}{2} (1 + sign(t)) denote the Heaviside step function. The integral representation (26) allows the kernel of covariance functions to be positive definite, and that the dependence between x y z, can be written as:

K^{(l)} (x, z) = \frac{1}{π} ∥x∥^{l}∥z∥^{l}J_{l}(θ)  (27)

where all the angular dependence is captured by J_{l}(θ). The angular dependence is given by:

J_{l}(θ) = (−1)^{l} (sin θ)^{2l+1} (\frac{1}{sin θ} \frac{∂}{∂θ} )^{l} (\frac{π − θ}{sin θ} ) .  (28)

In particular, when l = 0, we have the angle between x and z. When l > 0, the angular dependence is more complicated. Some terms of J_{l}(θ), are shown:

J^{0}(0) = θ, J^{1}(θ) = sin θ + (π − θ) cos θ,
J^{2}(θ) = 3 sin θ cos θ + (π − θ) (1 + 2 cos^{2} θ) .

The arc-cosine kernel for l = 0, is represented by:

K^{(0)} = 1 − arccos^{−1} (\frac{xz}{∥x∥∥z∥}) .  (29)

Neural network models are strongly related to the kernel of functions defined in (26), when considering the inner product between the different outputs of the neural network as:

ϕ(x)ϕ(z) = ∑_{i=1}^{m} Θ(w_{i} .x)Θ (w_{i} .z) (w_{i} .x)^{l} (w_{i} .z)^{l}  (30)

where w_{i} denote i−th row of the weight matrix W and m is the number of output units. In the limit, it can be seen that the equation (30) is equivalent to (26), Cho and Saul (2009), [5]:

lim_{m→∞} \frac{2}{m}ϕ(x)ϕ(z) = K^{(l)}(x, z).  (31)

Also, Cho and Saul (2009) [5] proved that for a vector of inputs x = (1, x_{1}, …, x_{N_{l}}), their characteristics can be mapped by means of a nonlinear transformation ϕ(x), using kernel functions:

K^{(l)}(x, z) = ϕ (ϕ (… ϕ(x)) .ϕ (… ϕ(z))) .  (32)

The iterated equation (32) mimics a multilayer neural network, for example for a one-layer neural network, K (x, z) = ϕ(x)ϕ(z). Cho and Saul (2009) [5] define a recursive kernel through a new mapping of features through compositions such as ϕ (ϕ(x)). In the case of a linear kernel K (x, z) = xz, the composition is ϕ (ϕ(x)) = ϕ(x) = x, and for homogeneous polynomial kernels K (x, z) = (xz)^{N_{l}}; the composition is:

K(x, z) = ϕ (ϕ(x)) .ϕ (ϕ(z)) = (ϕ(x)ϕ(z))^{N_{l}} = ((x.z)^{N_{l}})^{N_{l}} = (x.z)^{N_{l}^{2}} . (33)

Then we consider the composition of l layers based on the iterate defined in (32), applying a mathematical induction the inductive step is given by:

K^{(l)}(x, z) = \frac{1}{π} [K^{(l−1)}(x, x)K^{(l−1)}(z, z)]^{\frac{l}{2}} J_{l} (θ^{(l)})  (34)

where θ^{(l)} is the angle between x and y in the feature space:

θ^{(l)} = cos^{−1} \{\frac{K^{(l)}(x, z)}{[K^{(l)}(x, x)K^{(l)}(z, z)]^{\frac{1}{2}}} \} .  (35)

Recently, a linear rectified function ReLU(t) = max(0; t), was successfully used in neural networks as it carries the neuron signal better. ReLUs have the desirable property that they do not require input normalization, Krizhevsky et al. (2012) [12]. To compute the given integral (21), a rectified linear sigmoidal can be used, Hazan et al. (2015) [9], which results in an analytical kernel given by

K_{ReLU} (x_{i}, x_{j}) = \frac{∥x_{i}∥∥x_{j}∥}{π} sin (cos^{−1} (ρ_{ij})) + (π − cos^{−1} (ρ_{ij})) ρ_{ij}

where

ρ_{ij} = \frac{⟨x_{i}, x_{j}⟩}{∥x_{i}∥∥x_{j}∥}, ⟨x_{i}, x_{j}⟩ = ∫ ∫ ϕ(x_{i})ϕ(x_{j})p (x_{i}, x_{j}) dx_{i}dx_{j}.

To compute the entries of K(., .), let ϕ(t) = ReLU(t),

(\begin{matrix}z′_{1} \\ z′_{2} \end{matrix}) ∼ N (0, K) , 0 = (\begin{matrix}0 \\ 0 \end{matrix}) K = (\begin{matrix}σ_{1}^{2} & ρσ_{1}σ_{2} \\ ρσ_{1}σ_{2} & σ_{2}^{2} \end{matrix})

then

𝔼_{z′_{1},z′_{2}} (ϕ (z′_{1}) ϕ (z′_{2})) = h (σ_{1}, σ_{2}, ρ)
= \frac{σ_{1}σ_{2}}{π} sin (cos^{−1}(ρ)) + ρ (π − cos^{−1}(ρ))

and K_{ReLU}^{(1)} is given by

(\begin{matrix}h (\sqrt{1+α}∥x_{i}∥,√ ; \overline{1 +α}∥x_{i}∥, \frac{α}{1+α} ) & h (\sqrt{1 + α}∥x_{i}∥,\sqrt{1+α}∥x_{j}∥, \frac{α}{1+α} \frac{⟨x_{i},x_{j} ⟩}{∥x_{i}∥∥x_{j} ∥} ) \\ h (\sqrt{1+α}∥x_{i}∥, \sqrt{1+α}∥x_{j}∥, \frac{α}{1+α} \frac{⟨x_{i},x_{j} ⟩}{∥x_{i}∥∥x_{j} ∥} ) & h (\sqrt{1+α}∥x_{j}∥,\sqrt{1+α}∥x_{j}∥, \frac{α}{1+α} ) \end{matrix}).

Iterating successively, we obtain:

K_{ReLU}^{(2)} (x_{i}, x_{j}) = 𝔼_{z_{1},z_{2}} (ϕ (z_{1}) ϕ (z_{2})) , (\begin{matrix}z_{1} \\ z_{2} \end{matrix}) ∼ N (0, K)

is a recursive equation of h with the appropriate parameters, see Hazan et al. (2015) [9].

On the other hand, using GP prior functions allows Bayesian inference to be made precisely to obtain predictions and estimate the uncertainty in deep neural network models. The estimation of the parameters is not required through training based on a gradient-type algorithm.

Under the Bayesian statistical approach, the weights and biases of the network are generated following a probability distribution p(W|θ), where (W = (W, b), represents the weights and bias, and θ = (α, β, ν, R) ∼ p(θ) represents the hyperparameters, which can be integrated.

p (W) = ∫ p(W|θ)p(θ)dθ.  (36)

Given a training sample D = {(x_{1}, y_{1}), (x_{2}, y_{2}), …, (x_{n}, y_{n})}, where x = (x_{1}, …, x_{n})^{T} and y = (y_{1}, …, y_{n})^{T} denote the input and output of the network, respectively. Suppose now that we want to make a prediction with a test data x_{new} using priors over funtions rather than weights z(x) = (z_{1}, …, z_{n}) restricted to input values x. Then

p (z_{new}|D, x_{new}) = ∫ p (z_{new}|z, x, x_{new}) p (z|D) dz
= \frac{1}{p(y)} ∫ p (z_{new}, z|x_{new}, x) p (y|z) dz  (37)

where y = (y_{1}, …, y_{n})^{T} are the targets on the training set, p(y) is a marginal likelihood, and p (y|z) corresponds to observation noise. We will assume a noise consisting of a Gaussian (y|z ∼ N(0, σ_{ϵ}^{2})).

The importance of choosing priors over functions implies that z_{1}, …, z_{n}, z_{new} are generated from a Gaussian processes

( \begin{matrix}z \\ z_{new} \end{matrix} ) ∼ GP [( \begin{matrix}0 \\ 0 \end{matrix} ) , ( \begin{matrix}K (D, D) & K (x_{new}, D) \\ K (D, x_{new}) & K (x_{new}, x_{new}) \end{matrix} )] .

Then the integral in (37) can be obtained exactly, since the marginal distribution is:

z_{new}|D, x_{new} ∼ N (μ_{post}, K_{post})  (38)

where

μ_{post} = K (x_{new}, D) (K (D, D) + σ_{ϵ}^{2}I_{n})^{−1} y

and

K_{post} = K (x_{new}, x_{new}) − K (x_{new}, D) (K (D, D) + σ_{ϵ}^{2}I_{n})^{−1} K^{T} (x_{new}, D)

where I_{n} is the n × n matrix identity. The predicted distribution of z_{new}|D, x is clearly determined. Deep neural network training works using a Bayesian approach. The covariance function is determined by choosing a prior Gaussian process. In this case, the model depends on the depth, non-linearity of the transfer function, weights and biases.

Finally, the results obtained through the estimation method by a neural network approximated by a Gaussian process are compared with those obtained by a neural network. The weights are optimized using Adam’s variant of stochastic gradient descent (SGD) algorithm, using a loss function MSE, Lee et. al. (2018) [14].

5 Experimental analysis

The experimental goal is to validate the performance of the proposed method in diagnosing Alzheimer’s disease. Different experiments were executed, varying the size of the input set, to view the accuracy difference when the input size increased. As a reference point, this section compares the results of the proposed method with multi-layer neuronal networks.

5.1 Data description

The data set was used by a previous investigation of Alzheimer’s disease (AD) detection. AD is a progressive and irreversible brain disorder that causes memory problems, slowly destroying it until losing the ability to develop simple tasks. A non-invasive method to study AD uses Electroencephalograms (EEGs) that register the brain’s electrical activity. However, the EEG raw signals are challenging to classify directly. Martínez et . (2021) proposed to extract five significant features provided by public EGG experiments from https:// osf. io/ jbysn/ . In addition, for this study, we only considered the five critical features mentioned in Table 2.

Table 2: Features extracted from the EEGs to identify the AD effects.

Notation  Feature Name
HFD Higuchi Fractal Dimension
LZC  Lempel Ziv Complexity
PSD  Power Spectral Density
RP  Relative Power
SE  Sample Entropy

Each data set has 16 dimensions corresponding to each channel in an EGG test registration. However, for Relative Power feature (plots in section 5.3), the input size is 64. Also, there is an extra column for the label class in all files.

In this experiment, we cleaned up the data of significant features by randomly removing data to leave the same number of labels for each class in the data set. In total, 514 observations were considered, of which 85% were for training data and the rest 15% for testing.

5.2 Experiment description

The main goal is to use the above data set to classify the positive and negative cases of Alzheimer’s disease. We use the classification models in Table 3.

Figure 2: Experiment description diagram.

Each file was used to train the GP and NN models independently. That is, in total, results are obtained from 7 classification models (GP, NN = 5, NN = 10, NN = 50, NN = 100, NN = 500, and NN = 1000). Comparing the performance of GPs against neural networks changes parameters such as network width, data set size and activation function. In each model, the training data size for the learning process was changed; that is, n = 75, n = 100, n = 250 and n = 436 observations are taken to train the model for each of the data-set sizes. Similarly, all models are tested with the ReLu and Tanh activation functions. Also, in neural networks, the width of the network is changed to compare its performance against the Gaussian process. In Figure 2 there is a diagram that explains the experiment for each data set.

Lastly, in the following repository there is a Python code implementation of

this work, and the data sets used, https:// github. com/ KrishnaRoman/
Deep-Gaussian-Processes

Table 3: Notation of the classification models used in this work.

Notation Model Name
NN-5  Neural Network with width 5
NN-10  Neural Network with width 10
NN-50  Neural Network with width 50
NN-100 Neural Network with width 100
NN-500 Neural Network with width 500
NN-1000 Neural Network with width 1000

5.3 Performance results

The results for the Relative Power feature are shown in Figure 3. Clearly, the model that uses the proposed Gaussian Process has the best test accuracy. This figure shows that the Neural Network models behave better with Tanh activation function than Relu on this feature. While in the case of the GP there is no significant difference. Additionally, the performance results on the other features are shown in Table 4.

Figure 3: Relative power performance results.

5.4 Error results

The predictive performance of Deep Gaussian Process (DGP) and Neural Networks (NN) models are evaluated using the following metric, Mean Squared Error (MSE), defined by:

MSE = \frac{1}{n} ∑_{i=1}^{n} (y_{i} − ŷ_{i})^{2}  (39)

where ŷ_{i} indicates the predicted value, y_{i} indicates the actual values, and n is the number of predictions. The metric range is in [0, ∞) and lower values indicate better performance. It can be seen in Figure 4 for Relative Power that in terms of error, the GP with the Relu activation function is the best, whereas the Neural Network with 5-width is better in Tanh. Also, the error results on the other features are shown in Table 4.

Figure 4: Relative power error results.

The metrics for all features and models were recorded in Table 4. In this case only was considered the results with the highest dataset size that is 436.

6 Discussion and conclusions

Using a nonstochastic gradient-based training, we use EEG Signals for diagnosing Alzheimer’s Disease to test our GP’s behavior, which uses a hierarchical Bayesian structure to model the weights and biases of a neural network, and compare it with neural networks varying their widths. A general formula was derived to evaluate the resulting integrals of Gaussian processes with non-linear

Table 4: Performance and error results for all features and models.

Dataset Model  Test  Test Model  Test  Test
(ReLU) accuracy MSE (tanh)  accuracy MSE
GP  0.935 0.0564 GP  0.948 0.1092
NN-10  0.935 0.0508 NN-10  0.883 0.0784
LZC  NN-1000 0.961 0.0404 NN-1000 0.948 0.059
GP  0.987 0.0215 GP  0.974 0.0221
NN-10  0.935 0.0414 NN-10  0.961 0.0248
RP  NN-1000 0.987 0.0243 NN-1000 0.974 0.0299
GP  1.000 0.0125 GP  0.987 0.0153
NN-10  0.961 0.0243 NN-10  1.000 0.0152
HFD  NN-1000 1.000 0.0114 NN-1000 1.000 0.0178
GP  0.974 0.0246 GP  0.987 0.0245
NN-10  0.961 0.0244 NN-10  0.961 0.0215
PSD  NN-1000 0.987 0.0157 NN-1000 0.974 0.0207
GP  1.000 0.0125 GP  0.987 0.0153
NN-10  0.961 0.0243 NN-10  1.000 0.0152
SE  NN-1000 1.000 0.0114 NN-1000 1.000 0.0178

transfer functions, and obtained a kernel to update the covariance functions. The proposed methodology was applied to the classification of EEG signals for diagnosing Alzheimer’s disease, considering five data sets and estimating the models varying size of the samples. This study showed that DGP can be used in supervised learning and classification tasks. As mentioned in Lee et al. (2018), we check that the GP behaves as a neural network with an infinite number of neurons in the hidden layers. In general, all test accuracies were very high because the data set consisted on a limited number of patients and the work done by Martínez-Arias (2020) [17] on feature extraction and selection.

In our experiments, we realize that inverse gamma distribution hyperparameters fixed well for the model learning process. The bias and weight variances are randomly chosen using these hyper-parameters, and it contributes to evaluating the model performance. The results show that GP classifies well and sometimes even better than neural network models. In this sense, we proved that Deep Gaussian Model working with regression data is a Gaussian Process when the neural network width tends to be infinite. Also, this model is more effective when the data size is more extensive. Another big difference between DGP against traditional neural network models is the use of GP properties like precise uncertainty estimates.

For future research, we expect to prove our DGP model to classify other data types, such as crude data or images. Crude ECG signal data could be an interesting scenario to use DGP model to organize positive cases of Alzheimer’s disease. Also, analyze what role distribution hyper-parameters play in the model with other data types.Finally, the integral I_{z} computed in section 4 may be used to estimate the kernel using a different method that we proposed. This could be another research to see which alternative method works better and to compare the GP performance results.

Acknowledgements

We want to express our special thanks to the journal’s anonymous reviewers for their suggestions to improve the manuscript. In the same way, we thank the authors of the previous works who provided us with the data for our research.

Funding

We also thank the project Functional Data Analysis: Methods and Applications, with registration REG-INV-19-04054, Yachay Tech.

References

[1] D. Abásolo, J. Escudero, R. Hornero, C. Gómez, P. Espino, Approximate entropy and auto mutual information analysis of the electroencephalogram in Alzheimer’s disease patients, Medical & Biological Engineering & Computing 46 (2008), 1019–1028. Doi: 10.1007/s11517-008-0392-1

[2] D. Agrawal, T. Papamarkou, J. Hinkle, Wide neural networks with bottlenecks are deep Gaussian processes, Journal of Machine Learning Research 21 (2020), no. 175, 1–66. Available from: Link

[3] N. Cedeño, G. Carillo, M.J. Ayala, S. Lalvay, S. Infante, Analysis of chaos and predicting the price of crude oil in Ecuador using deep learning models, in: T. Guarda, F. Portela & M.F. Santos (Eds.) Advanced Research in Technologies, Information, Innovation and Sustainability, part of Communications in Computer and Information Science series 1485, Springer, Cham, 2021, 318–332. Doi: 10.1007/978-3-030-90241-4_25 [4] N. Cedeño, S. Infante, Estimation of ordinary differential equations

solutions with Gaussian processes and polynomial chaos expansion, in: J.P. Salgado Guerrero, J. Chicaiza Espinosa, M. Cerrada Lozada, S. Berrezueta-Guzman (Eds.) Information and Communication Technologies, part of Communications in Computer and Information Science series 1456, Springer, Cham, 2021, pp. 3–17. Doi: : 10.1007/978-3-030-89941-7_1

[5] Y. Cho, L. Saul, L. Kernel methods for deep learning, in: Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams & A. Culotta (Eds.) Advances in Neural Information Processing Systems 22 (2009), 342–350. Available from: Link

[6] A. Damianou, N. Lawrence, Deep Gaussian processes, Proceedings of the 16th International Conference on Artificial Intelligence and Statistics, Scottsdale AZ, U.S.A. (2013), 207–215. Available from: Link

[7] A. Garriga-Alonso, C.E. Rasmussen, L. Aitchison, Deep convolutional networks as shallow Gaussian procesesses, arXiv, 2019. 1808.05587 [stat.ML]. Doi: 10.48550/arXiv.1808.05587

[8] I. Goodfellow, Y. Bengio, A. Courville (2016). Deep Learning. The MIT Press, Cambridge MA, U.S.A. http://www.deeplearningbook.org

[9] T. Hazan, T. Jaakkola, Steps toward deep kernel methods from infinite neural networks, arXiv, 2015. 1508.05133 [cs.LG]. Doi: 10.48550/arXiv.1508.05133

[10] S. Infante, J. Ortega, F. Cedeño, Estimación de datos faltantes en estaciones meteorológicas de Venezuela vía un modelo de redes neuronales, Revista de Climatología 8, 51–70. Available from: Link

[11] M.E. Khan, A. Immer, E. Abedi, M. Korzepa, Approximate inference turns deep networks into Gaussian processes, arXiv, 2019. 1906.01930 [stat.ML]. Doi: 10.48550/arXiv.1906.01930

[12] A. Krizhevsky, I. Sutskever, G. Hinton, Image net classification with deep convolutional neural networks, in: F. Pereira, C.J. Burges, L. Bottou & K.Q. Weinberger (Eds.) Advances in Neural Information Processing Systems 25 (2012), 1097–1105. Available from: Link

[13] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (2015), no. 7553, 436–444. Doi: 10.1038/nature14539 [14] J. Lee, Y. Bahri, R. Novak, S.S. Schoenholz, J. Pennington, J. Sohl-

Dickstein, Deep neural networks as Gaussian processes, arXiv, 2018. 1711.00165v3 [stat.ML]. Doi: 10.48550/arXiv.1711.00165

[15] M. Luca, G. Barlacchi, B. Lepri, L. Pappalardo, Deep learning for human mobility: a survey on data and models, arXiv, 2020. 2012.02825 [cs.LG]. Doi: 10.48550/arXiv.2012.02825

[16] D.J.C. MacKay, Probable networks and plausible predictions. A review of practical Bayesian methods for supervised neural networks, Network: Computation in Neural Systems, 6 (1992), no. 3, 469–505. Doi: 10.1088/0954-898X_6_3_011

[17] P. Martínez-Arias, R. Fonseca-Delgado, R. Salum, I. Amaro-Martín, Alzheimer’s disease diagnosis system using electroencephalograms and machine learning models, Iberian Journal of Information Systems and Technologies. (2020) pp. 275–288. Available from: Link

[18] A.G. Matthews, M. Rowland, J. Hron, R.E. Turner, Z. Ghahramani, Gaussian process behaviour in wide deep neural networks, arXiv, 2018. 1804.11271v2 [stat.ML]. Doi: 10.48550/arXiv.1804.11271

[19] D. Ming, D. Williamson, S. Guillas, Deep Gaussian process emulation using stochastic imputation, arXiv, 2021. 2107.01590 [stat.ML] Doi: 10.48550/arXiv.2107.01590

[20] A. Mosavi, S. Ardabili, A.R. Varkonyi-Koczy, (2020) List of deep learning models, in: A.R. Várkonyi-Kóczy (Ed.) Engineering for Sustainable Future, Springer Nature, Cham Switzerland, 2020, pp. 202–214. Doi: 10.1007/978-3-030-36841-8_20

[21] R. Neal, Bayesian Learning for Neural Networks, Lecture Notes in Statistics, Springer Verlag, New York, 1996. Doi: 10.1007/978-1-4612-0745-0

[22] R. Novak, L. Xiao, J. Lee, Y. Bahri, G. Yang, J. Hron, …, J. Sohl-Dickstein, Bayesian deep convolutional networks with many channels are Gaussian processes, arXiv, 2019. 1810.05148v4 [stat.ML]. Doi: 10.48550/arXiv.1810.05148

[23] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, the MIT Press, Cambridge MA, U.S.A., 2006. Available from: http://gaussianprocess.org/gpml/ [24] H. Salimbeni, M. Deisenroth, Doubly stochastic variational inference

for deep Gaussian processes, arXiv, 2017. 1705.08933v2 [stat.ML]. Doi: 10.48550/arXiv.1705.08933

[25] A. Sauer, R.B. Gramacy, D. Higdon, Active learning for deep Gaussian process surrogates, arXiv, 2021. 2012.08015 [stat.ME]. Accepted in Technometrics. Doi: 10.48550/arXiv.2012.08015

[26] S. Schoenholz, J. Gilmer, S. Ganguli, J. Sohl-Dickstein, Deep information propagation, arXiv, 2017. 1611.01232v2 [stat.ML]. Doi: 10.48550/arXiv.1611.01232

[27] C.K.I. Williams, Computing with infinite networks, in: M.C. Mozer and M. Jordan & T. Petsche (Eds.) Advances in Neural Information Processing Systems 9 (1997), 295–301. Available from: Link

[28] A.G. Wilson, P. Izmailov, Bayesian deep learning and a probabilistic perspective of generalization, arXiv, 2020. 2002.08791v3 [cs.LG]. Doi: 10.48550/arXiv.2002.08791

[29] I.H. Witten, E. Frank, M.A. Hall, C.J. Pal, Data Mining: Practical Machine Learning Tools and Techniques, 4th Edition. (2016)

[30] Z. Zhao, M. Emzir, S. Särkkä, Deep state-space Gaussian processes Statistics and Computing. (2021), 31–75. Doi: 10.1007/s11222-021-10050-6