{'fraction_non_alphanumeric': 0.04042106773308466, 'fraction_numerical': 0.021752775137636068, 'mean_word_length': 4.693108728943338, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Challenges in Markov chain Monte Carlo for Bayesian neural networks

Markov chain Monte Carlo (MCMC) methods have not been broadly adopted in Bayesian neural networks (BNNs). This paper initially reviews the main challenges in sampling from the parameter posterior of a neural network via MCMC. Such challenges culminate to lack of convergence to the parameter posterior. Nevertheless, this paper shows that a non-converged Markov chain, generated via MCMC sampling from the parameter space of a neural network, can yield via Bayesian marginalization a valuable posterior predictive distribution of the output of the neural network. Classification examples based on multilayer perceptrons showcase highly accurate posterior predictive distributions. The postulate of limited scope for MCMC developments in BNNs is partially valid; an asymptotically exact parameter posterior seems less plausible, yet an accurate posterior predictive distribution is a tenable research avenue.


MOTIVATION
The universal approximation theorem (Cybenko, 1989) and its subsequent extensions (Hornik, 1991;Lu et al., 2017) state that feedforward neural networks with exponentially large width and widthbounded deep neural networks can approximate any continuous function arbitrarily well. This universal approximation capacity of neural networks along with available computing power explain the widespread use of deep learning nowadays.
Bayesian inference for neural networks is typically performed via stochastic Bayesian optimization, stochastic variational inference (Polson and Sokolov, 2017) or ensemble methods (Ashukha et al., 2020;Wilson and Izmailov, 2020). MCMC methods have been explored in the context of neural networks, but have not become part of the Bayesian deep learning toolbox.
Department of Mathematics, The University of Manchester, Manchester, UK, and Computational Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA The slower evolution of MCMC methods for neural networks is partly attributed to the lack of scalability of existing MCMC algorithms for big data and for high-dimensional parameter spaces. Furthermore, additional factors hinder the adaptation of existing MCMC methods in deep learning, including the hierarchical structure of neural networks and the associated covariance between parameters, lack of identifiability arising from weight symmetries, lack of a priori knowledge about the parameter space, and ultimately lack of convergence.
The purpose of this paper is twofold. Initially, a literature review is conducted to identify inferential challenges in MCMC developments for neural networks. Subsequently, Bayesian marginalization based on MCMC samples of neural network parameters is used for attaining accurate posterior predictive distributions of the respective neural network output, despite the lack of convergence of the MCMC samples to the parameter posterior.
An outline of the paper layout follows. Section 2 reviews the inferential challenges arising from the application of MCMC to neural networks. Section 3 provides an overview of the employed inferential framework, including the multilayer perceptron (MLP) model and its likelihood for binary and multiclass classification, the MCMC algorithms for sampling from MLP parameters, the multivariate MCMC diagnostics for assessing convergence and sampling effectiveness, and the Bayesian marginalization for attaining posterior predictive distributions of MLP outputs. Section 4 showcases Bayesian parameter estimation via MCMC and Bayesian predictions via marginalization by fitting different MLPs to four datasets. Section 5 posits predictive inference for neural networks, among else by combining Bayesian marginalization with approximate MCMC sampling or with ensemble training.

PARAMETER INFERENCE CHALLENGES
A literature review of inferential challenges in the application of MCMC methods to neural networks is conducted in this section thematically, with each subsection being focused on a different challenge.

Computational cost
Existing MCMC algorithms do not scale with increasing number of parameters or of data points. For this reason, approximate inference methods, including variational inference (VI), are preferred in highdimensional parameter spaces or in big data problems from a time complexity standpoint (MacKay, 1995;Blei, Kucukelbir and McAuliffe, 2017;Blier and Ollivier, 2018). On the other hand, MCMC methods are better than VI in terms of approximating the log-likelihood (Dupuy and Bach, 2017).
Literature on MCMC methods for neural networks is limited due to associated computational complexity implications. Sequential Monte Carlo and reversible jump MCMC have been applied on two types of neural network architectures, namely MLPs and radial basis function networks (RBFs), see for instance Andrieu, de Freitas and Doucet (1999);de Freitas (1999); Andrieu, de Freitas and Doucet (2000); de Freitas et al. (2001). For a review of Bayesian approaches to neural networks, see Titterington (2004).
Many research developments have been made to scale MCMC algorithms to big data. The main focus has been on designing Metropolis-Hastings or Gibbs sampling variants that evaluate a costly loglikelihood on a subset (minibatch) of the data rather than on the entire data set (Welling and Teh, 2011;Chen, Fox and Guestrin, 2014;Ma, Foti and Fox, 2017;Mandt, Hoffman and Blei, 2017;De Sa, Chen and Wong, 2018;Nemeth and Sherlock, 2018;Robert et al., 2018;Seita et al., 2018;Quiroz et al., 2019).
Among minibatch MCMC algorithms to big data applications, there exists a subset of studies applying such algorithms to neural networks (Chen, Fox and Guestrin, 2014;Gu, Ghahramani and Turner, 2015;Gong, Li and Hernández-Lobato, 2019). Minibatch MCMC approaches to neural networks pave the way towards data-parallel deep learning. On the other hand, to the best of the authors' knowledge, there is no published research on MCMC methods that evaluate the log-likelihood on a subset of neural network parameters rather than on the whole set of parameters, and therefore no reported research on model-parallel deep learning via MCMC.
Minibatch MCMC has been studied analytically by Johndrow, Pillai and Smith (2020). Their theoretical findings point out that some minibatching schemes can yield inexact approximations and that minibatch MCMC can not greatly expedite the rate of convergence.

Model structure
A neural network with ρ layers can be viewed as a hierarchical model with ρ levels, each network layer representing a level (Williams, 2000). Due to its nested layers and its non-linear activations, a neural network is a non-linear hierarchical model.
MCMC methods for non-linear hierarchical models have been developed, see for example Bennett, Racine-Poon and Wakefield; Gilks and Roberts; Daniels and Kass (1998) ;Sargent, Hodges and Carlin (2000). However, existing MCMC methods for non-linear hierarchical models have not harnessed neural networks due to time complexity and convergence implications.
Although not designed to mirror the hierarchical structure of a neural network, recent hierarchical VI (Ranganath, Tran and Blei, 2016;Esmaeili et al., 2019;Huang et al., 2019;Titsias and Ruiz, 2019) provides more general variational approximations of the parameter posterior of the neural network than mean-field VI. Introducing a hierarchi-cal structure in the variational distribution induces correlation among parameters, in contrast to the mean-field variational distribution that assumes independent parameters. So, one of the Bayesian inference strategies for neural networks is to approximate the covariance structure among network parameters. In fact, there are published comparisons between MCMC and VI in terms of speed and accuracy of convergence to the posterior covariance, both for linear or mixture models (Giordano, Broderick and Jordan, 2015;Mandt, Hoffman and Blei, 2017;Ong, Nott and Smith, 2018) and for neural networks (Zhang et al., 2018a).

Weight symmetries
The output of a feedforward neural network given some fixed input remains unchanged under a set of transformations determined by the the choice of activations and by the network architecture more generally. For instance, certain weight permutations and sign flips in MLPs with hyperbolic tangent activations leave the output unchanged (Chen, Lu and Hecht-Nielsen, 1993).
If a parameter transformation leaves the output of a neural network unchanged given some fixed input, then the likelihood is invariant under the transformation. In other words, transformations, such as weight permutations and sign-flips, render neural networks non-identifiable (Pourzanjani, Jiang and Petzold, 2017).
It is known that the set of linear invertible parameter transformations that leaves the output unchanged is a subgroup T of the group of invertible linear mappings from the parameter space R n to itself (Hecht-Nielsen, 1990). T is a transformation group acting on the parameter space R n . It can be shown that for each permutable feedforward neural network, there exists a cone H ⊂ R n dependent only on the network architecture such that for any parameter θ ∈ R n there exist η ∈ H and τ ∈ T such that τ η = θ. This relation means that every network parameter is equivalent to a parameter in the proper subset H of R n (Hecht-Nielsen, 1990). Neural networks with convolutions, max-pooling and batchnormalization contain more types of weight symmetries than MLPs (Badrinarayanan, Mishra and Cipolla, 2015).
In practice, the parameter space of a neural net-work is set to be the whole of R n rather than a cone H of R n . Since a neural network likelihood with support in the non-reduced parameter space of R n is invariant under weight permutations, sign-flips or other transformations, the posterior landscape includes multiple equally likely modes. This implies low acceptance rate, entrapment in local modes and convergence challenges for MCMC. Additionally, computational time is wasted during MCMC, since posterior modes represent equivalent solutions (Nalisnick, 2018). Such challenges manifest themselves in the MLP examples of section 4. For neural networks with higher number n of parameters in R n , the topology of the likelihood is characterized by local optima embedded in high-dimensional flat plateaus (Brea et al., 2019). Thereby, larger neural networks lead to a multimodal target density with symmetric modes for MCMC. Seeking parameter symmetries in neural networks can lead to a variety of NP-hard problems (Ensign et al., 2017). Moreover, symmetries in neural networks pose identifiability and associated inferential challenges in Bayesian inference, but they also provide opportunities to develop inferential methods with reduced computational cost (Hu, Zagoruyko and Komodakis, 2019) or with improved predictive performance (Moore, 2016). Empirical evidence from stochastic optimization simulations suggests that removing weight symmetries has a negative effect on prediction accuracy in smaller and shallower convolutional neural networks (CNNs), but has no effect in prediction accuracy in larger and deeper CNNs (Maddison et al., 2015).
Imposing constraints on neural network weights is one way of removing symmetries, leading to better mixing for MCMC (Sen, Papamarkou and Dunson, 2020). More generally, exploitation of weight symmetries provides scope for scalable Bayesian inference in deep learning by reducing the measure or dimension of parameter space. Bayesian inference in subspaces of parameter space for deep learning has been proposed before .
Lack of identifiability is not unique to neural networks. For instance, the likelihood of mixture models is invariant under relabelling of the mixture components, a condition known as the label switching problem (Stephens, 2000).
The high-dimensional parameter space of neu-ral networks is another source of non-identifiability. A necessary condition for identifiability is that the number of data points must be larger than the number of parameters. This is one reason why big datasets are required for training neural networks.

Prior specification
Parameter priors have been used for generating Bayesian smoothing or regularization effects. For instance, de Freitas (1999) develops sequential Monte Carlo methods with smoothing priors for MLPs and Williams (1995) introduces Bayesian regularization and pruning for neural networks via a Laplace prior.
When parameter prior specification for a neural network is not driven by smoothing or regularization, the question becomes how to choose the prior. The choice of parameter prior for a neural network is crucial in that it affects the parameter posterior (Lee, 2004), and consequently the posterior predictive distribution (Lee, 2005).
Neural networks are commonly applied to big data. For large amounts of data, practitioners may not have intuition about the relationship between input and output variables. Furthermore, it is an open research question how to interpret neural network weights and biases. As a priori knowledge about big datasets and about neural network parameters is typically not available, prior elicitation from experts is not applicable to neural networks.
It seems logical to choose a prior that reflects a priori ignorance about the parameters. A constantvalued prior is a possible candidate, with the caveat of being improper for unbounded parameter spaces, such as R n . However, for neural networks, an improper prior can result in an improper parameter posterior (Lee, 2005).
Typically, a truncated flat prior for neural networks is sufficient for ensuring a valid parameter posterior (Lee, 2005). At the same time, the choice of truncation bounds depends on weight symmetry and consequently on the allocation of equivalent points in the parameter space. Lee (2003) proposes a restricted flat prior for feedforward neural networks by bounding some of the parameters and by imposing constraints that guarantee layer-wise linear independence between activations, while Lee (2000) shows that this prior is asymptotically consistent for the posterior. Moreover, Lee (2003) demonstrates that such a restricted flat prior enables more effective MCMC sampling in comparison to alternative prior choices.
Objective prior specification is an area of statistics that has not infiltrated Bayesian inference for neural networks. Alternative ideas for constructing objective priors with minimal effect on posterior inference exist in the statistics literature. For example, Jeffreys priors are invariant to differentiable one-to-one transformations of the parameters (Jeffreys, 1962), maximum entropy priors maximize the Shannon entropy and therefore provide the least possible information (Jaynes, 1968), reference priors maximize the expected Kullback-Leibler divergence from the associated posteriors and in that sense are the least informative priors (Bernardo, 1979), and penalised complexity priors penalise the complexity induced by deviating from a simpler base model (Simpson et al., 2017).
To the best of the authors' knowledge, there are only two published lines of research on objective priors for neural networks; a theoretical derivation of Jeffreys and reference priors for feedforward neural networks by Lee (2007), and an approximation of reference priors via Monte Carlo sampling of a differentiable non-centered parameterization of MLPs and CNNs by Nalisnick (2018).
More broadly, research on prior specification for BNNs has been published recently (Pearce et al., 2019;Vladimirova et al., 2019). For a more thorough review of prior specification for BNNs, see Lee (2005).

Convergence
MCMC convergence depends on the target density, namely on its multi-modality and level of smoothness. An MLP with fewer than a hundred parameters fitted to a non-linearly separable dataset makes convergence in fixed MCMC sampling time challenging (see subsection 4.3).
Attaining MCMC convergence is not the only challenge. Assessing whether a finite sample from an MCMC algorithm represents an underlying target density can not be done with certainty (Cowles and Carlin, 1996). MCMC diagnostics can fail to detect the type of convergence failure they were designed to identify. Combinations of diagnostics are thus used in practice to evaluate MCMC convergence with re-duced risk of false diagnosis.
MCMC diagnostics were initially designed for asymptotically exact MCMC. Research activity on approximate MCMC has emerged recently. Minibatch MCMC methods (see subsection 2.1) are one class of approximate MCMC methods. Alternative approximate MCMC techniques without minibatching have been developed (Rudolf and Schweizer, 2018;Chen et al., 2019) along with new approaches to quantify convergence (Chwialkowski, Strathmann and Gretton, 2016).
Quantization and discrepancy are two notions pertinent to approximate MCMC methods. The quantization of a target density p by an empirical measurep provides an approximation to the target p (Graf and Luschgy, 2007), while the notion of discrepancy quantifies how well the empirical measurep approximates the target p (Chen et al., 2019). The kernel Stein discrepancy (KSD) and the maximum mean discrepancy (MMD) constitute two instances of discrepancy; for more details, see Chen et al. (2019) and Gretton et al. (2012), respectively. Rudolf and Schweizer (2018) provide an alternative way of assessing the quality of approximation of a target density p by an empirical measurep in the context of approximate MCMC using the notion of Wasserstein distance between p andp.

INFERENTIAL FRAMEWORK OVERVIEW
An overview of the inferential framework used in this paper follows, including the MLP model and its likelihood for classification, MCMC samplers for parameter estimation, MCMC diagnostics for assessing convergence and sampling effectiveness, and Bayesian marginalization for prediction.

Model definition.
An MLP is a feedforward neural network consisting of an input layer, one or more hidden layers and an output layer (Rosenblatt, 1958;Minsky and Papert, 1988;Hastie, Tibshirani and Friedman, 2016). Let ρ ≥ 2 be a natural number. Consider an index j ∈ {0, 1, . . . , ρ} indicating the layer, where j = 0 refers to the input layer, j = 1, 2, . . . , ρ − 1 to one of the ρ − 1 hidden layers and j = ρ to the output layer. Let κ j be the number of neurons in layer j and use κ 0:ρ = (κ 0 , κ 1 , . . . , κ ρ ) as a shorthand for the sequence of neuron counts per layer. Under such notation, MLP(κ 0:ρ ) refers to an MLP with ρ − 1 hidden layers and κ j neurons at layer j.
An MLP(κ 0:ρ ) with ρ − 1 ≥ 1 hidden layers and κ j neurons at layer j is defined recursively as for j = 1, 2, . . . , ρ. A data point x i ∈ R κ 0 corresponds to the input layer h 0 (x i ) = x i , yielding the sequence g 1 (x i , θ 1 ) = W 1 x i + b 1 in the first hidden layer. W j and b j are the respective weights and biases at layer j = 1, 2, . . . , ρ, which constitute the parameters θ j = (W j , b j ) at layer j. The shorthand θ 1:j = (θ 1 , θ 2 , . . . , θ j ) denotes all weights and biases up to layer j. Functions φ j , known as activations, are applied elementwise to their input g j .
The default recommendation of activation in neural networks is a rectified linear unit (ReLU), see for instance Jarrett et al. (2009);Nair and Hinton (2009) ;Goodfellow, Bengio and Courville (2016). Other activations are the ELU, leaky RELU, tanh and sigmoid (Nwankpa et al., 2018). If an activation is not present at layer j, then the identity function φ j (g j ) = g j is used as φ j in (3.2).
The weight matrix W j in (3.1) has κ j rows and κ j−1 columns, while the vector b j of biases has length κ j . Concatenating all θ j across hidden and output layers gives a parameter vector θ = θ 1:ρ ∈ R n of length n = ρ j=1 κ j (κ j−1 + 1). To define θ uniquely, the convention to traverse weight matrix elements row-wise is made. Apparently, each of g j in (3.1) and h j in (3.2) has length κ j .
The notation W j,k,l is introduced to point to the (k, l)-the element of weight matrix W j at layer j. Analogously, b j,k points to the k-th coordinate of bias vector b j at layer j.
Firstly, the sigmoid activation function φ ρ (g ρ ) = 1/(1 + exp (−g ρ )) is applied at the output layer of the MLP. So, the event probabilities Pr(y i = 1|x i , θ) are set to Assuming that the labels are outcomes of s independent draws from Bernoulli probability mass functions with event probabilities given by (3.3), the likelihood becomes ) of event probabilities for sample i = 1, 2, . . . , s. Furthermore, 1 denotes the indicator function, that is The negative value of log-likelihood (3.5) is known as the binary cross entropy (BCE). To infer the parameters θ of MLP(κ 0 , κ 1 , . . . , κ ρ = 1), the binary cross entropy or a different loss function is minimized using stochastic optimization methods, such as stochastic gradient descent (SGD).
Initially, a softmax activation function φ ρ (g ρ ) = exp (g ρ )/ κρ k=1 exp (g ρ,k ) is applied at the output layer of the MLP, where g ρ,j denotes the k-th coordinate of the κ ρ -length vector g ρ . Thus, the event probabilities Pr(y i = k|x i , θ) are . (3.6) It is assumed that the labels are outcomes of s independent draws from categorical probability mass functions with event probabilities given by (3.6), so the likelihood is The log-likelihood follows as The negative value of log-likelihood (3.8) is known as cross entropy, and it is used as loss function for stochastic optimization in multiclass classification MLPs.
An MLP(κ 0 , κ 1 , . . . , κ ρ = 2) with two neurons at the output layer, event probabilities given by softmax activation (3.6) and log-likelihood (3.8) can be used for binary classification. Such a formulation is an alternative to an MLP(κ 0 , κ 1 , . . . , κ ρ = 1) with one neuron at the output layer, event probabilities given by sigmoid activation (3.3) and log-likelihood (3.5). The difference between the two MLP models is the parameterization of event probabilities, since a categorical distribution with κ ρ = 2 levels otherwise coincides with a Bernoulli distribution.
The parameter posterior p(θ|x 1:s , y 1:s ) is alternatively denoted by p(θ|D 1:s ) for brevity. D 1:s = (x 1:s , y 1:s ) is a dataset of size s consisting of input x 1:s and output y 1:s .
This subsection provides an introduction to the MCMC algorithms and MCMC diagnostics used in the examples of section 4. Three MCMC algorithms are outlined, namely Metropolis-Hastings, Hamiltonian Monte Carlo, and power posterior sampling. Two MCMC diagnostics are described, the multivariate potential scale reduction factor (PSRF) and the multivariate effective sample size (ESS).

Metropolis-Hastings algorithm.
One of the most general algorithms for sampling from a posterior p(θ|D 1:s ) is the Metropolis-Hastings (MH) algorithm (Metropolis et al., 1953;Hastings, 1970). Given the current state θ, the MH algorithm initially samples a state θ * from a proposal density g θ and subsequently accepts the proposed state θ * with probability Typically, a normal proposal density g θ = N (θ, Λ) with a constant covariance matrix Λ is used. For such a normal g θ , the acceptance probability simplifies to min {p(θ * |D 1:s )/p(θ|D 1:s ), 1}, yielding the so called random walk Metropolis algorithm.

Hamiltonian Monte Carlo. Hamiltonian
Monte Carlo (HMC) draws samples from an augmented parameter space via Gibbs steps, by computing a trajectory in the parameter space according to Hamiltonian dynamics. For a more detailed review of HMC, see Neal (2011).

Power posterior sampling.
Power posterior (PP) sampling by Friel and Pettitt (2008) is a population Monte Carlo algorithm. It involves m + 1 chains drawn from tempered versions p t i (θ|D 1:s ) of a target posterior p(θ|D 1:s ) for a temperature schedule t i ∈ [0, 1], i ∈ {0, 1, . . . , m}, where t m = 1. At each iteration, the state of each chain is updated using an MCMC sampler associated with that chain and subsequently states between pairs of chains are swapped according to an MH algorithm. For the i-th chain, a sample j is drawn from a probability mass function p i with probability p i (j), in order to determine the pair (i, j) for a possible swap.
Power posteriors p t i (θ|D 1:s ), t i < t m , are smooth approximations of the target density p tm (θ|D 1:s ) = p(θ|D 1:s ), facilitating exploration of the parameter space via state transitions between chains of p t i (θ|D 1:s ) and of p(θ|D 1:s ). In this paper, a categorical probability mass function p i is used in PP sampling for determining candidate pairs of chains for state swaps (see Appendix A).
3.2.4 Multivariate PSRF. PSRF, commonly denoted byR, is an MCMC diagnostic of convergence conceived by Gelman and Rubin (1992) and extended to its multivariate version by Brooks and Gelman (1998). This paper uses the multivariate PSRF by Brooks and Gelman (1998), which provides a single-number summary of convergence across the n dimensions of a parameter, requiring a Monte Carlo covariance matrix estimator for the parameter.
To acquire the multivariate PSRF, the multivariate initial monotone sequence estimator (MINSE) of Monte Carlo covariance is employed (Dai and Jones, 2017). In a Bayesian setting, the MINSE estimates the covariance matrix of a parameter posterior p(θ|D 1:s ).
To compute PSRF, several independent Markov chains are simulated. Gelman et al. (2004) recommend terminating MCMC sampling as soon asR < 1.1. More recently, Vats and Knudson (2018) make an argument based on ESS that a cut-off of 1.1 for R is too high to estimate a Monte Carlo mean with reasonable uncertainty. Vehtari et al. (2019) recommend simulating at least m = 4 chains to computê R and using a threshold ofR < 1.01.

Multivariate ESS.
The ESS of an estimate obtained from a Markov chain realization is interpreted as the number of independent samples that provide an estimate with variance equal to the variance of the estimate obtained from the Markov chain realization. For a more extensive treatment entailing univariate approaches to ESS, see Vats and Flegal (2018); Gong and Flegal (2016); Kass et al. (1998).
R and its variants can fail to diagnose poor mixing of a Markov chain, whereas low values of ESS are an indicator of poor mixing. It is thus recommended to check bothR and ESS (Vehtari et al., 2019). For a theoretical treatment of the relation betweenR and ESS, see Vats and Knudson (2018).
Univariate ESS pertains to a single coordinate of an n-dimensional parameter. Vats, Flegal and Jones (2019) introduce a multivariate version of ESS, which provides a single-number summary of sampling effectiveness across the n dimensions of a parameter. Similarly to multivariate PSRF (Brooks and Gelman, 1998)

Bayesian marginalization for prediction
This subsection briefly reviews the notion of posterior predictive distribution based on Bayesian marginalization, posterior predictive distribution approximation via Monte Carlo integration, and associated binary and multiclass classification. Appendix B provides a derivation of (3.9). (3.10) states the posterior predictive distribution p(y|x, D 1:s ) as an expectation of the likelihood p(y|x, θ) evaluated at the test output y with respect to the parameter posterior p(θ|D 1:s ) learnt from the training set D 1:s . The expectation in (3.10) can be approximated via Monte Carlo integration. More specifically, a Monte Carlo approximation of the posterior predictive distribution is given by The sum in (3.11) involves evaluations of the likelihood across v iterations ω k , k = 1, 2, . . . , v, of a Markov chain realization ω 1:v obtained from the parameter posterior p(θ|D 1:s ).

Classification rule.
In the case of binary classification, the predictionŷ for the test label y ∈ {0, 1} is The classification rules (3.12) and (3.13) for binary and multiclass classification maximize the posterior predictive distribution. This way, predictions are made based on the Bayesian principle. The uncertainty of predictions is quantified, since the posterior predictive probability p(y|x, D 1:s ) of each predicted labelŷ is available.

EXAMPLES
Four examples of Bayesian inference for MLPs based on MCMC are presented. A different dataset is used for each example. The four datasets entail simulated noisy data from the exclusive-or (XOR) function, and observations collected from Pima Indians, penguins and hawks. Section 4.1 introduces the four datasets. Each of the four datasets is split into a training and a test set for parameter inference and for predictions, respectively. MLPs with one neuron in the output layer are fitted to the noisy XOR and Pima datasets to perform binary classification, while MLPs with three neurons in the output layer are fitted to the penguin and hawk datasets to perform multiclass classification with three classes. Table 1 shows the training and test sample sizes of the four datasets, and the fitted MLP models with their associated number n of parameters.

Datasets
An introduction to the four datasets used in this paper follows. The simulated noisy XOR dataset does not contain missing values, while the real datasets for Pima, penguins and hawks come with missing values. Data points containing missing values in the chosen variables have been dropped from the three real datasets. All features (input variables) in the three real datasets have been standardized. The four datasets, in their final form used for inference and prediction, are available at https: //github.com/papamarkou/bnn_mcmc_examples. A perceptron without a hidden layer can not learn the XOR function (Minsky and Papert, 1988). On the other hand, an MLP(2, 2, 1) with a single hidden layer of two neurons can learn the XOR function (Goodfellow, Bengio and Courville, 2016).
An MLP(2, 2, 1) has a parameter vector θ of length n = 9, as W 1 , b 1 , W 2 and b 2 have respective dimensions 2 · 2, 2 · 1, 2 · 1 and 1 · 1. Since the number s = 4 of data points defined by the exact XOR function is less than the number n = 9 of parameters in the fitted MLP(2, 2, 1), the parameters can not be fully identified.
To circumvent the lack of identifiability arising from the limited number of data points, a larger dataset is simulated by introducing a noisy version of XOR. Firstly, consider the auxiliary function ψ : ψ is presented in parametrized form, in terms of a constant c ∈ (0.5, 1) and a uniformly distributed random variable u ∼ U(0, 1). The noisy XOR function is then defined as the function composition f •ψ.
A training and a test set of noisy XOR points, generated using f • ψ and c = 0.55, are shown in figure 2a. 125 and 30 noisy XOR points per exact XOR point (x i , y i ), i = 1, 2, 3, 4, are contained in the training and test set, respectively. So, the training and test sample sizes are 500 and 120, as reported in table 1 and as visualized in figure 2a.

Pima dataset.
The Pima dataset contains observations taken from female patients of Pima Indian heritage. The binary output variable indicates whether or not a patient has diabetes. Eight features are used as diagnostics of diabetes, namely the number of pregnancies, plasma glucose concentration, diastolic blood pressure, triceps skinfold thickness, insulin level, body mass index, diabetes pedigree function and age.
For more information about the Pima dataset, see Smith et al. (1988). The original data, prior to removal of missing values and feature standardization, are available as the PimaIndiansDiabetes2 data frame of the mlbench R package.

Penguin dataset.
The penguin dataset consists of body measurements for three penguin species observed on three islands in the Palmer Archipelago, Antarctica. Adélie, Chinstrap and Gentoo penguins are the three observed species. Four body measurements per penguin are taken, specifically body mass, flipper length, bill length and bill depth. The four body measurements, sex and location (island) make up a total of six features utilized for deducing the species to which a penguin belongs. Thus, the penguin species is used as output variable.
Horst, Hill and Gorman (2020) provide more details about the penguin dataset. In their original form, prior to data filtering, the data are available at https://github.com/allisonhorst/ palmerpenguins.

Hawk dataset.
The hawk dataset is composed of observations for three hawk species collected from Lake MacBride near Iowa City, Iowa. Cooper's, red-tailed and sharp-shinned hawks are the three observed species. Age, wing length, body weight, culmen length, hallux length and tail length are the six hawk features employed in this paper for deducing the species to which a hawk belongs. So, the hawk species is used as output variable. Cannon et al. (2019) mention that Emeritus Professor Bob Black at Cornell College shared the hawk dataset publicly. The original data, prior to data filtering, are available as the Hawks data frame of the Stat2Data R package.

Experimental configuration
To fully specify the MLP models of table 1, their activations are listed. A sigmoid activation function is applied at each hidden layer of each MLP. Additionally, a sigmoid activation function is applied at the output layer of MLP(2, 2, 1) and of MLP(8, 2, 2, 1), conforming to log-likelihood (3.5) for binary classification. A softmax activation function is applied at the output layer of MLP(6, 2, 2, 3), in accordance with log-likelihood (3.8) for multiclass classification. The same MLP(6, 2, 2, 3) model is fitted to the penguin and hawk datasets.
A normal prior π(θ) = N (0, 10I) is adopted for the parameters θ ∈ R n of each MLP model shown in table 1. An isotropic covariance matrix 10I assigns relatively high prior variance, equal to 10, to each coordinate of θ, thus setting empirically a seemingly non-informative prior.
MH and HMC are run for each of the four examples of table 1. PP sampling incurs higher computational cost than MH and HMC; for this reason, PP sampling is run only for noisy XOR. Ten power posteriors are employed for PP sampling, and MH is used for within-chain moves. On the basis of pilot runs, the PP temperature schedule is set to t i = 1, i = 0, 1, . . . , 9; this implies that each power posterior is set to be the parameter posterior and consequently between-chain moves are made among ten chains realized from the parameter posterior. Empirical hyperparameter tuning for MH, HMC and PP is carried out. The chosen MH proposal variance, HMC number of leapfrog steps and HMC leapfrog step size for each example can be found in https: //github.com/papamarkou/bnn_mcmc_examples. m = 10 Markov chains are realized for each combination of training dataset shown in table 1 and of MCMC sampler. 110, 000 iterations are run per chain realization, 10, 000 of which are discarded as burn-in. Thereby, v = 100, 000 post-burnin iterations are retained per chain realization. MINSE computation, required by multivariate PSRF and multivariate ESS, is carried out using v = 100, 000 post-burnin iterations per realized chain. The multivariate PSRF for each dataset-sampler setting is computed across the m = 10 realized chains for the setting. On the other hand, the multivariate ESS is computed for each realized chain, and the mean across m = 10 ESSs is reported for each dataset-sampler setting.
Monte Carlo approximations of posterior predictive distributions are computed according to (3.11) for each data point of each test set. To reduce the computational cost, the last v = 10, 000 iterations of each realized chain are used in (3.11).
Predictions for binary and multiclass classifica-tion are made using (3.12) and (3.13), respectively. Given a single chain realization from an MCMC sampler, predictions are made for every point in a test set; the predictive accuracy is then computed as the number of correct predictions over the total number of points in the test set. Subsequently, the mean of predictive accuracies across the m = 10 chains realized from the sampler is reported for the test set. Table 2 shows numerical summaries for each set of m = 10 Markov chains realized by an MCMC sampler for a dataset-MLP combination of table 1. Multivariate PSRF and multivariate ESS diagnose the capacity of MCMC sampling to perform parameter inference. Predictive accuracy via Bayesian marginalization (3.11), based on classification rules (3.12) and (3.13) for binary and multiclass classification, demonstrates the predictive performance of MCMC sampling. The last column of table 2 displays the predictive accuracy via (3.11) with samples ω k , k = 1, 2, . . . , v, drawn from the prior π(θ) = N (0, 10I), thus providing an approximation of the expected posterior predictive probability (4.1) E θ [p(y|x, θ)] = p(y|x, θ)π(θ)dθ with respect to prior π(θ). PSRF is above 1.01 (Vehtari et al., 2019), indicating lack of convergence, in three out of four datasets. ESS is low considering the post-burnin length of v = 100, 000 of each chain realization, indicating slow mixing. MCMC sampling for Pima data is the only case of attaining PSRF less than 1.01, yet the ESS values for Pima are the lowest among the four datasets. Overall, simultaneous low PSRF and high ESS are not reached in any of the examples.

Numerical summaries
The predictive accuracy is high in multiclass classification, despite the lack of convergence and slow mixing. Bayesian marginalization based on HMC samples yields 100% and 98.03% predictive accuracy on the penguin and hawk test datasets, despite the PSRF values of 1.6082 and 1.4421 on the penguin and hawk training datasets.
PP sampling for the binary classification problem of noisy XOR leads to higher predictive accuracy (87.58%) than MH (75.92%) or HMC (74.75%). The 87.58% predictive accuracy is attained by PP sampling despite the associated PSRF value of 2.2885. Bayesian marginalization based on MCMC sampling outperforms prior beliefs or random guesses in terms of predictive inference, despite MCMC diagnostic failures. For instance, Bayesian marginalization via non-converged HMC chain realizations yields 74.75%, 100% and 98.03% predictive accuracy on the noisy XOR, penguin and hawk datasets. Approximating the posterior predictive distribution with samples from the parameter prior yields 48.33%, 36.45% and 28.85% predictive accuracy on the same datasets. It is noted that 48.33% is close to a 50/50 random guess for binary classification, while 36.45% and 28.85% are close to a 1/3 random guess for multiclass classification with three classes.

Visual summaries for parameters
Visual summaries for MLP parameters are presented in this subsection. In particular, Markov chain traceplots and a comparison between MCMC sampling and ensemble training are displayed. models introduced in table 1. These traceplots visually demonstrate entrapment in local modes, mode switching and more generally lack of convergence.
All 110, 000 iterations per realized chain, which include burnin, are shown in the traceplots of figure 1. The vertical dotted lines delineate the first 10, 000 burnin iterations.
Two realized MH chains for parameter θ 8 of the MLP(2, 2, 1) model fitted to the noisy XOR training data are plotted. The traces in orange and in blue gravitate during burnin towards modes in the vicinity of 8 and −8, respectively, and then get entrapped for the entire simulation time in these modes. Parameter θ 8 corresponds to a weight connecting a neuron in the hidden layer with the neuron of the output layer of MLP(2, 2, 1). The two realized chains for θ 8 explore two regions symmetric about zero associated with symmetries of weight θ 8 .
Two realized MH chains for parameter θ 18 of the MLP(6, 2, 2, 3) model fitted to the penguin training data are plotted, one shown in orange and one in blue. Each of these two traces initially explore a mode, transit to a seemingly symmetric mode about halfway through the simulation time (post-burnin) and explore the symmetric mode in the second half of the simulation.
One HMC chain traceplot for parameter θ 23 and one HMC chain traceplot for parameter θ 26 of the MLP(6, 2, 2, 3) model fitted to the penguin and hawk training data, respectively, are shown. The traces of these two parameters exhibit similar behaviour, each of them switching between two symmetric regions about zero.
Switching between symmetric modes, as seen in the displayed traceplots, manifests weight symmetries. These traceplots exemplify how computational time is wasted during MCMC to explore equivariant parameter posterior modes of a neural network (Nalisnick, 2018). Consequently, the realized chains do not converge.

MCMC sampling vs ensemble training.
An exemplified comparison between MCMC sampling and ensemble training for neural networks follows. To this end, the same noisy XOR training data and the same MLP(2, 2, 1) model, previously used for MCMC sampling, are used for ensemble training.
To recap, the noisy XOR dataset is introduced in subsection 4.1 and is displayed in figure 2a; a sigmoid activation function is applied to the hidden and output layer of MLP(2, 2, 1), and the BCE loss function is employed, which is the negative value of log-likelihood (3.5).
Ensemble learning is conducted by training the MLP(2, 2, 1) model on the noisy XOR training set multiple times. At each training session, SGD is used for minimizing the BCE loss. SGD is initialized by drawing a sample from π(θ) = N (0, 10I), which is the same density used as prior for MCMC sampling. 2, 000 epochs are run per training session, with a batch size of 50 and a learning rate of 0.002. The SGD solution from the training session is accepted if its predictive accuracy on the noisy XOR test set is above 85%, otherwise it is rejected. Ensemble learning is terminated as soon as 1, 000 SGD solutions with the required level of accuracy are obtained.
(a) Noisy XOR training set (left) and test set (right).

Fig 2:
Comparison between MH sampling and ensemble training of an MLP(2, 2, 1) model fitted to noisy XOR data. SGD is used for ensemble training. Each accepted SGD solution has predictive accuracy above 85% on the noisy XOR test set.   2, 2, 1). The green histogram represents all 1, 000 SGD solutions for θ 3 obtained from ensemble training based on noisy XOR. These 1, 000 modes cluster in two regions approximately symmetric about zero. The orange histogram belongs to one of ten realized MH chains for θ 3 based on noisy XOR. This realized chain is entrapped in a local mode in the vicinity of 5, where the orange histogram con-centrates its mass. The overlaid green and orange histograms show that MH sampling explores a region of the marginal posterior of θ 3 also explored by ensemble training.
The blue histogram in figure 2c comes from a chain realization for θ 3 using MH sampling to apply MLP(2, 2, 1) to the four exact XOR data points. The pink line in figure 2c shows the marginal prior π(θ 3 ) = N (0, σ 2 = 10). Four data points are not sufficient to learn from them, given that MLP(2, 2, 1) has nine parameters. For this reason, the blue histogram coincides with the pink line, which means that the marginal posterior p(θ 3 ) obtained from exact XOR via MH sampling and the marginal prior π(θ 3 ) coincide.

Visual summaries for predictions
Visual summaries for MLP predictions and for MLP posterior predictive probabilities are presented in this section. MLP posterior predictive probabilities are visually shown to quantify predictive uncertainty in classification. Figure 3 shows boxplots of predictive accuracies, hereinafter referred to as accuracies, for the examples introduced in table 1. Each boxplot summarizes m = 10 accuracies associated with the ten chains realized per sampler for a test set. Accuracy computation is based on Bayesian marginalization, as outlined in subsections 3.3 and 4.2. Horizontal red lines represent accuracy medians. Figure 3 and table 1 provide complementary summaries, as they present respective quartiles and means of accuracies across chains per sampler.

Predictive accuracy.
Boxplot medians show high accuracy on the penguin and hawk test sets. Moreover, narrow boxplots indicate accuracies with small variation on the penguin and hawk test sets. Thereby, Bayesian marginalization based on non-converged chain realizations attains high accuracy with small variability on the two multiclass classification examples. Figure 3 also displays boxplots of accuracies based on expected posterior predictive distribution approximation (4.1) with respect to the prior. For all four test sets and regardless of Markov chain convergence, Bayesian marginalization outperforms agnostic prior-based baseline (4.1).
The PP boxplot has more elevated median and is narrower than its MH and HMC counterparts for the noisy XOR test set. This implies that PP sampling attains higher accuracy with smaller variation than MH and HMC sampling on the noisy XOR test set.  5, 1.5]. D 1:500 refers to the noisy XOR training dataset of size s = 500 introduced in subsection 4.1. (3.11) is used for approximating p(y = 1|(x 1 , x 2 ), D 1:500 ). Previously acquired Markov chain realizations (subsection 4.3) via MCMC sampling of MLP(2, 2, 1) parameters, using the noisy XOR training dataset D 1:500 , are passed to (3.11).
The approximation p(y = 1|(x 1 , x 2 ), D 1:500 ) = c at the center (x 1 , x 2 ) of a square cell determines the colour of the cell in figure 4. If c is closer to 1, 0, or 0.5, the cell is plotted with a shade of red, blue or white, respectively. So, darker shades of red indicate that y = 1 with higher certainty, darker shades of blue indicate that y = 0 with higher certainty, and shades of white indicate high uncertainty about the binary label of noisy XOR.
Two posterior predictive distribution approximations based on two HMC chain realizations learn different regions of the exact posterior predictive distribution. Each of the two HMC chain realizations uncover about half of the ground truth of grid la- bels, while it remains highly uncertain for the other half of grid labels. Moreover, both HMC chain realizations exhibit higher uncertainty closer to the decision boundaries of ground truth. These decision boundaries are the vertical straight line x 1 = 0.5 and horizontal straight line x 2 = 0.5.
A posterior predictive distribution approximation based on a PP chain realization is displayed. PP sampling uncovers larger regions of the ground truth of grid labels than HMC sampling in the considered grid of noisy XOR features (x 1 , x 2 ). Although HMC and PP samples do not converge to the parameter posterior of MLP(2, 2, 1), approximate Bayesian marginalization using these samples predicts a subset of noisy XOR labels. Figures 5 and 6 show approximations of predictive prosterior probabilities for a binary classification (noisy XOR) and a multiclass classification (hawks) example. Two posterior predictive probabilities are (a) Scatterplot of noisy XOR features (x 1 , x 2 ).

Uncertainty quantification on a test set.
(b) Posterior predictive probabilities for noisy XOR. interpreted contextually in each example to quantify predictive uncertainty. Figure 5a visualizes the noisy XOR test set of subsection 4.1. This is the same test set shown in figure 2a, but with test points coloured according to their labels. Figure 5b shows the posterior predictive probability p(y = c|(x 1 , x 2 ), D 1:500 ) of true label c ∈ {0, 1} for each noisy XOR test point ((x 1 , x 2 ), y = c) given noisy XOR training set D 1:500 of subsection 4.1. The posterior probabilities p(y = c|(x 1 , x 2 ), D 1:500 ) of predicting true class c are ordered within class c. Moreover, each p(y = c|(x 1 , x 2 ), D 1:500 ) is coloured as red or pale green depending on whether the resulting prediction is correct or not. One of the ten MH chain realizations for MLP(2, 2, 1) parameter inference from noisy XOR data is used for approximating p(y = c|(x 1 , x 2 ), D 1:500 ) via (3.11) and for making predictions via (3.12).
Two points in the noisy XOR test set are marked in figure 5 using a square and a rhombus. These two points have the same true label c = 1. Given posterior predictive probabilities 0.5269 and 0.9750 for the rhombus and square-shaped test points, the label c = 1 is correctly predicted for both points. However, the rhombus-shaped point is closer to the decision boundary x 2 = 0.5 than the square-shaped point, so classifying the former entails higher uncertainty. As 0.5269 < 0.9750, Bayesian marginalization quantifies the increased predictive uncertainty associated with the rhombus-shaped point despite using a nonconverged MH chain realization. Figure 6a shows a scatterplot of weight against tail length for the hawk test set of subsection 4.1. Blue, red and green test points belong to Cooper's, red-tailed and sharp-shinned hawk classes. Figure 6b shows the posterior predictive probabilities p(y = c|x, D 1:596 ) for a subset of 100 hawk test points, where c ∈ {Cooper's, red-tailed, sharp-shinned} denotes the true label of test point (x, y = c) and D 1:596 denotes the hawk training set of subsection 4.1. These posterior predictive probabilities are shown ordered within each class, and are coloured red or pale green depending on whether they yield correct or wrong predictions. One of the ten MH chain realizations for MLP(6, 2, 2, 3) parameter inference is used for approximating p(y = c|x, D 1:596 ) via (3.11) and for making predictions via (3.13).
Two points in the hawk test set are marked in figure 6 using a square and a rhombus. Each of these two points represents weight and tail length measurements from a red-tailed hawk. The red-tailed hawk class is correctly predicted for both points. The squared-shaped observation belongs to the main cluster of red-tailed hawks in figure 6a and it is predicted with high posterior predictive probability (0.9961). On the other hand, the rhombus-shaped observation, which falls in the cluster of Cooper's hawk, is correctly predicted with a lower posterior predictive probability (0.5271). Bayesian marginalization provides approximate posterior predictive probabilities that signify the level of uncertainty in predictions despite using a non-converged MH chain realization.

Source code
The source code for this paper is split into three Python packages, namely eeyore, kanga and (a) Scatterplot of hawks' weight against tail length.
(b) Posterior predictive probabilities for hawks. eeyore is available via pip, via conda and at https://github.com/papamarkou/eeyore. eeyore implements the MLP model, as defined by (3.1)-(3.2), using PyTorch. An MLP class is set to be a subclass of torch.nn.Module, with log-likelihood (3.5) for binary classification equal to the negative value of torch.nn.BCELoss and with log-likelihood (3.8) for multiclass classification equal to the negative value of torch.nn.CrossEntropyLoss. Each MCMC algorithm takes an instance of torch.nn.Module as input, with the logarithm of the target density being a log target method of the instance. Logtarget density gradients for HMC are computed via the automatic differentiation functionality of the torch.autograd package of PyTorch. The MLP class of eeyore provides a predictive posterior method, which implements the posterior predictive distribution approximation (3.11) given a realized Markov chain.
kanga is available via pip, via conda and at https://github.com/papamarkou/kanga. kanga is a collection of MCMC diagnostics implemented using numpy. MINSE, multivariate PSRF multivariate ESS are available in kanga.
bnn mcmc examples organizes the examples of this paper in a package. bnn mcmc examples relies on eeyore for MCMC simulations and posterior predictive distribution approximations, and on kanga for MCMC diagnostics. For more details, see https: //github.com/papamarkou/bnn_mcmc_examples.
Optimization via SGD for the example involving MLP(2, 2, 1) and noisy XOR data (figure 2) is run using PyTorch. The loss function for optimization is computed via torch.nn.BCELoss. This loss function corresponds to the negative log-likelihood function (3.5) involved in MCMC, thus linking the SGD and MH simulations shown in figure 2c. SGD is coded manually instead of calling an optimization algorithm of the torch.optim package of PyTorch. Gradients for optimization are computed calling the backward method. The SGD code related to the example of figure 2 is available at https://github. com/papamarkou/bnn_mcmc_examples.

Hardware
Pilot MCMC runs indicated an increase in speed by using CPUs instead of GPUs; accordingly, computations were performed on CPUs for this paper. The GPU slowdown is explained by the overhead of copying PyTorch tensors between GPUs and CPUs for small neural networks, such as the ones used in section 4.
The computations for section 4 were run on Google Cloud Platform (GCP). Eleven virtual machine (VM) instances with virtual CPUs were created on GCP to spread the workload.
Setting aside heterogeneities in hardware configuration between GCP VM instances and in order to provide an indication of computational cost, MCMC simulation runtimes are provided for the example of applying an MLP(6, 2, 2, 3) to the hawk training dataset. The mean runtimes across the ten realized chains per MH and HMC are 0 : 42 : 54 and 1 : 10 : 48, respectively (runtimes are formatted as 'hours : minutes : seconds').

PREDICTIVE INFERENCE SCOPE
Bayesian marginalization can attain high predictive accuracy and can quantify predictive uncertainty using non-converged MCMC samples of neural network parameters. Thus, MCMC sampling elicits some information about the parameter posterior of a neural network and conveys such information to the posterior predictive distribution. It is possible that MCMC sampling learns about the statistical dependence among neural network parameters. Along these lines, groups of weights or biases can be formed, with strong within-group and weak between-group dependence, to investigate scalable block Gibbs sampling methods for neural networks.
Another possibility of MCMC developments for neural networks entails shifting attention from the parameter space to the output space, since the latter is related to predictive inference directly. Approximate MCMC methods that measure the discrepancy or Wasserstein distance between neural network predictions and output data (Rudolf and Schweizer, 2018) can be investigated.
Bayesian marginalization provides scope to develop predictive inference for neural networks. For instance, Bayesian marginalization can be examined in the context of approximate MCMC sampling from a neural network parameter posterior, regardless of convergence to the parameter posterior and in analogy to the workings of this paper. Moreover, the idea of Wilson and Izmailov (2020) to interpret ensemble training of neural networks from a viewpoint of Bayesian marginalization can be studied using the notion of quantization of probability distributions.

APPENDIX A: POWER POSTERIORS
This appendix provides the probability mass function p i (j) for proposing a chain j for a possible swap of states between chains i and j in PP sampling. Assuming m+1 power posteriors, a neighbouring chain j of i is chosen randomly from the categorical probability mass function p i = C(α i (0), α i (1), . . . , α i (i − 1), α i (i + 1), . . . , α i (m)) with event probabilities where i ∈ {0, 1, . . . , m}, j ∈ {0, 1, . . . , m} \ {i}, β is a hyperparameter and γ i is a normalizing constant. The hyperparameter β is typically set to β = 0.5, a value which makes a jump to j = i ± 1 roughly three times more likely than a jump to j = i ± 3 (Friel and Pettitt, 2008). The normalizing constant γ i is given by Starting from the fact that the event probabilities α i (j) add up to one, γ i is derived as follows: