Bayesian neural networks via MCMC: a Python-based tutorial
Bayesian inference provides a methodology for parameter estimation and
uncertainty quantification in machine learning and deep learning methods.
Variational inference and Markov Chain Monte-Carlo (MCMC) sampling methods are
used to implement Bayesian inference. In the past three decades, MCMC sampling
methods have faced some challenges in being adapted to larger models (such as
in deep learning) and big data problems. Advanced proposal distributions that
incorporate gradients, such as a Langevin proposal distribution, provide a
means to address some of the limitations of MCMC sampling for Bayesian neural
networks. Furthermore, MCMC methods have typically been constrained to
statisticians and currently not well-known among deep learning researchers. We
present a tutorial for MCMC methods that covers simple Bayesian linear and
logistic models, and Bayesian neural networks. The aim of this tutorial is to
bridge the gap between theory and implementation via coding, given a general
sparsity of libraries and tutorials to this end. This tutorial provides code in
Python with data and instructions that enable their use and extension. We
provide results for some benchmark problems showing the strengths and
weaknesses of implementing the respective Bayesian models via MCMC. We
highlight the challenges in sampling multi-modal posterior distributions for
the case of Bayesian neural networks and the need for further improvement of
convergence diagnosis methods.
Rohitash Chandra, Joshua Simmons
# Bayesian neural networks via MCMC: a Python-based tutorial

###### Abstract

Bayesian inference provides a methodology for parameter estimation and uncertainty quantification in machine learning and deep learning methods. Variational inference and Markov Chain Monte-Carlo (MCMC) sampling techniques are used to implement Bayesian inference. In the past three decades, MCMC methods have faced a number of challenges in being adapted to larger models (such as in deep learning) and big data problems. Advanced proposals that incorporate gradients, such as a Langevin proposal distribution, provide a means to address some of the limitations of MCMC sampling for Bayesian neural networks. Furthermore, MCMC methods have typically been constrained to use by statisticians and are still not prominent among deep learning researchers. We present a tutorial for MCMC methods that covers simple Bayesian linear and logistic models, and Bayesian neural networks. The aim of this tutorial is to bridge the gap between theory and implementation via coding, given a general sparsity of libraries and tutorials to this end. This tutorial provides code in Python with data and instructions that enable their use and extension. We provide results for some benchmark problems showing the strengths and weaknesses of implementing the respective Bayesian models via MCMC. We highlight the challenges in sampling multi-modal posterior distributions in particular for the case of Bayesian neural networks, and the need for further improvement of convergence diagnosis.

keywords: Bayesian neural networks, MCMC, Langevin dynamics, Deep learning, Convolutional Neural Networks, time series prediction +
Footnote †: journal: Elsevier

## 1 Introduction

Bayesian inference provides a probabilistic approach for parameter estimation in a wide range of models used across the fields of machine learning, econometrics, environmental and Earth sciences [1; 2; 3; 1; 4; 5]. The term 'probabilistic' refers to representation of unknown parameters as probability distributions rather than using fixed point estimates as in conventional optimisation methods and machine learning models where gradient-based methods are prominent [6]. The probability representation of unknown parameters requires a different approach to optimisation and, from the machine learning and computational statistics point-of-view, this is known as sampling.

Markov Chain Monte-Carlo (MCMC) sampling methods have been prominent for inference (estimation) of model parameters via the posterior distribution (probability distribution). In other words, Bayesian methods attempt to quantify the uncertainty in model parameters by marginalising over the predictive posterior distribution. Hence, in the case of neural networks, MCMC methods can be used to implement Bayesian neural networks that represent weights and biases as probability distributions [7; 8; 9; 10; 11; 7]. Probabilistic machine learning provides natural way of providing uncertainty quantification in predictions [12] since the uncertainties can be obtained by probabilistic representation of parameters. This inference procedure can be seen as a form of learning (optimisation) applied to the model parameters [8]. In this tutorial, we employ linear models and simple neural networks to demonstrate the use
of MCMC sampling methods. The probabilistic representation of weights and biases in the respective models allows uncertainty quantification on model predictions.

We note that MCMC refers to a family of algorithms for implementing Bayesian inference for parameter and uncertainty estimation in models. The range of models to which Bayesian inference is applied, including statistical, graphical, and machine learning models; has led to the existence of a wide range of MCMC sampling algorithms. Some of the prominent ones are Metropolis Hastings algorithm [13; 14; 15], Gibbs sampling [16; 17; 18], Langevin MCMC [19; 20; 21], rejection sampling [22; 23], sequential MCMC [24], adaptive MCMC [25], parallel tempering (tempered) MCMC [26; 27; 28; 29], reversible-jump MCMC [30; 31], specialised MCMC methods for discrete time series models[32; 33; 34], constrained parameter and model settings [35; 36], and likelihood free MCMC [37]. MCMC sampling methods have also been used for data augmentation [38; 39], model fusion [40], model selection [41; 42], and interpolation [43]. Apart from this, we note that MCMC methods have been prominent in a wide range of applications that include geophysical inversions [44; 45; 46], geo-scientific models [5; 47; 48], environmental and hydrological modelling [49; 50], bio-systems modelling [51; 52; 53], and quantitative genetics [54; 55].

In the case of Bayesian neural networks, the large number of model parameters that emerge from large neural network architectures and deep learning models pose challenges for MCMC sampling methods. Hence, progress in the application of Bayesian approaches to big data and deep neural networks has been slow. Research in this space has included a number of methods that have been fused with MCMC such as gradient-based methods [56; 21; 57; 58; 59], and evolutionary (meta-heuristic) algorithms which include differential evolution, genetic algorithms, and particle swarm optimisation [60; 61; 62; 63].

This use of gradients in MCMC was initially known as Metropolis-adjusted Langevin dynamics [21] and has shown promising performance for linear models [58] and has also been extended to Bayesian neural networks [59]. Another direction has been the use of better exploration features in MCMC sampling such as parallel tempering MCMC with Langevin-gradients and parallel computing [59], providing a competitive alternative to stochastic gradient-descent [64] and Adam optimizers [6] with the addition of uncertainty quantification in predictions. These methods have also been applied to Bayesian deep learning models such as Bayesian autoencoders [65] and Bayesian graph convolutional neural networks (CNNs) [66] which require millions of trainable parameters to be represent as posterior distributions. Recently, Kapoor et al. [63] combined tempered MCMC with particle swarm optimisation-based proposal distribution in an parallel computing environment and showed that the hybrid approach samples posterior distribution more effectively when compared with the conventional approach.

Hamiltonian Monte Carlo (HMC) approach that uses gradient-based proposal distributions [57] and has been effectively applied to Bayesian neural networks [67]. In similar way, Langevin dynamics can be used to incorporate gradient-based stepping with Gaussian noise into the proposal distribution [58]. HMC avoids random walk behaviour using an auxiliary momentum vector and implementing Hamiltonian dynamics where the momentum samples are discarded after later. The samples are hence less correlated and tend converge to the target distribution more rapidly.

Variational inference (VI) provides an alternative approach to MCMC methods to approximate Bayesian posterior distribution [68; 69]. _Bayes by backpropagation_ is a VI method that showed competitive results when compared to stochastic gradient descent and dropout methods used as approximate Bayesian methods [70]. Dropout is a regularisation technique that is relatively simple to implement and involves randomly dropping selected weights in forward-pass operation of backpropagation. This improves the generalization performance of neural networks and has been widely adopted [71]. Gal and Ghahramani [72] presented an approximate Bayesian methodology based on dropout-based regularisation which has been used for other deep learning models such as CNNs [73]. Later, Gal and Ghahramani [74] presented VI-based dropout technique for recurrent neural networks (RNNs); particularly, long-short term memory (LSTM) and gated recurrent unit (GRU) models for language modelling and sentiment analysis tasks.

We argue that the use of dropouts used for regularisation cannot be seen as a an alternative of MCMC sampling directly from the posterior distribution. In this case, we do not know the priors nor know much about the posterior distribution and there is little theoretical rigour, only computational efficacy or capturing noise and uncertainty during model training. Furthermore, in the Bayesian methodology, a probabilistic representation using priors is needed which is questionable in the dropout methodology for Bayesian computation. Given that VI methods are seen as approximate Bayesian methods, we need to invest more effort in directly sampling from the posterior distribution in case of Bayesian deep learning models. This can only be possible if both communities (i.e., statistics and machine learning) are aware about the strengths and weaknesses of MCMC methods for sampling Bayesian neural networks that span hundreds to thousands of parameters, and go orders of magnitude higher when looking at Bayesian deep
learning models. The progress of MCMC for deep learning has been slow due to lack of implementation details, libraries and tutorials that provide that balance of theory and implementation.

In this paper, we present a Python-based tutorial for MCMC methods that covers simple Bayesian linear models and Bayesian neural networks. We provide code in Python with data and instructions that enable their use and extension. We provide results for some benchmark problems showing the strengths and weaknesses of implementing the respective Bayesian models via MCMC with detailed instructions for every code sample in a related Github repository so that it is easy to clone and run. Our tutorial is simple, and relies on basic Python libraries such as numpy. We do not use advanced libraries as the goal is to serve as a go-to document for beginners who have basic knowledge of machine learning models, and need to get hands on experience with MCMC sampling. Hence, this is a code-based computational tutorial with theoretical background. Finally, we highlight the challenges in sampling multi-modal posterior distributions in the case of Bayesian neural networks and shed light on the use of convergence diagnostics.

The rest of the paper is organised as follows. In Section 2, we present a background and literature review of related methods. Section 3 presents the proposed methodology, followed by experiments and results in Section 4. Section 5 provides a discussion and Section 6 concludes the paper with directions of future work.

## 2 Background

### Bayesian inference

We recall that Bayesian methods account for the uncertainty in prediction and decision making via the posterior distribution [75]. Note that the posterior is the conditional probability determined after taking into account the prior distribution and the relevant evidence or data via sampling methods. Thomas Bayes (1702 - 1761) presented and proved a special case of the Bayes' theorem [76; 77] which is the foundation of Bayesian inference. However, it was Pierre-Simon Laplace (1749 - 1827) who introduced a general version of the theorem and used it to approach problems [78]. Figure 1 gives an overview of the Bayesian inference framework that uses data with a prior and likelihood to construct to sample from the posterior distribution. This is the building block of the rest of the lessons that will feature Bayesian logistic regression and Bayesian neural networks.

Bayesian inference estimates unknown parameters using prior information or belief about the variable. Prior information is captured in the form of a distribution. A simple example of a prior belief is a distribution that has a positive real-valued number in some range. This essentially would imply a belief that our result or posterior distribution would likely be a distribution of positive numbers in some range which would be similar to the prior but not the same. If the posterior and prior both follow the same type of distribution, this is known as a conjugate prior [79]. If the prior provides useful information about the variable, rather than very loose constraints, it is known as informative prior. The prior distribution is based on expert knowledge (opinion) and is also dependent on the domain for different types of models [80; 81].

The need for efficient sampling methods to implement Bayesian inference has been a significant focus of research in computational statistics. This is especially true in the case of multi-modal and irregular posterior distributions [82; 83; 29] which tend to dominate in Bayesian neural network problems [10; 84]. MCMC sampling methods are used to update the probability for a hypothesis (proposal \(\Theta\)) as more information becomes available. The hypothesis

Figure 1: The relationship of the likelihood with data and prior distribution is shown for sampling the posterior distribution.
is given by a prior probability distribution that expresses one's belief about a quantity (or free parameter in a model) before some data (**d**) are observed. MCMC methods use sampling to construct the posterior distribution (\(P(\Theta|\textbf{x})\)), iteratively using a proposal distribution, prior distribution \(P(\Theta)\), and likelihood function.

\[P(\Theta|\textbf{x})=\frac{\textbf{P}(\textbf{d}|\boldsymbol{\Theta})\times \textbf{P}(\boldsymbol{\Theta})}{\textbf{P}(\textbf{d})} \tag{1}\]

We note that \(P(\textbf{d}|\boldsymbol{\Theta})\) could be seen as the likelihood distribution in disguise. \(P(\textbf{d})\) is the marginal distribution of the data and is often seen as a normalising constant and ignored. Hence, ignoring it, we can also express the above in this way

\[P(\Theta|\textbf{x})\propto\textbf{P}(\textbf{x}|\boldsymbol{\Theta})\times \textbf{P}(\boldsymbol{\Theta}) \tag{2}\]

The likelihood function is a function of the parameters of a given model provided specific observed data [85]. The likelihood function can be seen as a measure of fit to the data for the proposals which are drawn from the proposal distribution. Hence, from an optimisation perspective, the likelihood function can be seen as a fitness or error function. The posterior distribution is constructed after taking into account the relevant evidence (data) and prior distribution, with the likelihood that considers the proposal and the model. MCMC methods essentially implement Bayesian inference via a numerical approach that marginalizes or integrates over the posterior distribution [86]. Note that probability and likelihood are not the same in the field of statistics, while in everyday language they are used as if they are the same. The term "probability" refers to the possibility of something happening in relation to a given distribution of data. The likelihood refers to the likelihood function that provides a measure of fit in relation to a distribution. The likelihood function indicates which parameter (data) values are more likely than others in relation to a distribution. Further and detailed explanations regarding Bayesian inference and MCMC sampling is given here [87; 88].

### Probability distributions

#### 2.2.1 Gaussian (Normal) Distribution

A normal probability density or distribution, also known as the Gaussian distribution, is described by two parameters, mean (\(\mu\)) around which the distribution is centered and the standard deviation (\(\sigma\)) which describes the spread (sometimes described instead by the variance, \(\sigma^{2}\)). Using these two parameters, we can fit a probability (normal) distribution to data from some source. In similar way, given a probability distribution, we can generate data and this process is known as sampling from the distribution. In sampling the distribution we simply present random data points (uniform) to the distribution and get data that is, in a way, transformed by the distribution. These parameters determine the shape of the probability distribution, e.g., if it is peaked or spread. Note that the normal distribution is symmetrical in nature and caters for negative and positive numbers of real data.

Equation 3 presents the Gaussian distribution probability density function (PDF) for parameters \(\mu\) and \(\sigma\).

\[f(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}exp\left(-\frac{1}{2}\left(\frac{x-\mu}{ \sigma}\right)^{2}\right) \tag{3}\]

We will sample from this distribution in Python via the NumPy library [89] to sample from the various distributions discussed in this tutorial and SciPy library [90] to get a representation of the PDF i.e., the probability distribution. The associated github repository1 contains the code to generate Figures 2 to 5 using the Seaborn and Matplotlib Python libraries.

Footnote 1: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/01-Distributions.ipynb](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/01-Distributions.ipynb)

```
1importnumpyasnp
2fromnumpyimportrandom
3fromscipyimportstats
4#10randomdrawsfromastandardGaussiandistribution-N(0,1)
5samples=random.randn(10)
* Gaussian distribution PDF with mean of 4 and standard deviation of 0.5 x = np.linspace(0, 8, 100) p pdf = stats.norm.pdf(x, loc=4, scale=0.5) ```

Listing 1: Random number generation for a Gaussian distribution

We note that the mean and standard deviation are purely based on the data and it will change depending on the dataset. Let us visualise what happens when standard deviation changes and mean remains the same for the distribution as shown in Figure 2.

We can review some more examples where changes to the mean and standard deviation gives us different shapes of the PDF as shown in Figure 3.

#### 2.2.2 Multivariate Normal distribution

The multivariate normal distribution or joint normal distribution generalises univariate normal distribution to more variables or higher dimensions as shown in the PDF in Equation 4.

\[f(x_{1},\ldots,x_{M})=\frac{1}{\sqrt{(2\pi)^{M}|\Sigma|}}\exp\left(-\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}}\Sigma^{-1}(\mathbf{x}-\boldsymbol {\mu})\right) \tag{4}\]

where \(\mathbf{x}\) is a real \(M\)-dimensional column vector and \(|\Sigma|\) is the determinant of symmetric covariance matrix which is positive definite.

#### 2.2.3 Gamma distribution

A gamma distribution is defined by the parameters shape (\(\alpha\)) and rate (\(\beta\)) as shown below.

Figure 3: Normal distributions with different parameters, i.e mean and the standard deviation.

Figure 2: Normal distributions with the same mean and different standard deviations.
\[f(x;\alpha,\beta)=\frac{\beta^{\alpha}x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)} \tag{5}\]

for \(x>0\quad\alpha,\beta>0\); where \(\Gamma(n)=(n-1)!\). Figure 4 presents the Gamma distribution for various parameter combinations, with corresponding Python code given below.

```
10(size)randomdrawsfromaGammadistributionwithashapeparameterof2andscale(1/rate)of2samples=random.gamma(2,scale=2,size=10)
#GammadistributionPDFwithashapeparameterof2andscale(1/rate)of2
#gammadistributionPDFwithashapeparameterof2andscale(1/rate)of2
#x=np.linspace(0,10,100)
#pdf=stats.gamma.pdf(x,a=2,scale=2)

#InverseGammadistributionPDFwithashapeparameterof2andscaleof1
#x=np.linspace(0,5,100)
#pdf=stats.invgamma.pdf(x,a=2,scale=1)
```

Listing 2: Random number generation for the Gamma and inverse Gamma distributions

The corresponding inverse-Gamma (IG) distribution takes the same parameters with examples given in Figure 5 and is more appropriate for real positive numbers.

#### 2.2.4 Binomial distribution

So far, we have only addressed real numbers with respective probability distributions; however, we need to also consider discrete numbers. The Bernoulli distribution is a discrete probability distribution typically used for modelling binary classification problems. We begin with an example where a variable \(x\) takes the value 1 with probability \(p\) and

Figure 4: Gamma distributions with different shape and rate parameters (\(\alpha\) and \(\beta\)).

Figure 5: Inverse gamma distributions with different shape and rate parameters (\(\alpha\) and \(\beta\)).
the value 0 with probability \(q=1-p\). The probability mass function for this distribution over the possible outcomes (\(x\)) is given by Equation 6.

\[f(x;p)=p^{x}(1-p)^{1-x} \tag{6}\]

for \(x\in{0,1}\). The probability of getting exactly \(k\) successes (\(x=1\)) in \(n\) independent Bernoulli trials (\(f(k,n,p)\)) is given as \(\Pr(k;n,p)\) in Equation 7.

\[\Pr(k;n,p)=\binom{n}{k}p^{k}(1-p)^{n-k} \tag{7}\]

for \(k=0,1,2,...,n\), where \(\binom{n}{k}=\frac{n!}{k!(n-k)!}\). Listing 3 below shows the implementation in Python.

```
10(size)randomdrawsfromaBinomialdistributionwithashapeparameterof2andscale(1/rate)of2
2samples=random.gamma(2,scale=2,size=10)
#BinomialdistributionPDFwithashapeparameterof2andscale(1/rate)of2
3x=np.linspace(0,10,100)
5pdf=stats.gamma.pdf(x,a=2,scale=2)
```

Listing 3: Random number generation for the Binomial distributions

#### 2.2.5 Multinomial distribution

Previously, we catered for the case of two outcomes; however, we can consider the case of more than two outcomes. Suppose a single trial can result in \(k\) (\(k\geq 2\)) possible outcomes numbered \(1,2,\ldots,k\) and let \(p_{i}=\mathbb{P}\)(a single trial results in outcome i) (\(\sum_{i=1}^{k}p_{i}=1\)). For \(n\) independent trials, let \(X_{i}\) denote the number of trials resulting in outcome \(i\) (then \(\sum_{i=1}^{k}X_{i}=n\)). Then we say that the distribution of \((X_{1},X_{2},\ldots,X_{k})\sim\text{Multinomial}(n;p_{1},p_{2},\ldots,p_{k})\) and it holds

\[\mathbb{P}(X_{1}=x_{1},X_{2}=x_{2},\ldots,X_{k}=x_{k})=\frac{n!}{x_{1}!x_{2}! \ldots x_{k}!}p_{1}^{x_{1}}p_{2}^{x_{2}}\ldots p_{k}^{x_{k}},\ 0<p_{i}<1,\ \sum_{i=1}^{k}p_{i}=1. \tag{8}\]

## 3 Mcmc

We begin by noting that a Markov process is uniquely defined by its transition probabilities \(P(x^{\prime}|x)\) which defines the probability of transitioning from any given state \(x\) to another given state \(x^{\prime}\). The Markov process has a unique stationary distribution \(\pi(x)\) given the following two conditions are met.

1. There must exist a stationary distribution \(\pi\) which solves the detailed balance equations, and therefore requires that each transition \(x\to x^{\prime}\) is reversible. This implies that for every pair of states \(x,x^{\prime}\), the probability of being in state \(x\) and moving to state \(x^{\prime}\), must be equal to the probability of being in state \(x^{\prime}\) and moving to state \(x\); hence, \(\pi(x)P(x^{\prime}\mid x)=\pi(x^{\prime})P(x\mid x^{\prime})\).
2. The stationary distribution must be unique which is guaranteed by ergodicity of the Markov process [91; 92; 93]. Ergodicity is guaranteed when every state is aperiodic (i.e., the system does not return to the same state at fixed intervals) and positive recurrent (i.e., the expected number of steps for returning to the same state is finite). An ergodic system is one that mixes well, in other words you get the same result whether you average its values over time or over space.

Given that \(\pi(x)\) is chosen to be \(P(x)\), the condition of detailed balance becomes \(P(x^{\prime}\mid x)P(x)=P(x\mid x^{\prime})P(x^{\prime})\) which is re-written as shown in Equation 9.

\[\frac{P(x^{\prime}\mid x)}{P(x\mid x^{\prime})}=\frac{P(x^{\prime})}{P(x)} \tag{9}\]

Algorithm 1 presents a basic MCMC sampler with random-walk proposal distribution that runs until a maximum number of samples (\(N_{max}\)) is reached for training data, \(\mathbf{d}\).
```
Data: Training data, \(\mathbf{d}\) Result:\(N_{max}\) samples from the posterior distribution - Initialise \(x_{0}\) ; for\(i=1\) until\(N_{max}\)do  1. Propose a value \(x^{\prime}|x_{i}\sim q(x_{i})\), where \(q(.)\) is the proposal distribution;  2. Given \(x^{\prime}\), execute the model \(f(x^{\prime},\mathbf{d})\) to compute the predictions (output \(y\)) and the likelihood;  3. Calculate the acceptance probability \(\alpha=\min\left(1,\frac{P(x^{\prime})}{P(x)}\frac{q(x|x^{\prime})}{q(x^{ \prime}|x_{0})}\right)\)  4. Generate a random value from a uniform distribution \(u\sim U(0,1)\);  5. Accept or reject proposed value \(x^{\prime}\); if\(u<\alpha\)then  1 accept the sample, \(x_{i}=x^{\prime}\) else  2 reject current and retain previous sample, \(x_{i}=x_{i-1}\)  end if  end for
```

**Algorithm 1**A basic MCMC sampler leveraging the Metropolis-Hastings algorithm

Algorithm 1 proceeds by proposing new values of the parameter \(x\) (Step 1) from the selected proposal distribution \(q(.)\), in this case a uniform distribution between 0 and 1. Conditional on these proposed values, the model \(f(x^{\prime},\mathbf{d})\) computes or predicts an output using proposal x' and data \(\mathbf{d}\) (Step 2). We computer the likelihood using the prediction, and employ a Metropolis-Hasting criterion (Step 3) to determine whether to accept or reject the proposal (Step 5). We compare the acceptance ratio \(\alpha\) with \(u\sim U(0,1)\), this enforces that the proposal is accepted with probability \(\alpha\). If the proposal is accepted, the chain moves to this proposed value. If rejected, the chain stays at the current value. The process is repeated until the convergence criterion is met, which in this case is the maximum number of samples (\(N_{max}\)) defined by the user.

### Priors

The prior distribution is generally based on belief, expert opinion or other information without viewing the data [8; 94]. Information to construct the prior can be based on past experiments or the posterior distribution of the model for related datasets. There are no hard rules for how much information should be encoded in the prior distribution; hence, we can take multiple approaches.

An _informative prior_ gives specific and definite information about a variable. If we consider the prior distribution for the temperature tomorrow evening, it would be reasonable to use a normal distribution with an expected value (as mean) of today's evenings temperature with a standard deviation of the temperature each evening for the entire season. A _weakly informative prior_ expresses partial information about a variable. In the case of the prior distribution of evening temperature, a weakly informative prior would consider day time temperature of the day (as mean) with a standard deviation of day time temperature for the whole year. An _uninformative prior_ or _diffuse prior_ expresses vague information about a variable, such as the variable is positive or has some limit range. A number of studies have been done regarding priors for linear models [95; 96] and Bayesian neural networks and deep learning models [97]. Hobbs et al. [98] presented a study for Bayesian priors in generalised linear models for clinical trials. We note that incorporation of prior knowledge in deep learning models [99], is different from selecting or defining priors in Bayesian deep learning models. Due to the similarity of terms, we caution the readers that these can be often confused and mixed up.

### MCMC sampler in Python

We begin with an example with a deliberately simple example in which we are only sampling one parameter from a binomial distribution in order to demonstrate a simple MCMC implementation in Python. Looking at a simple binomial (e.g., coin flipping) likelihood (we will explore the likelihood later) problem and given the data of \(k\) successes in \(n\) trials, we will calculate the posterior probability of the parameter \(p\) which defines the chance of success for any given trial. MCMC sampling requires a prior distribution along with a likelihood function that is used to evaluate a
set of parameters proposed for the given data and model. In other words, the likelihood is a measure of the quality of proposals obtained from a defined proposal distribution using MCMC sampling. We first need to define our prior.

Listing 4 presents an implementation 2 of this simple MCMC sampling exercise in Python of Algorithm 1.

Footnote 2: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/02-Basic-MCMC.ipynb](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/02-Basic-MCMC.ipynb)

In this example, we adopt a uniform distribution as an uninformative prior, only constraining the \(p\) to be between the values of 0 and 1 (\(p\in[0,1]\)).

```
1#Firstdefineourlikelihoodfunctionwhichwillbedependentonprovided"data'
2#inthiscassewwillchoosek=50,n=100
3deflikelihood(query_prob):
4,'
5Giventhedatofksuccessesinntrials,returnalikelihoodfunctionwhich
6evaluatestheprobabilitythatasingsuccess(p)isquery_probforanygiven
7query_prob(between0and1).
8,'
9k=50
10n=100
11returnstats.binnom.pmf(k,n,query_prob)
12
13##MCMCSettingsandSetup
14n_samples=10000#numberofsamplestodrawfromtheposterior
15burn_in=500#numberofsamplestodiscardbeforerecordingdrawfromtheposterior
16
17x=random.uniform(0,1)#initialiseavalueofx0
18
19#createanarrayofNaNstoffillwithoursamples
20p_posterior=np.full(n_samples,np.nan)
21
22print('Generating{}MCMCsamplesfromtheposterior:'.format(n_samples))
23
24#nowweecanstarttheMCMCsamplingloop
25foriiinnp.arange(n_samples):
26#Samplevalueuniformlyfrom0tod1asapprosal
27x_new=random.uniform(0,1)
28
29#CalculateMetropolis-Hastingsacceptanceprobabilitybasedontheprior
30#(canbeignoredinthiscase)andlikelihood
31prior_ratio=1#forthissimpleexamplesdiscussedabove
32likelihood_ratio=likelihood_function(x_new)/likelihood_function(x)
33alpha=np.min([1,likelihood_ratio*prior_ratio])
34
35#Hereweusearandomdrawfromanuniformdistributionbetween0and1asa
36#methodofacceptingthemewproposalwithaprobabilityofalpha
37#(i.e.,acceptifu<alpha)
38u=random.uniform(0,1)
39ifu<alpha:
40x=x_new#thenupdatethecurrentsampletotheproposalforthenextiteration
41
42#Storethecurrentsample
43p_posterior[ii]=x
```

Listing 4: Python implementation of Algorithm 1

In this example, we adopt a uniform distribution as an uninformative prior, only constraining the \(p\) to be between the values of 0 and 1 (\(p\in[0,1]\)).

```
1#Firstdefineourlikelihoodfunctionwhichwillbedependentonprovided"data'
2#inthiscassewwillchoosek=50,n=100
3deflikelihood(query_prob):
4,'
5Giventhedatofksuccessesinntrials,returnalikelihoodfunctionwhich
6evaluatestheprobabilitythatasingsuccess(p)isquery_probforanygiven
7query_prob(between0and1).
8,'
9k=50
10n=100
11returnstats.binnom.pmf(k,n,query_prob)
12
13##MCMCSettingsandSetup
14n_samples=1000#numberofsamplestodrawfromtheposterior
15burn_in=500#numberofsamplestodiscardbeforerecordingdrawfromtheposterior
16
17x=random.uniform(0,1)#initialiseavalueofx0
18
19#createanarrayofNaNstoffillwithoursamples
20p_posterior=np.full(n_samples,np.nan)
21
22print('Generating{}MCMCsamplesfromtheposterior:'.format(n_samples))
23
24#nowecanstarttheMCMCsamplingloop
25foriinnp.arange(n_samples):
26#Samplevalueuniformlyfrom0tod1asapprosal
27x_new=random.uniform(0,1)
28
29#CalculateMetropolis-Hastingsacceptanceprobabilitybasedontheprior
30#(canbeignoredinthiscase)andlikelihood
31prior_ratio=1#forthissimpleexamplesdiscussedabove
32likelihood_ratio=likelihood_function(x_new)/likelihood_function(x)
33alpha=np.min([1,likelihood_ratio*prior_ratio])
34
35#Hereweusearandomdrawfromanuniformdistributionbetween0and1asa
36#methodofacceptingthemewproposalwithaprobabilityofalpha
37#(i.e.,acceptifu<alpha)
38u=random.uniform(0,1)
39ifu<alpha:
40x=x_new#thenupdatethecurrentsampletotheproposalforthenextiteration
41
42#Storethecurrentsample
43p_posterior[ii]=x
```

Listing 4: Python implementation of Algorithm 1

An important note, however, in MCMC, a certain portion of the initial samples are discarded. The discarded samples are known as the burn-in or warmup period. The burn-in can range depending on the sampling problem, but here we will use 25 % to 50 % depending on the complexity of the model. If you use MCMC for large neural network architectures, 50 % burn-in will likely be required. Note that burn-in could be seen as an optimisation stage. Essentially you are discarding material that is not part of the posterior distribution. Your posterior distribution should feature good predictions, and that is what you get after your sampler goes towards convergence.

In order to visualise the MCMC results, typically histograms of the posterior distribution and the trace plot are used. The histogram of the posterior distribution allows us to examine the mean and variance visually, while the
trace plot shows the value of samples at each iteration, allowing us to examine the behaviour and convergence of the MCMC.

Although it is necessary to exclude the burn-in samples in the posterior distribution, it can be helpful to include them in the trace plots so that we can examine where the model started and how well it converged. The following plots provide a visualisation of the results from Listing 4, where a normal distribution is produced by a simple MCMC. Since this is a relatively simple model, a small burn-in proportion is used. The histogram of posterior shows a normally distributed shape, and the trace plot shows that the samples are distributed around the convergence value, as well as the burn-in samples which are in red. We also note that the value of the posterior is usually taken as the mean of the distribution, and in this case our mean value is 0.0413.

## 4 Bayesian linear models via MCMC

We give details of implementing Bayesian logistic regression that uses Metropolis-Hastings MCMC with random-walk proposal distribution.

We wish to model a univariate timeseries consisting of observed data points with inputs \(X=(\mathbf{x}_{1},\ldots,\mathbf{x}_{S})^{\prime}\). We assume that the relationship between inputs and outputs is a signal plus noise model where the signal depends upon a set of parameters \(\theta\), denoted by \(f(\mathbf{x},\theta)\). We assume the noise to be Gaussian with a mean of zero and a variance of \(\tau^{2}\), so that

\[y=f(\mathbf{x},\theta)+e\qquad e\sim\mathcal{N}(0,\tau^{2}) \tag{10}\]

or equivalently,

\[p(y|x,\theta,\tau^{2})\sim\mathcal{N}\left(f(\mathbf{x},\theta),\tau^{2}\right) \tag{11}\]

We have adopt above form to generalise to the neural network model later; hence, Equation 12 expresses the case of a linear model.

\[f(\mathbf{x},\theta)=\theta\mathbf{X}^{T} \tag{12}\]

Figure 6: Posterior and trace plot for the basic MCMC sampler example.
In this linear regression problem, we then have the parameters \(\theta\) (that represents all weights and biases) and \(\tau\) (that represents a single noise parameter). We sample the parameters using MCMC to find their posterior distributions. Note that from an optimisation perspective, the case of sampling can be seen as a form of optimisation, eg. using gradient-based methods [100] for learning the parameters of linear models or neural networks, i.e similar to the case of training a perceptron model [101] in the machine learning and neural networks literature. The additional aspect of a MCMC sampler is their ability to sample a posterior probability distribution that represents the parameters of a model rather than a fixed point estimate given by optimisation methods.

### Likelihood

As detailed above in Section 2.1, our Bayesian treatment of the problem of finding the posterior \(p(\theta\mid y)\) requires the definition of both a likelihood \(p(\theta\mid\mathbf{x})\) and prior distribution \(p(\theta)\). We begin by defining the likelihood. The likelihood, i.e probability of the data given the model is given by the product of the likelihood for each of our data points (\(S\) total) as shown in Equation 13.

\[p(y\mid x,\theta,\tau^{2})=\prod_{t=1}^{S}p(y_{t}\mid x_{t},\theta,\tau^{2}) \tag{13}\]

We note that for MCMC implementation, we would be using log-likelihood (i.e taking the log of the likelihood function) to eliminate numerical instabilities which can occur since we multiple probabilities together which grows with the size of the data. It is also more convenient to maximize the log of the likelihood function since the logarithm is monotonically increasing function of its argument, maximization of the log of a function is equivalent to maximization of the function itself. In order to transform a likelihood function into a log-likelihood, we will use the log product rule as given below.

\[log_{b}(x\times y)=log_{b}(x)+log_{b}(y) \tag{14}\]

The log-likelihood simplifies the subsequent mathematical analysis and also helps avoid numerical instabilities due to the product of a large number of small probabilities. In the log-likelihood, Equation 13 is much simplified by computing the sum of the log probabilities as given in Equation 15.

\[\ln p(y\mid x,\theta,\tau^{2})=\sum_{t=1}^{S}\ln p(y_{t}\mid x_{t},\theta, \tau^{2}) \tag{15}\]

In order to construct the likelihood function, we use our definition of the probability for each data point given the model as shown in Equation 11, and the form of the Gaussian distribution as defined in Equation 3. We use a set of weights and biases as the model parameters \(\theta\) in our model \(f(x,\theta)\) for \(S\) training data instances and variance \(\tau^{2}\). Our assumption of normally distributed errors leads to a likelihood given in Equation 16.

\[p(y\mid x,\theta,\tau^{2})=\frac{1}{(2\pi\tau^{2})^{S/2}}\times\exp\left(- \frac{1}{2\tau^{2}}\sum_{t=1}^{S}(y_{t}-f(x_{t},\theta))^{2}\right) \tag{16}\]

### Prior

We note that a linear model transforms into a Bayesian linear model with the use of a prior distribution, and a likelihood function to sample the posterior distribution via MCMC. In Section 3.2, we discussed the need to define a prior distribution for our model parameters \(\theta\) and \(\tau\). In the case where the prior distribution comes from the same probability distribution family as the posterior distribution, the prior and posterior are then called conjugate distributions [102; 103]. The prior is called a conjugate prior for the likelihood function of the Bayesian model.

To implement conjugate priors in our linear model, we will assume a multivariate Gaussian prior for _theta_ (Equation 17) and an inverse Gamma distribution (IG) for \(\tau^{2}\) (Equation 18). We use an IG prior distribution since we require positive real numbers to represent \(\tau^{2}\). We use the multivariate Gaussian distribution to represent the parameters such as weights and bias of the linear models. These parameters feature negative and positive real number values and we
are involving a model with more than one parameter, hence the multivariate Gaussian distribution is most appropriate for the prior.

\[\theta\sim\mathcal{N}(0,\sigma^{2}) \tag{17}\]

\[\tau^{2}\sim IG(\nu_{1},\nu_{2}) \tag{18}\]

In this example, we choose parameter values of \(\sigma=5\), \(\nu_{1}=0\), and \(\nu_{2}=0\) to define the prior distributions. These values are based on expert opinion who consider trained models.

First, lets revisit multivariate normal distribution from Equation 4, used here to define the prior distribution for our weights and biases. Suppose that our \(\mathbf{x}\) is our set of \(M\) weights and biases, \((\theta_{1},\ldots,\theta_{M})\). Suppose our \(\mu\) is a vector of zeros (as defined above), then we get:

\[f(\theta_{1},\ldots,\theta_{M})=\frac{1}{\sqrt{(2\pi)^{M}[\Sigma]}}\exp\left(- \frac{1}{2}(\mathbf{\theta})^{\mathrm{T}}\Sigma^{-1}(\mathbf{\theta})\right) \tag{19}\]

The covariance matrix \(\mathbf{\Sigma}\) in this case is just a diagonal matrix with all values equal to \(\sigma^{2}\) (scalar), \(\mathbf{\Sigma}^{-1}\) will become \(I/\sigma^{2}\) where \(I\) is an identity matrix (diagonal elements which are all ones). Hence, the numerator in above equation

\[(\mathbf{\theta})^{\mathrm{T}}\mathbf{\Sigma}^{-1}(\mathbf{\theta}) \tag{20}\]

becomes

\[\frac{(\mathbf{\theta})^{\mathrm{T}}\mathbf{I}(\mathbf{\theta})}{\sigma^{2}} \tag{21}\]

We note that multiplying identity matrix with any other matrix is the matrix itself, hence finally we get \(\theta^{2}\) in numerator. We can now move to the IG distribution used to define the prior for our model variance (\(\tau^{2}\)).

\[f(\tau^{2})=\frac{\nu_{1}^{\nu_{2}^{\nu_{2}}}}{\Gamma(\nu_{1})}\left(\frac{1} {\tau^{2}}\right)^{\nu_{1}+1}\exp\left(\frac{-\nu_{2}}{\tau^{2}}\right) \tag{22}\]

We note that \(\nu_{1}^{\nu_{2}}/\Gamma(\nu_{1})\) is a constant which can be dropped considering proportionality.

Equation 23 takes into account the product of all our MCMC parameters to define the overall prior.

\[p(\mathbf{\theta})\propto\frac{1}{(2\pi\sigma^{2})^{M/2}}\times\exp\bigg{\{}-\frac {1}{2\sigma^{2}}\bigg{(}\sum_{i=1}^{M}\theta^{2}\bigg{)}\bigg{\}}\times\tau^{ -2(1+\nu_{1})}\exp\left(\frac{-\nu_{2}}{\tau^{2}}\right) \tag{23}\]

### Python Implementation

The following code presented in Listing 13 implements a model in the form outlined above. First, we will define our simple linear model as in Equation 12.

```
1classLinearModel:
2'''
3Simplelinearmodelwithasingoutput(y)giventhecovariatesx_1...x_Moftheform:
4y=w_1*x_1+...+w_M*x_M+b
5whereM=numberoffeatures,waretheweights,andbisthebias.
6'''
7#Initialisevaluesofmodelparameters
8def_init...(self):
9self.w=None
10self.b=None
11
* Functiontotakeindataandparametersampleandreturntheprediction
*default_proposal(self,data,theta):
*"
*Encodetheproposedparametersandthenusethemodeltoppredict
*Input:
*data:(NxM)arrayofdata
*theta:(M+1)vectorofparameters.Thelastelementoftheatsomitesthebisterm(givingM+1elements)
*";"
*self.encode(theta)#methodotencodewandb
*prediction=self.predict(data)#predictandreturn
*returnprediction
*"Linearmodelprediction
*defpredict(self,x_in):
*y_out=x_in.dot(self.w)+self.b
*returny_out
*#Helperfunctiontsplittheparametervectorintowandbandstoreinthemodel
*defencode(self,theta):
*self.w=theta[0:-1]
*self.b=theta[-1]

```
12#Functiontotakeindataandparametersampleandreturntheprediction
*defevaluate_proposal(self,data,theta):
*"
*Encodetheproposedparametersandthenusethemodeltoppredict
*Input:
*data:(NxM)arrayofdata
*theta:(M+1)vectorofparameters.Thelastelementoftheatsomitesthebisterm(givingM+1elements)
*log_likelihood:loglikelihoodofthedatagiventheparameters
*model_prediction:predictionofthemodelgiventheparameters
*accuracy:accuracy(RMSE)ofthemodelgiventheparameters
*","
*#firstmakeapplicationwithparameterstheta
*model_prediction=self.model.evaluate_proposal(self.x_data,theta)
*accuracy=self.rmse(model_prediction,self.y_data)#RMSEerrormetric
*nowcalculatetheloglikelihood
*log_likelihood=np.sum(-0.5*np.log(2*np.pi*tausq)-0.5*np.square(self.y_data-model_prediction)/tausq)
*return[log_likelihood,model_prediction,accuracy]
*Definetheprior
*defprior_likelihood(self,sigma_squared,nu_1,nu_2,theta,tausq):
*"
*Calculatethepriorlikelihoodoftheparameters
*
```
24Input:
25sigma_squared: varianceofnormalpriorfortheta
26nu_1:parameternu_1ofttheinversegammapriorfortau^2
27nu_2:parameternu_2ofttheinversegammapriorfortau^2
28theta:(M+1)vectorofparameters.Thelastelementofthetaconstitutesthebias term(givingM+1elements)
29tausq:varianceoftheerrorterm
30Output:
31log_prior:logpriorlikelihood
32,,
33n_params=self.theta_size#numberofparametersinmodel
34part1=-1*(n_params/2)*np.log(sigma_squared)
35part2=1/(2*sigma_squared)*(sum(np.square(theta)))
36log_prior=part1-part2-(1+nu_1)*np.log(tausq)-(nu_2/tausq)
37returnlog_prior
```

Listing 6: Python implementation of likelihood and prior functions for linear regression model to be incorporated into the MCMC sampling class

Before running MCMC, we need to set up the sampler. First we generate an initial sample for our parameters, and initialise arrays to capture the samples that form the posterior distribution, the accuracy, and the model predictions. Then we proceed with sampling as per the MCMC sampling algorithm detailed in Algorithm 1. This algorithm uses a Gaussian random walk for the parameter proposals (\(\theta_{p}\) and \(\tau_{p}^{2}\)), perturbing the previous proposed value with Gaussian noise as shown in Equations 26 and 27, respectively.

\[\theta_{p}\sim\theta_{p-1}+\mathcal{N}(0,\sigma_{\theta}) \tag{26}\]

\[\tau_{p}^{2}\sim\tau_{p-1}^{2}+\mathcal{N}(0,\sigma_{\tau^{2}}) \tag{27}\]

We then accept/reject the proposed value according to the Metropolis-Hastings acceptance ratio. Note that the log-likelihood is used and hence the ratio of previous and current likelihood will need to consider log laws (rules), i.e we note the log product rule in Equation 28 and the quotient rule in Equation 29.

\[log_{b}(x\times y)=log_{b}(x)+log_{b}(y) \tag{28}\]

\[log_{b}(x/y)=log_{b}(x)-log_{b}(y) \tag{29}\]

```
1#MCMCsampler
2defsampler(self):
3'''
4Runthesamplerforadefinedlinearmodel
5'''
6##Definemeptryarraystotorethesampledposteriorvalues
7#posteriorofallweightsandbiasoverallsamples
8pos_theta=np.ones((self.n_samples,self.theta_size))
9#posteriordefiningthevarianceofthenoiseinpredictions
10pos_tau=np.ones((self.n_samples,1))
11
12#recordoutputf(x)overallsamples
13pred_y=np.ones((self.n_samples,self.x_data.shape[0]))
14#recordtheRMSEofeachsample
15rmse_data=np.zeros(self.n_samples)
16
17#Initialisation
18#initialistheta-themodelparameters
19theta=np.random.randn(self.theta_size)
20#makeinitialprediction
21pred_y[0,]=self.model.evaluate_proposal(self.x_data,theta)
22
23#initialiseeta-wesampleetaasagaussianrandomwalkinthelogspaceoftau^2
24eta=np.log(np.var(pred_y[0,]-self.y_data))
tausq_proposal=np.exp(eta)
2
3#calculatethepriorlikelihoodprior_likelihood=self.prior_likelihood(self.sigma_squared,self.nu_1,self.nu_2,theta,tausq_proposal)
4#calculatethelikelihoodconsideringobservations
5[likelihood,pred_y[0,],_,rmse_data[0]]=self.likelihood_function(theta,tausq_proposal)
6
7n_accept=0
8##RuntheMCMCsampleforn_samples
9foriiminp.arange(1,self.n_samples):
10#SamplenewvaluesfortheatandtauusingaGaussianrandomwalk
11theta_proposal=theta+np.random.normal(0,self.step_theta,self.theta_size)
12eta_proposal=eta+np.random.normal(0,self.step_eta,1)#sampletau^2inlogspace
13tausq_proposal=np.exp(eta_proposal)
14
15#calculatethepriorlikelihoodprior_proposal=self.prior_likelihood(self.sigma_squared,self.nu_1,self.nu_2,theta_proposal,tausq_proposal)
16)
17#calculatethelikelihoodconsideringobservations
18[likelihood_proposal,pred_y[ii,],_,rmse_data[ii]]=self.likelihood_function(theta_proposal,tausq_proposal)
19#Notingthatlikelihood_functionandprior_likelihoodreturnloglikelihoods,
20#wecanuseglaustocalculatethemacetometaprobolity
21diff_likelihood=likelihood_proposal-likelihood
22diff_priorlikelihood=prior_proposal-prior_likelihood
23
24mh_prob=min(1,mp.exp(diff_likelihood+diff_priorlikelihood))
25
26#sampletocaptorrejecttheproposalaccordingtotheacceptanceprobability
27u=np.random.uniform(0,1)
28ifu<mh_prob:
29#acceptandupdatethevalues
30n_accept+=1
31likelihood=likelihood_proposalprior_likelihood=prior_proposaltheta=theta_proposal
32eta=eta_proposal
33#storetomakeuptheposterior
34pos_theta[ii,]=theta_proposal
35pos_theta[ii,]=theta_proposal
36pos_theta[ii,]=theta_proposal
37pos_tau[ii,]=theta_proposal
38else:
39#rejectmoveandstoretheoldvalues
40pos_theta[ii,]=pos_theta[ii-1,]
41pos_tau[ii,]=pos_tau[ii-1,]
42
43#calculatetheacceptancratesascheck
44accept_rate=(n_accept/self.n_samples)*100 print('{:.3f}%wereaccepted'.format(accept_rate))
45
46
47#storetheposterior(samplesafterburnin)inapandasdataframeandreturn
48self.pos_theta=pos_theta[self.n_burnin;,]
49self.pos_tau=pos_tau[self.n_burnin;,]
40self.rmse_data=rmse_data[self.n_burnin;]
41
42#splitthetaintowandb
43results_dict={'w()'.format(_):self.pos_theta[:,_].squeeze()for_inrange(self.theta_size-1)}
44results_dict['b']=self.pos_theta[:,-1].squeeze()
45results_dict['tau']=self.pos_tau.squeeze()
results_dict['rmse']=self.rmse_data.squeeze()
*results_df=pd.DataFrame.from_dict(
*results_dict
*)
*returnresults_df ```
Listing 7: Python implementation of an MCMC sampler for the linear model

Now that we have the sampler, we can create an MCMC class that brings together the model, data, hyperparameters, and sampling algorithm.
```
1classMCMC:
2def__init__(self,n_samples,n_burnin,x_data,y_data):
3self.n_samples=n_samples#numberofMCMCsamples
4self.n_burnin=n_burnin#numberofburn-insamples
5self.x_data=x_data#(NxM)
6self.y_data=y_data#(Nx1)
7self.theta_size=x_data.shape[1]+1#weightsforeachfeatureandabiasterm(M+1)
8
9#MCMCparameters-definesthevarianceterminourGaussianrandomwalk
10self.step_theta=0.02;
11self.step_eta=0.01;#noteetaisusedastuinthesamplertoconsiderlog scale.
12
13#modelhyperparameters
14#consideredbylookingatdistributionofsimilartrainedmodels-i.edistributionofweightsandbias
15self.sigma_squared=5
16self.nu_1=0
17self.nu_2=0
18
19#initialisethelinearmodelclass
20self.model=LinearModel()
21
22#storeoutput
23self.pos_theta=None
24self.pos_tau=None
25self.rmse_data=None
26
27#functionsdefinedabove-thisispoorpractice,butdoneforreadability
28#andclarity
29self.likelihood_function=MethodType(likelihood_function,self)
30self.prior_likelihood=MethodType(prior_likelihood,self)
31self.sampler=MethodType(sampler,self)
32
33defmodel_draws(self,num_samples=10):
34'''
35Simulatenewmodelpredictions(mu)undertheassumptionthatourposteriorsare
36Gaussian.
37,
38#num_samplesxnum_data_points
39pred_y=np.zeros((num_samples,self.x_data.shape[0]))
40sim_y=np.zeros((num_samples,self.x_data.shape[0]))
41
42foriiinrange(num_samples):
43theta_drawn=np.random.normal(self.pos_theta.mean(axis=0),self.pos_theta.std(axis=0),self.theta_size)
44tausq_drawn=np.random.normal(self.pos_tau.mean(),self.pos_tau.std())
45
46[_,pred_y[ii,:],sim_y[ii,:],_]=self.likelihood_function(theta_drawn,tausq_drawn
47)
48returnpred_y,sim_y
* Additional error metric
* targets) ** 2).mean()) ```
Listing 8: Python implementation of an MCMC class for fitting a linear regression model

We can then run the MCMC sampler.
```
1#MCMCSettingsandSetup
2n_samples=20000#numberofsamplestodrawfromtheposterior
3burn_in=int(n_samples*0.25)#numberofsamplestodiscardbeforerecordingdrawsfromtheposterior
4
5#Generatedotadata
6n_data=100
7n_features=1
8x_data=np.repeat(np.expand_dims(np.linspace(0,1,n_data),axis=-1),n_features,axis=1)
9y_data=3*x_data[:,0]+4+np.random.randn(n_data)*0.5
10
11#InitialistheMCMCclass
12mcmc=MCMC(n_samples,burn_in,x_data,y_data)
13#Runthesampler
14results=mcmc.sampler()
15#Drawsamplemodelsfromtheposterior
16pred_y_sim_y=mcmc.model_draws(num_samples=100) ```

Listing 9: Code to call the MCMC sampling class and fit a model to some toy data

After the code runs, we can see the predictions of the trained Bayesian linear model and simulated values (model fit with added variance according to \(\tau^{2}\)), such as below in Figure 7.

## 5 Bayesian neural networks via MCMC

### Neural networks

We utilise a simple neural network, also known as a multilayer perceptron to demonstrate the process of training a Bayesian neural network via MCMC. A neural network model \(f(x)\) is made up of a series of lower level computations, which can be used to transform inputs to their corresponding outputs \(\{\mathbf{\tilde{x}_{t}},\mathbf{y_{t}}\}\). Neural networks are structured as layers of neurons whose value is determined based on a linear combination of inputs from the previous layer, with an activation function used to introduce nonlinearity. In this case, we consider a neural network with one hidden layer and one output, as shown in Figure 8.

Figure 7: Bayesian linear regression model fitting a toy dataset.
As an example, we can calculate the output value of the jth neuron in the first hidden layer of a network \((h,j)\) using a weighted combination of the \(m\) inputs (\(\bar{x}_{t}\)) as shown in Equation 30.

\[g\Big{(}\delta_{h,j}+\sum_{i=1}^{m}w_{i,j}\bar{x}_{t,i}\Big{)} \tag{30}\]

where, the bias (\(\delta_{h,j}\)) and weights (\(w_{i}\) for each of the \(m\) inputs) are parameters to be estimated (trained), and \(g(.)\) is the activation function that can be used to perform a nonlinear transformation. In our case, the function \(g(.)\) is the sigmoid activation function, used for the hidden and output layers of the neural network model with a single hidden layer as shown in Figure 8.

We train the model to approximate the function \(f\) such that \(f(\mathbf{\bar{x}_{t}})=\mathbf{y_{t}}\) for all pairs. We extend our previous calculation a single neuron in the hidden layer to calculate the output \(f(\mathbf{\bar{x}_{t}})\) as shown in Equation 31.

\[f(\mathbf{\bar{x}_{t}})=g\Big{(}\delta_{o}+\sum_{j=1}^{H}v_{h}\times g\Big{(} \delta_{h}+\sum_{i=1}^{m}w_{i,j}\bar{x}_{t,i}\Big{)}\Big{)}, \tag{31}\]

where \(H\) is the number of neurons in the hidden layer, \(\delta_{o}\) is the bias for the output, and \(v_{h}\) are the weights from the hidden layer to the output neuron. The complete set of parameters for the neural network model, shown in Figure 8, is made up of \(\theta=(\mathbf{\bar{w}},\mathbf{\bar{v}},\boldsymbol{\delta})\), where \(\boldsymbol{\delta}=(\boldsymbol{\delta}_{o},\boldsymbol{\delta}_{h})\). \(\mathbf{\bar{w}}\) are the weights transforming the input to hidden layer. \(\mathbf{\bar{v}}\) are the weights transforming the hidden to output layer. \(\boldsymbol{\delta}_{h}\) is the bias for the hidden layer, and \(\boldsymbol{\delta}_{o}\) is the bias for the output layer.

Figure 8: Simple neural network with a single hidden layer. The information is passed and processed from the input to hidden and then finally to the output layer as shown.
### Bayesian neural networks

A Bayesian neural network is a probabilistic implementation of a standard neural network with the key difference being that the weights and biases are represented via the posterior probability distributions rather than single point values as shown in Figure 9. Similar to canonical neural networks [104], Bayesian neural networks also have universal continuous function approximation capabilities. However, the posterior distribution of the network parameters allows uncertainty quantification on the predictions.

The task for MCMC sampling is to learn (sample) the posterior distributions representing the weights and biases of the neural network that best fit the data. As in the previous examples, we begin inference with prior distributions over the weights and biases of the network and use a sampling scheme to find the posterior distributions given training data. Since non-linear activation functions exist in the network, the conjugacy of prior and posterior is lost, and therefore we must employ an MCMC sampling scheme and make assumptions about the distribution of errors.

We specify the model similar to the Bayesian linear regression, assuming a Gaussian error as given in Equation 32.

\[y=f(\mathbf{x},\theta)+e\qquad e\sim\mathcal{N}(0,\tau^{2}) \tag{32}\]

This leads to the same likelihood function as presented in logarithmic form in Equation 24. As in Section 4.2, we adopt Gaussian priors for all parameters of the model (\(\theta\)), with zero mean and a user-defined variance (\(\sigma^{2}\)), and an IG distribution for the variance of the error model (\(\tau^{2}\)), with parameters \(\nu_{1}\) and \(\nu_{2}\). The _likelihood function_ and _prior likelihood_ function, therefore remain unchanged from their definition in Listing 6.

### Training neural networks via backpropagation

We note that typically random-walk proposal distributions are used for small scale-models; however, neural network models features a large number of parameters. The choice of a proposal distribution is essential for models with large number of parameters. To better guide the sampling, we can incorporate gradients to our proposals and we will start by examining how gradients are incorporated in non-probabilistic deep learning problems.

Gradient based methods are widely used in deep learning alongside backpropagation to train models with various neural network architectures [106]. Stochastic gradient descent (SGD) involves stepping through the parameter space

Figure 9: Bayesian neural network and MCMC sampling. Adapted from [105].
iteratively, guided by gradients, to optimize a differentiable objective function. The method has prominently featured in the training of various neural networks architectures including deep learning models. SGD is combined with backpropagation when training neural networks [107; 108; 109]. Backpropagation involves a forward pass which propagates information forward to get prediction (decision) at the output layer and a backward pass to compute the local gradients for each of the parameters. These gradients are then used to inform the update of the model parameters in an iterative process, where the weights and biases are updated at each step. Training neural networks also can be considered as solving a non-convex optimization problem \(\text{argmin}\,L(\theta)\); where, \(\theta\in R^{n}\) is the set of parameters and \(L\) is the loss function. The parameter (weight) update for an iteration (epoch) of SGD is shown in Equation 33.

\[\theta_{k}=\theta_{k-1}-a_{k-1}\nabla L(\theta_{k-1}) \tag{33}\]

where, \(\theta_{k}\) denotes the \(k^{th}\) iteration, \(a_{k}\) is the learning rate, and \(\nabla L(\theta_{k})\) denotes the gradient.

We note that the learning rate is user defined hyperparameter which depends on the problem and data at hand. It is typically determined through tuning using cross-validation or trial and error. A number of extensions of the backpropagation algorithm employing SGD have been proposed to address limitations. These include use of weight decay regularization during training to improve generalisation ability [110], a momentum mechanism for faster training [111; 112], adaptive learning rate [112], and second order gradient methods [113] which, although efficient, have problems in scaling up computationally with larger models. In the last decade, with the deep learning revolution, further attempts have been made to ensure that enhanced backpropagation algorithms, not only improve training accuracy, but can also scale well computationally for large deep learning models. Hence, the development of methods such as the adaptive gradient algorithm (AdaGrad) [114], ada-delta [115], and Adam (adaptive moment estimation) [116]. These advanced algorithms are generally based on the idea of adapting the learning rate automatically during the training, taking into account the recent history of the optimisation process.

#### 5.3.1 Langevin gradient-based proposal distribution

We mentioned earlier that random-walk proposal distributions are suited for small scale models and better proposal distributions would be required for neural network models. Although simple neural networks have much lower number of parameters, when compared to deep learning models, training simple neural networks with MCMC sampling is a challenge with random-walk proposal distribution. We need to utilise the properties of backpropagation algorithm and the mechanism of weight update using gradients. Hence, we utilise _stochastic gradient Langevin dynamics_[58] for the proposal distribution which features the addition of noise to the stochastic gradients. The method has shown to be effective for linear models [58] which motivated its use in Bayesian neural networks. In our previous work, Langevin-gradient MCMC has been very promising for simple and deep neural networks [59; 66; 65]. Hence, we draw the proposed values for the parameters (\(\mathbf{\theta}^{p}\)) according to a one-step (epoch) gradient as shown in Equation 34.

\[\mathbf{\theta}^{p}\sim\mathcal{N}(\mathbf{\bar{\theta}}^{[s]},\Sigma_{\theta}) \tag{34}\]

A Gaussian distribution with a standard deviation of \(\Sigma_{\theta}\), and mean (\(\mathbf{\bar{\theta}}^{[s]}\)) calculated using a gradient based update (Equation 35 of the parameter values from the previous step (\(\mathbf{\theta}^{[s]}\)).

\[\mathbf{\bar{\theta}}^{[s]}=\mathbf{\theta}^{[s]}+r\times\nabla E(\mathbf{\theta}^{[s]}) \tag{35}\]

with learning rate \(r\) (a user-defined hyperparameter) and gradient update (\(\nabla E(\mathbf{\theta}^{[s]})\)) according to the model residuals.

\[E(\mathbf{\theta}^{[s]})=\sum_{t\in\mathcal{T}}(y_{t}-F(\mathbf{x}_{i},\mathbf{\theta} ^{[s]}))^{2} \tag{36}\]

\[\nabla E(\mathbf{\theta}^{[s]})=\left(\frac{\partial E}{\partial\theta_{1}}, \cdots,\frac{\partial E}{\partial\theta_{L}}\right).\]

Hence, the Langevin gradient-based proposal consists of 2 parts:

1. Gradient descent-based weight update
2. Addition of Gaussian noise from \(\mathcal{N}(0,\Sigma_{\theta})\)
We need to ensure that the detailed balance is maintained since the Langevin-gradient proposals are not symmetric. We note that MCMC implementation with relaxed detailed balance condition for some applications also exist [117]. Therefore, a combined update is used as a proposal in a Metropolis-Hastings step, which accepts the proposal \(\mathbf{\theta}^{p}\) for a position \(s\) with the probability \(\alpha\) as shown in Equation 37.

\[\alpha=\min\left\{1,\frac{p(\mathbf{\theta}^{p}|\mathbf{y})q(\mathbf{\theta}^{|s|}|\mathbf{ \theta}^{p})}{p(\mathbf{\theta}^{|s|}|\mathbf{y})q(\mathbf{\theta}^{p}|\mathbf{\theta}^{|s |})}\right\} \tag{37}\]

where, \(p(\mathbf{\theta}^{p}|\mathbf{y})\) and \(p(\mathbf{\theta}^{|s|}|\mathbf{y})\) can be computed using the likelihood from Equation (16), and priors given by Equation (23). The ratio of the proposed and the current \(q(\mathbf{\theta}^{p}|\mathbf{\theta}^{|s|})\), is given by Equation 38 which is based on a one-step (epoch) gradient \(\nabla E_{\mathbf{y}}[\mathbf{\theta}^{|s|}\) with user-defined learning rate \(r\) given in Equation 39.

\[q(\mathbf{\theta}^{|s|}|\mathbf{\theta}^{p})\sim N(\bar{\mathbf{\theta}}^{|s| },\Sigma_{\mathbf{\theta}}) \tag{38}\] \[\bar{\mathbf{\theta}}^{|s|}=\mathbf{\theta}^{|s|}+r\times\nabla E_{ \mathbf{y}}[\mathbf{\theta}^{|s|}] \tag{39}\]

Thus, this ensures that the detailed balance condition holds and the sequence \(\mathbf{\theta}^{|s|}\) converges to draws from the posterior \(p(\mathbf{\theta}|\mathbf{y})\). Since our implementation is in the log-scale, the log posterior is given by 40.

\[\log\left(p(\mathbf{\theta}|\mathbf{y})\right)=\log\left(p(\mathbf{\theta})\right)+ \log\left(p(\mathbf{y}|\mathbf{\theta})\right)+log(q(\mathbf{\theta}|\mathbf{\theta}^{ \star})) \tag{40}\]

Algorithm 2 gives a full description of the Langevin-gradient MCMC sampling scheme, bringing together all the steps outlined above. User defined parameters are present for the maximum number of samples (\(S_{max}\)), rate of Langevin-gradient proposals (\(L_{prob}\)) compared to simple random walk proposals, and learning rate \(r\) for Langevin-gradient proposals. We note that in a standard Langevin MCMC approach, the rate of it being used (\(L_{prob}\)) is 1 and Gaussian noise is already part of Langevin gradients. However, in our implementation, we use a combination of random-walk proposal distribution with Langevin gradients, as this seem to be computationally more efficient than Langevin gradients alone which consume extra computational time when computing gradients, especially in larger models.

We begin by drawing initial values for the \(\theta\) from the prior distribution given in Equation (23) (Stage 1.1). A new proposal is drawn for \(\theta^{p}\) (which incorporates the model weights and biases and \(\tau^{2}\) from either a Langevin-gradient or random-walk proposal distribution (Stage 1.2). The proposal is evaluated using the Bayesian neural network (BNN) model with the likelihood function in Equation 16 (Stage 1.4). We also evaluate the likelihood for the prior given in Equation (23) (Stage 1.3). Using these likelihoods, we can check if the proposal should be accepted using Metropolis-Hastings condition (Stage 1.5 and 1.6). If accepted, the proposal becomes part of the chain, else we retain previous the last accepted state of the chain. We repeat the procedure until the maximum samples are reached (\(S_{max}\)). Finally, we execute post-sampling stage where we obtain the posterior distribution by concatenating the history of the samples in the chain.

### Python Implementation

We first define and implement the simple neural network module (class), and implement methods (functions) for the forward and backward pass to calculate the output of the network given a set of inputs. We need to compute the gradients, and update the model parameters given a model prediction and observations, respectively. Listings 10, 11 and 12 present an implementation4 of the Bayesian Neural Network and associated Langevin-gradient MCMC sampling scheme.

Footnote 4: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/04-Bayesian-Neural-Network.ipynb](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/blob/main/04-Bayesian-Neural-Network.ipynb)

```
#NNprediction

defwardpass(self,X):
```
*TakeaninputXandreturntheoutputofthenetwork
*Input:
*-X:(Nxnum_features)arrayofinputdata
*Output:
Data: Dataset Result: Posterior distribution of model parameters (weights and biases) - Stage 1.0: Metropolis Transition - 1.1 Draw initial values \(\theta_{0}\) from the prior for each s until Smaxdo 1.2 Draw \(\kappa\) from a Uniform-distribution [0,1] if\(\kappa\leq L_{prob}\)then  Use Langevin-gradient proposal distribution: \(\mathbf{\theta}^{p}\sim\mathcal{N}(\bar{\mathbf{\theta}}^{[s]},\Sigma_{\theta})\)  end if else  Use random-walk proposal distribution: \(\mathbf{\theta}^{p}\sim\mathcal{N}(\mathbf{\theta}^{[s]},\Sigma_{\theta})\)  end if 1.3 Evaluate prior given in Equation 23 1.4 Evaluate log-likelihood given in Equation 16 1.5 Compute the posterior probability for Metropolis-Hastings condition - Equation 40 1.6 Draw \(u\) from a Uniform-distribution [0,1] if\(\log(u)\leq log(p(\mathbf{y}|\theta))\)then  Accept replica state: \(\theta^{[s+1]}\leftarrow\theta^{p}\)  end if else  Reject and retain previous state: \(\theta^{[s+1]}\leftarrow\theta^{[s]}\)  end if  end while ```
**Algorithm 2**Bayesian neural network framework using Langevin-based MCMC sampling.
``` -self.12output:(N)arrayofoutputdataf(x)whichcanbe -comparedtoobservations(Y)... #Hiddenlayer 11_z=np.dot(X,self.11_weights)+self.11_biases self.11_output=self.sigmoid(11_z)#activationfunctiong(.) #Outputlayer 12_z=np.dot(self.11_output,self.12_weights)+self.12_biases self.12_output=self.sigmoid(12_z) returnself.12_output
18
19defbackward_pass(self,X,Y):
20'''
21ComputehegradientsusingabackwardpassandundertakeLangevin-gradientupdatingofparameters
22Input:
23-X:(Nxnum_features)arrayofinputdata
24-Y:(N)arrayoftargetdata
25,...
26#dE/dtheta
2712_delta=(Y-self.l2_output)*(self.l2_output*(1-self.l2_output))
28l2_weights_delta=np.outer( self.l1_output,
29l2_delta
30)
31)
32#backpropofl2_deltaandsamesasabove
33l1_delta=np.dot(l2_delta,self.l2_weights.T)*(self.l1_output*(1-self.l1_output))
34l1_weights_delta=np.outer( X,
35l1_delta
36) ```

**Algorithm 3**Bayesian neural network framework using Langevin-based MCMC sampling.
* #updateforoutputlayer
* self.12_weights+=self.rate*12_weights_delta
* self.12_biases+=self.rate*12_delta
* updateforhiddenlayer
* self.11_weights+=self.rate*11_weights_delta
* self.11_biases+=self.rate*11_delta

Listing 10: Python implementation of Neural Network forward and backward passes

Next, we implement the model for a single hidden layer neural network with multiple input neurons and a single output neuron (one step ahead prediction).

```
1classNeuralNetwork:
2'''
3NeuralNetworkmodelwithasinglehiddenlayerandasingleoutput(y)
4'''
5def_init__(self,layer_sizes,learning_rate=0.01):
6'''
7Initializethemodel
8Input:
9-layer_sizes(input,hidden,output):arrayspecifyingthenumberof
10nodesineachlayer
11-learning_rate:learningrateforthegradientupdate
12'''
13#Initialvaluesofmodelparameters
14self.input_num=layer_sizes[0]
15self.hidden_num=layer_sizes[1]
16self.output_num=layer_sizes[2]
17
18#totalnumberofparametersfromweightsandbiases
19self.n_params=(self.input_num*self.hidden_num)+(self.hidden_num*self.output_num)+\
20self.hidden_num+self.output_num
21#learningparams
22self.rate=learning_rate
23
24#Initializetworkstructure
25self.initialise_network()
26
27#functionsdefinedabove-thisispoorpractice,butdoneforreadability
28#andclarity
29self.forward_pass=MethodType(forward_pass,self)
21self.backward_pass=MethodType(backward_pass,self)
23
24definitialize_network(self):
31'''
32Initializetworkstructure-weightsandbiasesforthehiddenlayer
33andoutputlayer
34'''
35#hiddenlayer
36self.11_weights=np.random.normal(
37loc=0,scale=1/np.sqrt(self.input_num),
38size=(self.input_num,self.hidden_num))
39self.11_biases=np.random.normal(
40loc=0,scale=1/np.sqrt(self.hidden_num),
41size=(self.hidden_num,))
42#placeplaceforstoringthehiddenlayervalues
43self.11_output=np.zeros((1,self.hidden_num))
44
45#outputlayer
46self.12_weights=np.random.normal(
47loc=0,scale=1/np.sqrt(self.hidden_num),
48size=(self.hidden_num,self.output_num))
49self.12_biases=np.random.normal(
50loc=0,scale=1/np.sqrt(self.hidden_num),
51size=(self.output_num,))
#placeholderforstoringthemodeloutputs self.12_output=np.zeros((1,self.output_num)) defevaluate_proposal(self,x_data,theta): ''' Ahelperfunctiontotaketheinputdataandproposedparametersample andreturntheprediction Input: data:(Nxnum_features)arrayofdata theta:(w,v,b_h,b_o)vectorofparameterswithweightsandbiases ''' self.decode(theta)#methododecodewintoW1,W2,B1,B2. size=x_data.shape[0] fx=np.zeros(size) foriinrange(0,size):#tosewhatfxisproducedbyyourcurrentweight update fx[i]=self.forward_pass(x_data[i,]) returnfx deflangevin_gradient(self,x_data,y_data,theta,depth): ''' ComputetheLangevingradientbasedproposaldistribution Input: -x_data:(Nxnum_features)arrayofinputdata -y_data:(N)arrayoftargetdata -theta:(w,v,b_h,b_o)vectorofproposedparameters. -depth:SGDdepth Output: -theta_updated:Updatedparameterproposal ''' self.decode(theta)#methododecodewintoW1,W2,B1,B2. size=x_data.shape[0] #UpdateparametersbasedonLG for_inrange(0,depth): foriinrange(0,size): self.forward_pass(x_data[i1,]) self.backward_pass(x_data[i1,],y_data[i1]) theta_updated=self.encode() returntheta_updated
#Helperfunctions defsigmod(self,x): ''' Implementationofthesigmoidfunction ''' return1/(1+np.exp(-x)) defencode(self): ''' Encodethemodelparametersintoavector Output: -theta:vectorofparameters. ''' w1=self.11_weights.ravel() w2=self.12_weights.ravel() theta=np.concatenate([w1,w2,self.11_biases,self.12_biases]) returntheta defdecode(self,theta): ''' Decodethemodelparametersfromavector
```
11Input:
12 - theta:vectorofparameters.
13,,,
14w_layer1size=self.input_num*self.hidden_num
15w_layer2size=self.hidden_num*self.output_num
16
17w_layer1=theta[0:w_layer1size]
18self.11_weights=np.reshape(w_layer1,(self.input_num,self.hidden_num))
19
20w_layer2=theta[w_layer1size:w_layer1size+w_layer2size]
21self.12_weights=np.reshape(w_layer2,(self.hidden_num,self.output_num))
22self.11_biases=theta[w_layer1size+w_layer2size:w_layer1size+w_layer2size+self.hidden_num]
23self.12_biases=theta[w_layer1size+w_layer2size+self.hidden_num:w_layer1size+w_layer2size+self.hidden_num+self.output_num]
```

Listing 11: Python implementation of the Neural Network class

Finally, we implement an MCMC sampler following the Langevin-gradient MCMC sampling scheme. As discussed above, the prior and likelihood functions remain the same as in the linear regression model case, as presented in Listing 6.

```
1#MCMCsampler
2defsampler(self):
3,'
4RunthesamplerforadefinedNeuralNetworkmodel
5'''
6#define emptyarraystotorethesampledposteriorvalues
7#posteriorofallweightsandbiasoverallsamples
8pos_theta=np.ones((self.n_samples,self.theta_size))
9#posteriordefiningthevarianceofthenoiseinpredictions
10pos_tau=np.ones((self.n_samples,1))
11
12#recordoutputf(x)overallsamples
13pred_y=np.ones((self.n_samples,self.x_data.shape[0]))
14#recordsimulatedvaluesf(x)+erroroverallsamples
15sim_y=np.ones((self.n_samples,self.x_data.shape[0]))
16#recordtheRMSEofeachsample
17rmse_data=np.zeros(self.n_samples)
18
19#Initialisation
20#initialistheta
21theta=np.random.randn(self.theta_size)
22#makeinitialprediction
23pred_y[0,]=self.model.evaluate_proposal(self.x_data,theta)
24
25#initialiseeta
26eta=np.log(np.var(pred_y[0,]-self.y_data))
27tau_proposal=np.exp(eta)
28
29#Hyperpriors-consideredbylookingatdistributionofsimilartrainedmodels-i.e.distributionofweightsandbias
30sigma_squared=self.sigma_squared
31nu_1=self.nu_1
32nu_2=self.nu_2
33
34#calculatethepriorlikelihood
35prior_likelihood=self.prior_likelihood(sigma_squared,nu_1,nu_2,theta,tau_proposal)
36#calculatethelikelihoodconsideringobservations
37[likelihood,pred_y[0,],sim_y[0,],rmse_data[0]]=self.likelihood_function(theta,tau_proposal)
38
39n_accept=0
40n_langevin=0
#RuntheMCMCsampleforn_samples
* #Samplenewvaluesforthetaandtau

[MISSING_PAGE_POST]
```
10pos_tau[ii,]=tau_proposal
11else:
12#store
13pos_theta[ii,]=pos_theta[ii-1,]
14pos_tau[ii,]=pos_tau[ii-1,]
15
16#printthe%oftimestheproposalwasaccepted
17accept_ratio=(n_accept/self.n_samples)*100
18print(accept_ratio,'%wasaccepted')
19
20
21#storetheposteriorofthetaandtau,aswellastheRMSEoftheesamples
22self.pos_theta=pos_theta[self.n_burnin:,]
23self.pos_tau=pos_tau[self.n_burnin:,]
24self.rmse_data=rmse_data[self.n_burnin:]
25
26#Createapandasdataframetostoretheposteriorsamplesofthetaandtau,the
27associatedRMSE
28results_dict={'w()'.format(_):self.pos_theta[:,_].squeeze()for_inrange(self.theta_size-2)}
29results_dict['b0']=self.pos_theta[:,self.theta_size-2].squeeze()
30results_dict['bu']=self.pos_theta[:,self.theta_size-1].squeeze()
31results_dict['tau']=self.pos_tau.squeeze()
32results_dict['rmse']=self.rmse_data.squeeze()
33
34results_df=pd.DataFrameFrame.from_dict(
35results_dict
36)returnresults_df
```

Listing 12: Python implementation MCMC sampler - Langevin proposal distribution

### Results

We use the Sunspot time series5 data and Abalone6 datasets for regression problems. The Abalone dataset provides the ring age for Abalone based on eight features that represent physical properties such as length, width, and weight and associated target feature, i.e., the ring age. Determining the age of Abalone is difficult as requires cutting the shell and counting the number of rings using a microscope. However, other physical measurements can be used to predict the age and a model can be developed to use the physical features to determine the ring age. Sunspots are regions of reduced surface temperature in the Sun's photosphere caused by concentrations of magnetic field flux, and appear as spots darker than the surrounding areas. The Sunspot cycles are about every eleven years and over the solar cycle, the number of occurrence of Sunspot change more rapidly. Sunspot activities are monitored since they have an include on Earth's climate and weather systems. We obtain the Abalone dataset from the University of California (UCI) Machine Learning Repository7 and keep a processed version of all the datasets in our repository8. We also obtain datasets for classification problems from the same repository that features a number of datasets, but we only close Iris 9 and Ionosphere 10 classification data. The Iris classification dataset contains 4 features (sepal length, sepal width, petal length, petal with) of three types of Iris flower species, featuring 50 instances for each case. This dataset is one of the most prominent datasets used for machine learning. In the Ionosphere dataset, there are 34 continuous features and the target is either "good" or "bad", a binary classification task with 351 instances. The task is to filter the radio signals where the "good" radar refers certain structure in the ionosphere and "bad" imply that signals pass through the ionosphere.

Footnote 5: [https://www.sidc.be/silso/datafiles](https://www.sidc.be/silso/datafiles)

Footnote 6: [https://archive.ics.uci.edu/ml/datasets/abalone](https://archive.ics.uci.edu/ml/datasets/abalone)

Footnote 7: [https://archive-beta.ics.uci.edu/about](https://archive-beta.ics.uci.edu/about)

Footnote 8: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/tree/main/data](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/tree/main/data)

Footnote 9: [https://archive.ics.uci.edu/ml/datasets/iris](https://archive.ics.uci.edu/ml/datasets/iris)

Footnote 10: [https://archive.ics.uci.edu/ml/datasets/ionosphere](https://archive.ics.uci.edu/ml/datasets/ionosphere)

In the Sunspot time series problem, we employ a one-step ahead prediction and hence use one output neuron. We process the Sunspot dataset (univariate time series) using Taken's embedding theorem [118] to construct a state-space vector, i.e., the standard approach for using neural networks for time series prediction [59]. This is essentially
using a sliding window approach of size \(D\) overlapping \(T\) time lags. The window size \(D\) determines the number of input neurons in the Bayesian neural network and Bayesian linear model. We used \(D=4\) and \(T=2\) for our data reconstruction for the Sunspot time series, as these values have given good performance in our previous works [59].

We created a train and test set for all the datasets by selecting 2/3 of the dataset for training and remaining for testing. In the Bayesian linear model and neural network, we choose the number of samples to be 25,000 for all problems, distributed across 5 chains and excluding burn-in. In the Bayesian linear model, we choose the learning rate \(r=0.1\), and the step sizes for \(\theta=0.02\) and \(\tau=0.01\),respectively. Additionally, for the Gaussian prior distribution, we choose the parameters \(\sigma^{2}=5,\nu_{1}=0\) and \(\nu_{2}=0\), respectively. In the Bayesian neural network models, we choose the learning rate \(r=0.01\), and the step sizes for \(\theta=0.025\) and \(\tau=0.2\), respectively. In the Gaussian prior distribution, we choose the parameters \(\sigma^{2}=25,\nu_{1}=0\) and \(\nu_{2}=0\), respectively. We also use a burn in rate of 0.5 for both the Bayesian linear and the Bayesian neural network models. This indicates the proportion of the MCMC samples that are discarded. We report the mean and standard deviation (std) of the accuracy for respective problems (RMSE or classification performance), obtained from the posterior distribution after discarding the burn-in samples. We note that although only random walk proposal distribution can be used, we show results where both random walk and Langevin gradient proposals are used at a rate of 0.5.

We first present the results of Bayesian regression with the sunspot (time series) and Abalone (regression) datasets. We evaluate the model performance using the _root mean squared error_ (RMSE) which is a standard metric for time series prediction and regression problems. We present the results obtained by the Bayesian linear model and Bayesian neural network model for regression problems in Table 1. In Table 1, we observe that Bayesian neural network performs better for the Sunspot time series prediction problem as it achieves a lower RMSE on both the training and testing set. This can also be seen in Figures 1 and 1. In the case of Abalone problem, Table 1 shows that both models obtain similar classification performance, but Bayesian neural network has better test performance. However, we note that in both problems, Bayesian neural networks has a much lower acceptance rate; we prefer roughly a 23 % acceptance rate [119] that implies that the posterior distribution has been effectively sampled.

Table 2 presents results for the classification problems in the Iris and Ionosphere datasets. We notice that that both models have similar test and training classification performance for the Iris classification problem, and Bayesian neural networks give better results for the test dataset for the Ionosphere problem. The acceptance rate is much higher for Bayesian neural networks, the Bayesian linear model gets a suitable acceptance rate when compared to literature.

## 6 Convergence diagnosis

It is important to ensure that the MCMC sampling is adequately exploring the parameter space and constructing an accurate picture of the posterior distribution. One method of monitoring the performance of the adopted MCMC sampler is to examine convergence diagnostic to monitor the extent to which the Markov chains have become a stationary distribution. Practitioners routinely apply the Gelman-Rubin (GR) convergence diagnostic to this end [120]. This diagnostic is developed using sampling from multiple MCMC chains whereby the variance of each chain

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Method & Problem & Train (RMSE) & Test (RMSE) & Accept. rate \\ \hline \multirow{2}{*}{Bayesian linear model} & \multirow{2}{*}{Sunspot} & 0.025 & 0.022 & 13.5\% \\  & & (0.013) & (0.012) & \\ \multirow{2}{*}{Bayesian neural network} & \multirow{2}{*}{Sunspot} & 0.027 & 0.026 & 7.4\% \\  & & (0.007) & (0.007) & \\ \hline \multirow{2}{*}{Bayesian linear model} & \multirow{2}{*}{Abalone} & 0.085 & 0.086 & 5.8\% \\  & & (0.005) & (0.005) & \\ \multirow{2}{*}{Bayesian neural network} & \multirow{2}{*}{Abalone} & 0.080 & 0.080 & 3.8\% \\  & & (0.002) & (0.002) & \\ \hline \hline \end{tabular}
\end{table}
Table 1: Regression results using Bayesian linear model and Bayesian neural networks via MCMC. 100 samples are drawn from the parameter posteriors and the RMSE calculated with the mean and standard deviation (in brackets) reported.
\begin{table}
\begin{tabular}{c c c c c} \hline \hline Method & Problem & Train Accuracy & Test Accuracy & Accept. rate \\ \hline \multirow{2}{*}{Bayesian linear model} & \multirow{2}{*}{Iris} & 90.392\% & 90.844\% & 83.5\% \\  & & (2.832) & (3.039) & \\ \multirow{2}{*}{Bayesian neural network} & \multirow{2}{*}{Iris} & 97.377\% & 98.116\% & 97.0\% \\  & & (0.655) & (1.657) & \\ \hline \multirow{2}{*}{Bayesian linear model} & \multirow{2}{*}{Ionosphere} & 89.060\% & 85.316\% & 58.8\% \\  & & (1.335) & (2.390) & \\ \multirow{2}{*}{Bayesian neural network} & \multirow{2}{*}{Ionosphere} & 99.632\% & 92.668\% & 94.5\% \\  & & (0.356) & (1.890) & \\ \hline \hline \end{tabular}
\end{table}
Table 2: Classification accuracy with Bayesian linear model and Bayesian neural networks via MCMC. We show results for the posterior mean and standard deviation (in brackets).

Figure 10: Bayesian linear regression model - Sunspot dataset
Figure 11: Bayesian neural network model - Sunspot dataset

Figure 12: Distribution of \(\hat{R}\) values for our Bayesian linear model and Bayesian neural network.
is assessed independently (within-chain variance) and then compared to the variance between the multiple chains (between-chain variance) for each parameter. Large differences between these two variances would indicate the at the chains have not converged on the same stationary distribution.

For our examples above, we will therefore run a number of independent experiments and compare the MCMC chains using the convergence diagnostic. This measure of convergence is calculated as the potential scale reduction factor (PSRF or \(\hat{R}\)), where values of PSRF close to 1 indicate convergence. The PSRF is obtained by first calculating the between-chain variance

\[B=\frac{n}{m-1}\sum_{j=1}^{m}(\bar{\theta}_{j}-\bar{\bar{\theta}})^{2} \tag{41}\]

and within-chain variance

\[W=\frac{1}{m}\sum_{j=1}^{m}s_{j}^{2} \tag{42}\]

where \(n\) is the number of samples, \(m\) is the number of chains,

\[\bar{\theta}_{j}=\frac{1}{n}\sum_{i=1}^{n}\theta_{ij}\]

\[\bar{\bar{\theta}}=\frac{1}{m}\sum_{j=1}^{m}\bar{\theta}_{j}\]

and

\[s_{j}^{2}=\frac{1}{n-1}\sum_{j=1}^{n}\theta_{ij}-\bar{\theta}_{j}\]

After observing these estimates, we may then estimate the target variance \(\sigma^{2}\) with

\[\hat{\sigma}^{2}=\frac{n-1}{n}W+\frac{B}{n}\]

Then what is known about \(\theta\) can be estimated and the result is an approximate Student t's distribution for \(\theta\) with centre \(\bar{\bar{\theta}}\), scale \(\sqrt{\hat{V}}=\sqrt{\hat{\sigma}^{2}+\frac{B}{mn}}\) and degrees of freedom \(d=\frac{2\hat{V}}{\sqrt{\hat{\sigma}^{2}}(\hat{V})}\). Here

\[\dot{\mathrm{var}}(\hat{V}) =\left(\frac{n-1}{n}\right)^{2}\frac{1}{m}\dot{\mathrm{var}}(s_ {i}^{2})+\left(\frac{m+1}{mn}\right)^{2}\frac{2}{m-1}B^{2}\] \[+2\frac{(m+1)(n-1)}{mn^{2}}\] \[\cdot\frac{n}{m}[\mathrm{c}\dot{\mathrm{ov}}(s_{i}^{2},\bar{ \theta}_{i}^{2})-2\bar{\bar{\theta}}\mathrm{c}\dot{\mathrm{ov}}(s_{i}^{2},\bar {\theta}_{i})].\]

Finally, we calculate the PSRF given by \(\hat{R}\) as

\[\sqrt{\hat{R}}=\sqrt{\frac{\hat{V}}{W}\frac{d}{d-2}} \tag{43}\]

Listing 13 gives further details of the Python implementation 11.

Footnote 11: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/tree/main/convergence/convergence-GR.py](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial/tree/main/convergence/convergence-GR.py)
### Results

We show results for Gelman-Rubin diagnosis for Bayesian linear regression and Bayesian neural networks using the the Sunspot problem. We test five chains for each model, each with 5,000 samples of weights excluding burn-in samples and evaluating the \(\hat{R}\) values for each weight. Figure 12 shows the distribution of the \(\hat{R}\) values on a log-scale, and we observe that the \(\hat{R}\) values of the weights for Bayesian linear regression model are all much smaller than Bayesian neural network. We can conclude that the Bayesian neural network shows poor convergence based on the Gelman-Rubin diagnostics presented.

By closely examining each weight individually, we observe that the problem of non-convergence mainly arises from multi-modality of the posterior. In Figure 13 we look at samples from a single chain of 25,000 samples excluding 50% burn-in. In Figure 13a and 13b, we present a visualisation for a selected weight from the Bayesian linear model. In Figure 13c and 13d, we present a visualisation for a selected weight from the Bayesian neural network model. We observe potential multi-modal distributions in both cases, with a high degree of auto-correlation and poor convergence. To examine the impact of longer MCMC chains in achieving convergence, an additional test was run taking 400,000 samples excluding 20% burn-in and then thinning the chain by a factor of 50. These thinned results are shown for the Iris dataset in Figure 14. We can see that the chains exhibit more desirable properties and (particularly in the case of the linear model as shown in Figure 14a,b) convergence diagnostics were improved as expected when removing some of the auto-correlation through thinning. We can see in Figure 14c,d that the Bayesian neural network still exhibits a multi-model posterior for this parameter.

## 7 Discussion

Overall, this highlights the problem of applying Gelman-Rubin diagnostics on the weights of Bayesian neural networks. We observe that the Bayesian neural network performs better than the Bayesian linear regression in most
Figure 13: Posterior and trace-plot for a parameter in each of the Bayesian linear and Bayesian Neural Net models - Sunspot data. For this example, 25,000 samples were taken excluding 20% burn-in with the chain thinned by a factor of 50.
Figure 14: Posterior and trace-plot for a paramater in each of the Bayesian linear and Bayesian Neural Net models - Iris data. For this example, 400,000 samples were taken excluding 50% burn-in with no thinning.
prediction problems, despite showing no indications of convergence. This is due to the challenge of attaining MCMC convergence in the weights of Bayesian neural networks, as they are influenced by the problem of multi-modality. Hence, we conclude that for a Bayesian neural network, a poor performance in the Gelman-Rubin diagnostics does not necessarily imply a poor performance in prediction tasks. We revisit the principle of equifinality [121, 122] that states that in open systems, a given end state can be reached by many potential means. In our case, the system is a neural network model and many solutions exist that represent a trained model displaying some level of performance accuracy. Despite this, we note that the goal of Bayesian models are to offer uncertainty quantification in predictions, and proper convergence is required. The original Gelman-Rubin diagnosis [120] motivated a number of enhancements for different types of problems [123, 124, 125]; and we may need to develop better diagnosis for Bayesian neural networks. Nonetheless, in our case comparing Bayesian logistic regression (converged) and Bayesian neural networks (not converged), but achieved good accuracy, we can safely state that so our our Bayesian neural network can only provide a means for uncertainty quantification but is not mature enough to qualify as a robust Bayesian model.

We also revisit convergence issue in the case of Bayesian linear model as shown in Panel (a) - Figure 13, which shows a multi-modal distribution that has not well converged. In the case of linear models, we can have certain features that are not contributing much to the decision making process (predictions) and the weights associated (coefficients) with those features may not converge since they are not really contributing. This is similar to the case of neural networks, where certain weight links are not needed and can be pruned.

We note that we attach much higher acceptance rate in Bayesian neural networks for classification problems (Table 2) when compared to regression problems(Table 2). We note that multinomial likelihood is used for the classification, and Gaussian likelihood is used for the regression problems, and both use Langevin-gradient proposal distributions that account for detailed balanced condition. Hence, we need to fine tune the hyper-parameters associated with the proposal distribution to ensure we get higher acceptance rate for the regression problems. Although 23 % acceptance rate [119] has been prominently used as a "golden rule", the optimal acceptance rate depends on the nature of the problem. The number of parameters, Langevin-based proposal distribution, and type of model would raise questions about the established acceptance rate [126]. Hence, more work needs to be done to establish what acceptance rates are appropriate for simple neural networks and deep learning models.

A way to address the issue of convergence would be to develop an ensemble of linear models that can complete with the accuracy of neural networks or deep learning models. In ensemble methods such as bagging and boosting, we can use linear models that have attained convergence as per Gelman-Rubin diagnosis and then combine results of the ensemble using averaging and voting, as done in ensemble methods. Furthermore, we note that Langevin-based MCMC is a way ahead but not the only way ahead. A number of works have utilised the Hamiltonian MCMC which provides another approach to incorporate gradients in MCMC methods that can provide advantages over Langevin-based methods; however, we need comprehensive evaluation studies to check how both methods work for deep learning models such as CNN and LSTM.

Our previous work has shown that despite the challenges, the combination of Langevin-gradients with parallel tempering MCMC [59], presents opportunities for sampling larger neural network architectures [66, 65]. The need to feature a robust methodology for uncertainty quantification in CNNs will make them more suitable for applications where uncertainty in decision making poses major risks, such as medical image analysis [127] and human security [128]. Recently, CNNs have been considered for modelling temporal sequences, they have proven to be successful for time series classification [129, 130], and time series forecasting problems [131, 132]. Our recent work has shown that one dimensional CNNs provide better prediction performance than LSTM-based models for multi-step time series prediction problems [133]. Leveraging CNNs within a Bayesian framework can provide better uncertainty quantification in predictions and make them useful for cutting-edge real-world applications.

We envision that this tutorial will enable statisticians and machine learners to utilise MCMC sampling in a more effective way when developing new models and also developing a Bayesian framework for existing deep learning models. The tutorial has introduced basic concepts with code and provides an overview of challenges, when it comes to convergence of Bayesian neural networks.
## 8 Code and Data

All code (for implementation, results and figures) and data presented in this paper are available in the associated github repository12. This repository presents the implementations in separate Jupyter notebooks in the base directory, with sub-directories containing data, convenient functions and details of the environment setup.

Footnote 12: [https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial](https://github.com/sydney-machine-learning/Bayesianneuralnetworks-MCMC-tutorial)
