Bayesian Methods and Monte Carlo Simulations

Abstract

Bayesian methods provide the means for studying probabilistic models of linear as well as non-linear stochastic systems. They allow tracking changes in probability distributions by applying Bayes’s theorem and the chain rule for factoring the probabilities. However, an excessive complexity of resulting distributions often dictates the use of numerical methods when performing statistical and causal inferences over probabilistic models. In this chapter, the Bayesian methods for intractable distributions are first introduced as sampling, filtering, approximation, and likelihood-free methods. Their fundamental principles are explained, and the key challenges are identified. The concise survey of Bayesian methods is followed by outlining their applications. In particular, Bayesian experiment design aims at maximizing information gain or utility, and it is often combined with an optimum model selection. Bayesian hypothesis testing introduces optimality in the data-driven decision making. Bayesian machine learning assumes data labels to be random variables. Bayesian optimization is a powerful strategy for configuring and optimizing large-scale complex systems, for which conventional optimization techniques are usually ineffective. The chapter is concluded by examining Bayesian Monte Carlo simulations. It is proposed that augmented Monte Carlo simulations can achieve explainability and also provide much better information efficiency.

1. Introduction

Many real-world systems exhibit some level or form of randomness. This is reflected in their mathematical models, which are often probabilistic. There are two basic strategies to interpret the randomness captured by probabilistic models. The frequency of occurrence of a random phenomenon is a measurable quantity that can be used to predict how likely the phenomenon can be observed in future. The other approach quantifies the uncertainty about the occurrence of random phenomenon more subjectively, that is, as a degree of belief or expectation of observing the phenomenon, given domain knowledge and the past experiences. This latter approach led to a broad area of general Bayesian methods [1], Bayesian data analyses [2], Bayesian signal processing [3], Bayesian regression [4], Bayesian machine learning [5], and Bayesian optimization [6]. In the tasks of statistical inference, the assumption of knowing (or not) a prior distribution crucially affects the feasibility as well as the structure of estimators [7]. There is also an intimate connection between Bayesian probabilistic models and making causal inferences [8], as will be discussed in Section 3.

Bayesian methods found widespread applications in many probabilistic modeling frameworks. These methods are all rooted in the surprisingly simple Bayes’s theorem, that is,

$$p\left(\theta |x\right)=\frac{p\left(x|\theta \right)p\left(\theta \right)}{p\left(x\right)}$$E1

where $p\left(\cdot \right)$ denotes the (for continuous random variables) or the probability (for discrete random variables). Eq. (1) quantifies how our belief about the parameter $\theta$ changes from its prior, $p\left(\theta \right)$, to posterior, $p\left(\theta |x\right)$, after observing data $x$. The conditional term, $p\left(x|\theta \right)$, represents the likelihood of the parameter $\theta$, given observations, $x$. The scaling term, $p\left(x\right)$, in Eq. (1) is usually referred to as evidence. It should be noted that both the parameter and the data can be assumed in an arbitrary number of dimensions.

The majority of tasks in Bayesian inference involves explicitly or implicitly solving one of the following integrals (or sums, for discrete random variables), that is,

$$\text{marginalization}:p\left(x\right)=\int p\left(x,\theta \right)\mathit{d\theta }$$E2

$$\text{summarization}: E\left[f\left(x\right)|\theta \right]=\int f\left(x\right)p\left(x|\theta \right)\mathit{dx}$$E3

$$\text{prediction}: p\left(x_{t+1}\right)=\int p\left(x_{t+1}|x_t\right)p\left(x_t\right){\mathit{dx}}_t$$E4

$$p\left(x_{t+1}|x_t\right)=\int p\left(x_{t+1}|\theta \right)p\left(\theta |x_t\right)\mathit{d\theta }.$$

Unfortunately, in most but a few real-world scenarios, the expressions (2)–(4) are not mathematically tractable. In particular, the distributions often involve multiple sums and/or integrals, and their closed-form expressions cannot be obtained. The distributions are sometimes only known up to a scaling constant. Then, even numerically computing Eq. (1) can be rather challenging, since many complex distributions are multi-modal with a large number of local minima and maxima. Moreover, in online data processing, the distributions must be updated continuously as soon as the new data arrives.

The rest of this chapter is organized as follows. The strategies for performing Bayesian inferences with intractable distributions are outlined in Section 2 including sampling, filtering, approximation, and likelihood-free methods. The applications of Bayesian inferences are discussed in Section 3 including Bayesian experiment design, Bayesian hypothesis testing, Bayesian machine learning, and Bayesian optimization. Bayesian Monte Carlo simulations are reviewed in Section 4. Although the chapter mostly reviews known concepts and frameworks in Bayesian analysis of probabilistic models, Section 4 contributes a description of augmented Monte Carlo simulations. These simulations aim at providing explainability and improving information gain.

The references cited at the end of this chapter are by no means comprehensive; rather, they are suggestions of initial readings where to find further information about the topics discussed in this chapter.

2. Bayesian inference with intractable distributions

Bayesian inference assumes that prior knowledge is available, so it should be exploited to improve the estimation accuracy. Such knowledge is often represented as a prior distribution of the parameters to be estimated. It can be provided by the domain expertise or otherwise be subjectively selected. As already mentioned, non-linear and high-dimensional models are rarely tractable analytically. This section presents four basic strategies on how to deal with intractable distributions in Bayesian inference. In particular, the sampling methods define a proposal distribution, which is much easier to sample from. Most of these methods form a Markov chain of samples, which gradually converges to the desired distribution. The filtering methods update the estimates iteratively, so they are particularly suitable for streaming data. These methods are also more efficient in higher dimensions than the pure sampling methods. The approximate methods strive to improve the efficiency, even at the cost of estimation accuracy. The last group is the likelihood-free methods, which are particularly suited for parameter estimation using simulations.

2.1 Sampling methods

Consider, for example, the expectation (3). The integral (3) may be not only mathematically intractable but also numerically intractable, provided that the distribution $p\left(x\right)$ is difficult or even impossible to sample from. A common strategy then is to choose and sample from another, simpler distribution, $q\left(x\right)$, such that $q\left(x\right)>0$, whenever $f\left(x\right)p\left(x\right)>0$. The resulting method is known as importance sampling (IS), and $q\left(x\right)$ is referred to as the proposal distribution. Eq. (3) is then rewritten as

$$\int f\left(x\right)p\left(x\right)\mathit{dx}=\int f\left(x\right)\frac{p\left(x\right)}{q\left(x\right)}q\left(x\right)\mathit{dx}\approx \frac{1}{N}\sum _{i=1}^{N}f\left(x^{\left(i\right)}\right)\underset{w^{\left(i\right)}}{\underbrace{\frac{p\left(x^{\left(i\right)}\right)}{q\left(x^{\left(i\right)}\right)}}}$$E5

where the sum is an unbiased and consistent estimator of the expectation due to the law of large numbers.

Although the estimate (5) is unbiased and consistent, the number of samples required grows exponentially with the number of dimensions. It is also desirable to make the weights, $w^{\left(i\right)}$, more uniform; that is, $q\left(x\right)$ should approximate $p\left(x\right)$, in order to reduce the variance of the estimator in (5). Alternatively, the distribution $p\left(x\right)$ could be approximated by a histogram; the bins of such histogram can be optimized for a maximum efficiency. The histogram can be sampled directly using a uniform random number generator.

The sampling efficiency of the IS methods, particularly in high dimensions, can be improved by accounting for the likelihood of samples. The Markov chain Monte Carlo (MCMC) strategy represents a broad variety of sampling methods, such that the current sample is affected by the location of the previous sample. It allows learning complex target distributions in multiple dimensions. However, the MCMC sampling requires a burn-in period for achieving the convergence. Furthermore, the samples may get stuck in areas of small probability (requiring a random restart), and the correlation between consecutive samples increases the estimator variance.

Metropolis–Hastings sampling is the most commonly used MCMC algorithm. It can be succinctly described as the following three steps:

1. Given the current state $x_t$, sample a new state, $x_{t+1}$, from the proposal $q\left(x_{t+1}|x_t\right)$;

2. Calculate the acceptance probability $\alpha \left(x_{t+1},x_t\right)$ of the new state $x_{t+1}$;

3. If $q\left(x_t|x_{t+1}\right)=q\left(x_{t+1}|x_t\right)$ is symmetric, then accept the new state $x_{t+1}$ if $\pi \left(x_{t+1}\right)>\pi \left(x_t\right)$; otherwise remain at the current state, that is, $x_{t+1}\equiv x_t$.

The challenge is to find a good proposal, $q$, with a bounded support to achieve fast convergence to the target distribution, $\pi \left(x\right)$. The target distribution only needs to be known up to a proportionality constant. The generated samples tend to concentrate in the areas of large probability. The generated samples occasionally contain sub-sequences of the same value. More importantly, as with any other sampling methods, there is an exploration–exploitation trade-off. Large acceptance rate means slow exploration of the target distribution, $\pi$, as well as large correlations between the samples. On the other hand, a small acceptance rate means large jumps across the support of $\pi$. The acceptance rate can be controlled by assuming a parameterized proposal, $q_{\varphi }$. The acceptance rate of about $50\%$ is recommended for samples in a few dimensions, and it should be reduced to about $25\%$ for distributions in many dimensions.

In the literature, many variants of the Metropolis–Hastings sampling algorithm were proposed, for example, combining random walk sampling, adaptive rejection sampling, Langevin algorithm, and augmented estimator.

Gibbs sampling iteratively samples from the conditional densities

$$\begin{array}{ccc}{X}_1^{\left(t\right)}& \sim & p_1\left(x_1|x_2^{\left(t\right)},x_3^{\left(t\right)},\dots ,x_N^{\left(t\right)}\right)\\ X_2^{\left(t\right)}& \sim & p_2\left(x_2|x_1^{\left(t\right)},x_3^{\left(t\right)},\dots ,x_N^{\left(t\right)}\right)\\ & ⋮& \\ X_N^{\left(t\right)}& \sim & p_N\left(x_N|x_1^{\left(t\right)},x_2^{\left(t\right)},\dots ,x_{N-1}^{\left(t\right)}\right).\end{array}$$E6

It can be represented as a set of $N$ univariate Metropolis–Hastings samplers. The sampling dimensions can be chosen systematically or at random. Even though Gibbs sampling has $100\%$ acceptance rate, it can still experience a slow convergence to the target distribution, so it is often combined with other sampling methods.

Gibbs sampling motivates the Rao–Blackwell estimators

$$E\left[h\left(x_1\right)\right]\approx \frac{1}{T}\sum _{t=1}^{T}E\left[h\left(x_1\right)|x_2^{\left(t\right)},\dots ,x_N^{\left(t\right)}\right]$$E7

$$p_i\left(x_i\right)\approx \frac{1}{T}\sum _{t=1}^T{p}_i\left(x_i|x_j^{\left(t\right)},j\ne i\right)$$E8

of function expectation and of the marginal distribution, respectively. These estimators are unbiased and have lower variance, since

$$\mathit{var}\left[E\left[h\left(x_1\right)|x_2^{\left(t\right)},\dots ,x_N^{\left(t\right)}\right]\right]\le \mathit{var}\left[E\left[h\left(x_1\right)\right]\right].$$E9

An alternative strategy combines importance sampling with MCMC sampling assuming independent proposals, ${\prod }_{i=1}^N{q}_i\left(x_i\right)$ or ${\prod }_{i=1}^N{q}_t\left(x_i|x_i^{\left(t-1\right)}\right)$.

An interesting problem is how to design MCMC sampling when the target distribution is not stationary. The adaptive MCMC sampling can be trained online with the generation of new data. However, using past samples for adaptation invalidates the Markov assumption, so the convergence to the target distribution may be problematic, or the adaptation should cease after the burn-in period.

The Hamiltonian MCMC sampling tracks movements of a hypothetical ball under the potential and kinetic energy constraints. The new sample representing the updated ball location is obtained by integrating the ball’s speed. The integration can be approximated, for example, by discretization. More importantly, the interval of integration trades-off the acceptance rate and the sampling rate, so, also the amount of correlations between subsequent samples and the waiting time. The main advantage of Hamiltonian MCMC method is a fast mix-in; that is, the samples converge quickly to the target distribution.

Another sampling method motivated by models in statistical physics is a restricted Boltzmann machine (RBM). This method learns the target distribution as a configuration, $h$, of a bipartite graph. The target distribution is $p\left(h\right)=exp\left(-E\left(h\right)\right)/Z$, where $E\left(h\right)$ is the energy of the configuration $h$, and $Z={\sum }_{h}exp\left(-E\left(h\right)\right)$ is a partition (scaling) function.

The reversible-jump MCMC sampling can be used for distributions having an uncertain number of parameters (the model order). The idea is to increase the number of parameters by 1 with a certain probability every time the sample is taken from the proposal distribution.

The common pitfall of all sampling methods is how to recognize that the generated samples are not converging to the target distribution, for example, since the assumptions were violated, or the method has not been implemented correctly. These issues may not be easily detected by simply evaluating the samples or the simulation outputs.

2.2 Filtering methods

Processing time-series and streaming data must be often done recursively. The corresponding data model can be usually represented as a dynamic system, so the current state, $x_t$, is updated from the previous state, $x_{t-1}$, by an innovation random process, and the states, $x_t$, are observed as values, $z_t$. The states form a Markov chain, that is, $p\left(x_t|x_{t-1},\dots ,x_1\right)=p\left(x_t|x_{t-1}\right)$, and the observations are assumed to be independent, that is, $p\left(z_i,z_j|x_t\right)=p\left(z_i|x_t\right)p\left(z_j|x_t\right)$. Assuming the Bayes’s theorem

$$p\left(x_{1:t}|z_{1:t}\right)=\frac{p\left(z_{1:t}|x_{1:t}\right)p\left(x_{1:t}|z_{1:t-1}\right)}{p\left(z_t|z_{1:t-1}\right)},$$E10

the states can be estimated recursively from the incoming observations as a two-step procedure. The first step predicts the (distribution of) current state as

$$p\left(x_t|z_{1:t-1}\right)=\int p\left(x_t|x_{t-1}\right)p\left(x_{t-1}|z_{1:t-1}\right){\mathit{dx}}_{t-1}.$$E11

The predicted state is corrected in the second step after receiving the latest observation, $z_t$, as

$$p\left(x_t|z_{1:t}\right)=\frac{p\left(z_t|x_t\right)p\left(x_t|z_{1:t-1}\right)}{\int p\left(z_t|x_t\right)p\left(x_t|z_{1:t-1}\right){\mathit{dx}}_t}.$$E12

For linear systems with Gaussian inputs, the processes remain Gaussian, so only their mean values, $E\left[x_{1:t}|z_{1:t}\right]$, and covariances, $\mathit{cov}\left[x_{1:t}|z_{1:t}\right]$, need to be tracked. In such a case, the recursive estimator implementing (11) and (12) is known as Kalman filter. The Kalman filtering can be succinctly described as

$${\widehat{x}}_t=\left(\text{prediction}\text{of}{x}_t\right)+K_t\cdot \left[z_t-\left(\text{prediction}\text{of}{z}_t\right)\right]$$E13

where $K_t$ is the Kalman gain at time $t$. It can be shown that Kalman filter is unbiased, and it is optimum in a sense of having the smallest possible variance; it belongs to a class of the best linear unbiased estimators (BLUE). Practical implementations of Kalman filter are usually concerned with numerical stability and process observability. The posterior distribution of the model states can be estimated using kernel and other smoothing methods.

Kalman filter can be modified to work with non-linear models. In particular, extended Kalman filter performs local linearization at each iteration assuming the first-order Taylor expansion. This can, however, cause large estimation errors and stability issues. Unscented Kalman filter defines a set of sigma points, which are tracked through non-linearity. The weighted and transformed sigma points are used to reconstruct the mean and the covariance of Gaussian distribution, so canonical Kalman filter can be then used. Ensemble Kalman filter has been adopted for large but sparse systems. Expectation maximization (EM) filtering is used when there are also unknown model parameters; see Eqs. (16) and (17) below.

In practice, having an approximate solution for a complex model is often much more useful than obtaining the exact solution of a simplified model. This has motivated the development of particle filters to estimate random signals under non-linear and non-Gaussian conditions. The key idea is to represent the posterior distribution by a set of random weighted samples (called particles), $x_i$, that is,

$$p\left(x\right)\approx \sum _{i=1}^N{w}_i\delta \left(x-{\tilde{x}}_i\right)\Rightarrow E\left[f\left(x\right)\right]\approx \sum _{i=1}^N{w}_{i}f\left({\tilde{x}}_i\right)$$E14

where $\delta$ denotes the Dirac-delta function. However, the number of required particles grows quickly with the number of dimensions.

The particles should be chosen to represent high-density regions by many samples with large weights. Provided that the particles are generated using some MCMC sampling strategies, the resulting recursive Bayesian estimator is referred to as the sequential MC (SMC) method. Since only after a few iterations, most particles often have negligible weights (so-called particle degeneracy problem), particle resampling (with replacement and proportionally to their weight) is necessary in order to recover a good representation of the underlying distribution while maintaining a constant number of particles.

Expectation maximization (EM) method has been developed to deal with missing data and to filter data under parameter uncertainty. Assuming the observed but incomplete data, $x$, as well as the missing data, $y$, the data likelihood conditioned on an unknown parameter, $\theta$, can be computed as

$$p\left(x|\theta \right)=\int p\left(x,y|\theta \right)\mathit{dy}=\int p\left(y|x,\theta \right)p\left(x|\theta \right)\mathit{dy}.$$E15

The integral (15) is evaluated iteratively by first computing the Q-function (the E-step) as

$$Q\left(\theta ,{\theta }_{n-1}\right)=E\left[logp\left(x,y|\theta \right)|x,{\theta }_{n-1}\right]$$E16

followed by estimating the parameter $\theta$ (the M-step) as

$${\theta }_n={\text{argmax}}_{\theta }Q\left(\theta ,{\theta }_{n-1}\right).$$E17

This procedure guarantees that the likelihood is maximized, since always $p\left(x|{\theta }_n\right)>p\left(x|{\theta }_{n-1}\right)$. However, Eqs. (16) and (17) are normally analytically intractable, so they must be evaluated numerically.

2.3 Approximation methods

The sampling methods, which mostly involve MCMC sampling, for performing Bayesian inference are asymptotically exact. However, these methods are also numerically expensive, especially for high-dimensional models. When working with big data, or when there is a need to choose among many plausible models, having less accurate but numerically cheaper solution is highly desirable. Variational inference deterministically approximates the log-evidence $p\left(x\right)$ in Bayes’s theorem (1) as

$$logp\left(x\right)=\mathcal{L}+\mathit{KL}\left[q\left(\theta \right)\Vert p(\theta |x)\right]\ge \mathcal{L}$$E18

where $\mathcal{L}$ is the evidence lower-bound (ELBO), since the Kullback–Leibler (KL) divergence is always non-negative. Therefore, given the evidence $p\left(x\right)$, minimizing the KL divergence between the approximation, $q\left(\theta \right)$, and the posterior, $p\left(\theta |x\right)$, is equivalent to maximizing ELBO, $\mathcal{L}$. Alternatively, it is possible to approximate the prior, $p\left(\theta \right)$, instead of the posterior. These are functional optimization problems defined as $\frac{d}{\mathit{dq}}\mathcal{L}\left(q\right)=0$, s.t., $\int q\left(\theta \right)\mathit{d\theta }=1$, with the optimum solutions $q^{\ast }\left(\theta \right)=p\left(\theta |x\right)$ or $q^{\ast }\left(\theta \right)=p\left(\theta \right)$.

Any function, $q\left(\theta \right)$, can be assumed, but some choices are better than others. The mean-field approximation assumes that the dimensions are independent, that is, $q\left(\theta \right)={\prod }_{i=1}^M{q}_i\left({\theta }_i\right)$. Then, $q_i\left({\theta }_i\right)$ can be optimized in turn until a convergence is obtained, indicated by no more increases in the ELBO value. There is a performance penalty if the independence assumption is not satisfied. Another popular choice is to assume a class of distributions, $q_{\varphi }\left(\theta \right)$, parameterized by $\varphi$.

The mean-field assumption is implemented in the coordinated ascent variational inference (CAVI) algorithm. It is related to Gibbs sampling and message passing. However, this method is not computationally efficient and only works for distributions, which are conditional conjugates of the prior. Stochastic variational inference (SVI) improves the efficiency as well as accuracy of CAVI by assuming the stochastic approximation of the ELBO gradient and by optimizing the parameters using only a subset of data in each iteration (similarly to the EM method). Unfortunately, both CAVI and SVI require deriving the approximation components, $q_i$, analytically, which is often impractical.

Motivated by the limitations of conventional variational inference methods, the black box variational inference (BBVI) algorithm samples from a family of distributions, $q_{\varphi }\left(\theta \right)$, assuming that $p\left(x|\theta \right)p\left(\theta \right)=p\left(x,\theta \right)$ is known. The objective is to avoid the need for analytical derivations. The BBVI method can be used to estimate a wide range of linear and non-linear models. Since it is also the most efficient method for performing variational inference, it has become available as a library in many programming languages.

The divergence between prior or posterior and its approximation cannot be performed directly but only by maximizing the ELBO. Although, at each iteration, the ELBO is guaranteed not to decrease, the convergence can be slow. However, the efficiency of variational methods is still better than that of the sampling and EM methods. The KL divergence can be replaced with alternative measures such as $\alpha$ divergence, $f$ divergence, and mutual information. In addition, there are only a few theoretical guarantees, and the overall estimation accuracy of variational methods is crucially affected by not knowing how well the posterior (or prior) has been approximated and approximating asymmetric distributions is more challenging.

2.4 Likelihood-free inference

In many practical scenarios, a mathematical model can be available as a simulation. It is often easy to simulate the model for different parameter values. Thus, by generating and comparing many simulation outputs for different parameter values, the simulation run can be identified that best matches the observed data. The parameter values for this simulation run are then declared to be the estimate. The inference methods implementing this simple idea are known as being likelihood free. They can be found in the literature as indirect inference, bootstrap filter, synthetic likelihood, and other methods. However, the most popular among these methods is approximate Bayesian computation (ABC).

The ABC method only requires that a generative model to generate data is available. There is otherwise no need to assume, know, or calculate the model posterior, parameter likelihoods, or importance weights. The basic ABC rejection sampling is performed by repeating the following three steps:

1. Sample the parameter $\theta$ from its prior, $p_0\left(\theta \right)$;

2. Simulate or otherwise generate the data $x\left(\theta \right)$;

3. If the distance, $‖x-\tilde{x}‖<\epsilon$, between the simulated data, $x$, and the observed data, $\tilde{x}$, is sufficiently small, record the parameter $\theta$.

The posterior distribution, ${\pi }_{\epsilon }\left(\theta |x\right)$, can be estimated after enough samples, $\theta$ and $x$, have been collected. It is clear that the number of samples required to accurately estimate ${\pi }_{\epsilon }\left(\theta |x\right)$ grows quickly with the number of dimensions. The exact inference can be obtained in the limit ${lim}_{\epsilon \to 0}{\pi }_{\epsilon }\left(\theta \right)=p\left(\theta |x\right)$. On the other hand, no learning occurs if ${lim}_{\epsilon \to \infty }{\pi }_{\epsilon }\left(\theta \right)=p_0\left(\theta \right)$.

The efficiency of the ABC method can be improved by considering summary statistics, $S\left(x\right)$, rather than directly comparing the data, $x$. If $S\left(x\right)$ is the sufficient statistic for estimating the parameter $\theta$, then ${\pi }_{\epsilon }\left(\theta |S\left(x\right)\right)={\pi }_{\epsilon }\left(\theta |x\right)$; however, in practice, an informative statistic is often used.

Even though the ABC method can be readily parallelized, it can still be very inefficient, especially if the informative statistic has been poorly chosen and in high dimensions. Since the default ABC method does not exploit information about the accepted and rejected values of $\theta$, the efficiency could be improved by employing the MCMC sampling. The resulting ABC-MCMC algorithm performs the following four steps:

1. Augment the posterior distribution as $p\left(\tilde{\theta },\tilde{x}|x\right)$;

2. At the current state $\tilde{\theta }$, generate a new sample, $\theta \sim q\left(\theta |\tilde{\theta }\right)$;

3. Simulate the data, $x$, for the new state $\theta$;

4. Accept the new state, $\theta$, with some probability, $\alpha \left(\theta ,x\right)$.

This algorithm can be fine-tuned by gradually reducing $\epsilon$ in order to balance rejections with the estimation accuracy while the convergence requires $\epsilon$ to be sufficiently small. Unlike pure MCMC sampling, the MCMC-ABC method requires a small acceptance rate in order to achieve good accuracy.

3. Other topics in Bayesian analysis

The Bayesian frameworks find many real-world applications in data analysis and data-driven planning and decision making. Bayesian experiment design is about choosing experiments and models that can provide the largest information gain or utility. Bayesian hypothesis testing accounts for prior knowledge to bias the hypothesis rejection or acceptance. This framework can be also used to make statistically informed decisions. Bayesian machine learning has become increasingly popular in recent years. It considers data labels to be random variables. Both supervised and semi-supervised strategies are discussed. Bayesian optimization shows a remarkable promise in optimizing complex, difficult-to-evaluate systems. It considers the observed system responses to be samples of a random process.

3.1 Bayesian experiment design

In Bayesian experiment design, the task is to choose the optimum experiment and possibly also a data model [9]. This can be probabilistically represented as the augmented posterior (1), that is,

$$p\left(\theta |x,d,m\right)=\frac{p\left(x|\theta ,d,m\right)p\left(\theta ,m\right)}{p\left(x|d,m\right)}$$E19

where $d$ and $m$, respectively, represent the experiment and the model choice. More specifically, $\theta$ denotes uncontrolled inputs to the experiment as well as unknown model parameters, whereas $d$ are the experiment inputs that can be controlled, that is, designed by a selection. The objective is to specify the optimum design, $d$, in order to facilitate the estimation of $\theta$ as well as selection of the best model, $m$, for observations, $x$.

For every configuration of the experiment including the model selection, let $U\left(d,x,m,\theta \right)$ be the perceived utility, for instance, the amount of information that can be gained from the experiment. Since the observations are random and the data model is unknown, the optimum experiment maximizes the average utility,

$$d^{\ast }={\text{argmax}}_{d\in \mathcal{D}}\overline{U}\left(d\right)={\text{argmax}}_{d\in \mathcal{D}}{E}_{x,m,\theta }\left[U\left(d,x,m,\theta \right)\right].$$E20

Alternatively, the model selection can be formulated as a hypothesis-testing problem. The observed data and models are combined to obtain and then sample the predictive distribution for each model candidate. Additional experiments can be performed as needed in order to maximize reduction in the model uncertainty. A mismatch between the supports of the model prior and the model likelihood indicates that the model has either too few or too many parameters, and there is not enough evidence to make the model selection with a high confidence [10].

Instead of choosing one best model, the predictions from multiple models can be combined. The joint prediction has potentially much larger discriminatory power than individual predictions [10]. Moreover, how well the model can describe the majority of outcomes from many random experiments under different experimental conditions can be as important as the model likelihood.

The information value of an experiment can be quantified as the Fisher information matrix with the entries

$${\left[\mathcal{I;}\left(\theta \right)\right]}_{i,j}=E\left[\left(\frac{\partial }{\partial {\theta }_i}logp\left(X;\theta \right)\right)\left(\frac{\partial }{\partial {\theta }_j}logp\left(X;\theta \right)\right)|\theta \right].$$E21

Since performing experiments is often costly, the following optimization objectives have been defined in the literature for linear models:

• A-optimality: minimize the average variance of parameter estimates;

• C-optimality: obtain the BLUE of linearly combined model parameters;

• D-optimality: maximize the entropy or determinant of information matrix;

• E-optimality: maximize the minimum eigenvalue of information matrix;

• T-optimality: maximize the trace of information matrix.

These optimality objectives may include prior distributions about the model and experiment parameters and combine other objectives for model selection.

The main challenge in evaluating the experiment design (20) is the computational complexity involved in searching the whole design space. The search strategies include linearization, local search, discretization, enumeration, approximation by regression or surrogate models, random (e.g., MCMC) sampling, genetic algorithms, as well as using Bayesian optimization methods, as will be discussed in Section 3.4.

Furthermore, the design (20) yields the single best experiment. In batch-optimum experiment design, $N$ experiments are performed simultaneously. The data from different experiments are conditionally independent given the parameters, $\theta$. The expected utility from these $N$ experiments will be, in general, different from the sum of utilities of individual experiments. A simpler objective is to minimize the predicted variance of observations, $x$, when the experiment designed as a single optimum is repeated $N$ times.

In the sequential experiment design, the posterior, $p\left(\theta ,x_t|d_t\right)$, from the $t$th experiment is used as a prior for the $\left(t+1\right)$th experiment, and then, the experiment conditions, $d_{t+1}$, are optimized. This is a greedy (sub-optimum) approach; however, it may still outperform batch design due to an inherent adaptation to $x_t$ and $d_t$. The optimum sequential experiment design is more complex, and it leads to a problem of dynamic programming.

3.2 Bayesian hypothesis testing

Let ${\theta }_0$ be the critical parameter value, so that the null hypothesis, ${\mathcal{H}}_0$, can be accepted, if the parameter value $\theta <{\theta }_0$, and it is rejected in favor of the alternative hypothesis, ${\mathcal{H}}_a$, otherwise. Denoting the loss function as $L\left(\theta ;{\mathcal{H}}_0\right)$ and $L\left(\theta ;{\mathcal{H}}_a\right)$ under the respective hypotheses, the Bayesian risk for observation, $x$, is computed as

$$R\left(x;{\theta }_0\right)={\int }_{\theta <{\theta }_0}L\left(\theta ;{\mathcal{H}}_0\right)\mathit{Pr}\left(x|\theta \right)p\left(\theta \right)\mathit{d\theta }+{\int }_{\theta \ge {\theta }_0}L\left(\theta ;{\mathcal{H}}_a\right)\mathit{Pr}\left(x|\theta \right)p\left(\theta \right)\mathit{d\theta }.$$E22

The optimum decision threshold is then obtained by minimizing the average risk, that is, ${\theta }_0^{\ast }={\mathit{\text{argmin}}}_{\theta }{E}_x\left[R\left(x;\theta \right)\right]$. The statistical power of the hypothesis test is determined by the sample size [11].

3.3 Bayesian machine learning

The general objective of machine learning is to learn how to label unseen data. In supervised learning, this objective is accomplished by training the model with labeled data. For instance, consider the minimum mean square estimation (MMSE) of label $Y$ for data $\mathbf{X}$ [3], that is,

$$f^{\ast }\left(\mathbf{x}\right)={\mathit{\text{argmin}}}_{f}E\left[{\left(Y,-,f,\left(\mathbf{X}=\mathbf{x}\right)\right)}^2\right]=E\left[Y|\mathbf{X}=\mathbf{x}\right].$$E23

The training data, $\left({\mathbf{X}}_1,Y_1\right),\left({\mathbf{X}}_2,Y_2\right),\dots ,\left({\mathbf{X}}_n,Y_n\right)$, are samples from the distribution, $p\left(Y,\mathbf{X}\right)=p\left(\mathbf{X}|Y\right)p\left(Y\right)$, assuming the likelihood, $p\left(\mathbf{X}|Y\right)$, and the prior, $p\left(Y\right)$, of data labels, $Y$. The MMSE expectation (23) can be then approximated as

$$E\left[Y|\mathbf{X}=\mathbf{x}\right]\approx \frac{1}{n}\sum _{i=1}^n{Y}_i\cdot {\mathcal{I;}}_{{\mathbf{X}}_i=\mathbf{x}}$$E24

assuming the indicator function $\mathcal{I}$.

In semi-supervised learning, there are additional unlabeled data, ${\mathbf{X}}_{n+1},{\mathbf{X}}_{n+2},\dots$, which can be used to estimate the distribution, $\int p\left(\mathbf{X},Y\right)\mathit{dY}$. It allows formulating and improving the estimates of the conditional expectation, $E\left[Y|\mathbf{X}=\mathbf{x}\right]$, and cast it as a missing labels problem. More generally, automated data labeling involves the problem of finding the discriminative data model, $p\left(Y|\mathbf{X}\right)$, or $p\left(\mathbf{X}|Y\right)$. However, the full generative model, $p\left(Y|\mathbf{X}\right)p\left(\mathbf{X}\right)=p\left(\mathbf{X}|Y\right)p\left(Y\right)$, may be required by some of the Bayesian inference methods outlined in the previous section.

Assuming data labels as random variables with their prior and posterior distributions is useful in tolerating label errors as well as missing labels. However, the challenge is that the assumed distributions may be biased or even incorrect. In addition to automated data labeling, Bayesian machine learning can generate more training data and improve the quality of training data by correcting the labels.

There is also an interesting connection to causal machine learning. In particular, the label can be assumed to be a cause of observed data (features) representing the effects. Since the probability, $\mathit{Pr}\left(\text{effect}\right)$, and the conditional probability, $\mathit{Pr}\left(\text{cause}|\text{effect}\right)$, are not independent, the data samples can be used to estimate their distribution, $p\left(\mathbf{X}\right)$, which in turn facilitates an anti-causal learning of causes (i.e., the data labels, $Y$). On the other hand, the probabilities $\mathit{Pr}\left(\text{cause}\right)$ and $\mathit{Pr}\left(\text{effect}|\text{cause}\right)$ are independent, and so, intuitively, they cannot be exploited for causal learning of data from their labels.

3.4 Bayesian optimization

Consider the task of minimizing or maximizing a function $f\left(x\right)$ over some high-dimensional feasible set, $A\subseteq {\mathrm{ℛ}}^n$. It is typically assumed that evaluating the function is numerically or otherwise very expensive, whereas testing whether $x\in A$ is computationally cheap. Even though $f$ is normally assumed to be continuous, its derivatives are not known. Furthermore, the problem is non-convex, as the function has no special structure, it usually contains many local optima, and its observations may be noisy.

The basic idea to solve such a difficult optimization problem is to learn a surrogate approximation of $f$, which is cheap to evaluate. The surrogate approximation is constructed in the context of Bayesian inference [12]. In particular, Bayesian optimization starts by evaluating the function $f$ at a few randomly chosen points, $x$. The obtained values, $f\left(x\right)$, are assumed to be samples of a random process. A Gaussian process (GP) is most commonly assumed, since it can yield closed-form expressions for the extrapolated values, even though inverting large covariance matrices is often numerically rather problematic. Then, the following steps are repeated as many times as can be practically afforded.

1. The posterior of a candidate sample is obtained given the existing samples $\left\{x_1,x_2,\dots ,x_n\right\}$;

2. The new sample is chosen as the one maximizing so-called acquisition function (to be discussed below);

3. The function $f$ is evaluated at the new sample, and the posterior is updated accordingly.

Assuming Bayesian regression over a Gaussian process, the $n$ existing samples, $f\left(x_1\right),\dots ,f\left(x_n\right)$, are normally distributed and have the joint mean, ${\mu }_0\left(x_{1:n}\right)$, and the covariance matrix, ${\mathrm{\Sigma }}_0\left(x_{1:n}\right)$. The corresponding posterior distribution of a candidate sample, $f\left(x\right)$, is also normal with the mean, ${\mu }_n\left(x\right)$, and the variance, ${\sigma }_n^2\left(x\right)$.

The acquisition function predicts the value of a new sample after already knowing the values, $f\left(x_1\right),\dots ,f\left(x_n\right)$. The values of the acquisition function tend to be larger when the calculated credible intervals are wider and when the posterior mean is larger. Since the mean, ${\mu }_n\left(x\right)$, is also a point estimate of $f\left(x\right)$ after $n$ observations, the value of $f\left(x\right)$ is estimated by interpolating ${\mu }_n\left(x\right)$; this corresponds to Gaussian regression. In other words, the measured values, $f\left(x_1\right),\dots ,f\left(x_n\right)$, are interpolated to estimate the mean, ${\mu }_n\left(x\right)$. The covariance matrix, ${\mathrm{\Sigma }}_n\left(x_{1:n}\right)$, determines how fast the function values vary in between the already measured samples. It is usually required that the values de-correlate with their distance. The covariance matrix can be also estimated from the existing samples as $\widehat{\eta }={\mathit{\text{argmax}}}_{\eta }p\left(f\left(x_{1:n}\right)|\eta \right)$ (maximum likelihood estimation, MLE) or as $\widehat{\eta }={\mathit{\text{argmax}}}_{\eta }p\left(f\left(x_{1:n}\right)|\eta \right)p\left(\eta \right)$ (maximum a posterior estimation, MAP). It is also possible to treat some parameters, $\eta$, as nuisance parameters and marginalize them from the likelihood assuming their prior, $p\left(\eta \right)$, before employing one of the Bayesian methods for intractable distributions.

Expected improvement (EI) is the most commonly assumed acquisition function. Provided that $f$ is observed without measurement noise, the current best choice is $f_n^{\ast }=max\left\{f\left(x_1\right),\dots ,f\left(x_n\right)\right\}$. The expected improvement is then defined as

$${\mathrm{EI}}_n\left(x\right)=E\left[{\left|f\left(x\right)-f_n^{\ast }\right|}_{+}\right]$$E25

where ${\left|\cdot \right|}_{+}\ge 0$, and the expectation is taken over by the posterior of $f\left(x\right)$ conditioned on the values known so far. The next best choice for sampling the function $f$ is

$$x_{n+1}={\text{argmax}}_x{\mathrm{EI}}_n\left(x\right).$$E26

Unlike $f$, the function ${\mathrm{EI}}_n\left(x\right)$ is inexpensive to evaluate, and its derivatives can be used to effectively maximize Eq. (26). The best expected improvement occurs at points far away from the previously evaluated points (the points having a large posterior variance) and at points having large posterior means; this represents an exploration–exploitation trade-off.

Knowledge gradient (KG) acquisition function selects the sampling point $x$ having the largest posterior mean, ${\mu }_n^{\ast }={max}_x{\mu }_n\left(x\right)$. It allows considering the posterior across the full domain of $f$ and how it is changed by a new sample. The knowledge gradient is computed as

$${\mathrm{KG}}_n\left(x\right)=E\left[{\mu }_{n+1}^{\ast }-{\mu }_n^{\ast }|x_{n+1}=x\right].$$E27

The next best sample corresponds to the maximum, that is,

$$x_{n+1}={\text{argmax}}_x{\mathrm{KG}}_n\left(x\right).$$E28

An efficient practical implementation of KG can assume multi-start stochastic gradient descent (or ascent, when the target function is to be maximized). This leads to a two-step maximization procedure when the maximum is first searched among a collection of candidate functions, and then, the selected function is maximized separately, for example, by differentiation.

In general, the KG acquisition function tends to significantly outperform the EI method, especially when the function observations are noisy.

Entropy search (ES) acquisition function assumes that the global optimum, $x^{\ast }$, is a random variable implied by the Gaussian process, $f\left(x\right)$. It performs the search for a new sample in order to obtain the largest decrease in differential entropy corresponding to the largest reduction of uncertainty about the global optimum. The ES acquisition function is defined as

$${\mathrm{ES}}_n\left(x\right)=H\left(P_n\left(x^{\ast }\right)\right)-E_{f\left(x\right)}\left[H\left(P_n(x^{\ast }|x,f\left(x\right))\right)\right]$$E29

where $H$ denotes entropy, so that the posterior across a full domain of $f$ and how it is changed by a new sample are again accounted for. This is useful when the observations are noisy. However, unlike the KG method, stochastic gradients cannot be obtained for the ES acquisition function.

Finally, the predictive entropy search (PES) acquisition function rewrites Eq. (29) using mutual information. Conceptually, it is exactly the same as the ES method; however, numerical properties of these two algorithms are different.

The basic optimization problem can be augmented by considering a set of objectives, $f\left(x,s\right)$, indexed by “fidelity”, $s$, such that lower values of $s$ mean higher fidelity, and $f\left(x,0\right)\equiv f\left(x\right)$ represents the original objective. For example, fidelity can represent the model order or the modeling granularity. There is also a cost, $c\left(s\right)$, associated with fidelity, $s$; the larger the required fidelity, the larger the cost. These problems are then referred to as multi-fidelity source evaluation, in the literature. The task is to maximize the function, $f\left(x,0\right)$, by observing a sequence of values of $f\left(x,s\right)$ at $n$ points, that is, $\left(x_1,s_1\right),\dots ,\left(x_n,s_n\right)$, subject to the total available budget, ${\sum }_{i=1}^{n}c\left(s_i\right)\le C_{\text{total}}$.

This problem can be further generalized by assuming that neither $f\left(x,s\right)$ nor $c\left(s\right)$ is monotonic in $s$; for example, keeping $s$ constant, the observed $f\left(x,s\right)$ across different regions of $x$ can have varying accuracy.

The problem of random environmental conditions is closely related to multi-task Bayesian optimization. It maximizes the function $\int f\left(x,w\right)p\left(w\right)\mathit{dw}$ by evaluating $f\left(x,w\right)$ for each $w$ at multiple values of $x$. It is assumed that evaluating $f\left(x,w\right)$ at both $x$ and $w$ is expensive, but evaluating $p\left(w\right)$ is cheap. The values $w$ represent random environmental conditions, and they act as noise in observations of $f\left(x,w\right)$. For example, $f$ can represent the average performance of a machine learning model with $w$-fold cross-validation. Another example are Bayesian neural networks having random coefficients with a certain probability distribution.

Comparing Bayesian optimization with other optimization methods, the latter usually work well only for specific problems or under specific conditions, but have high computational cost, and suffer from the curse of dimensionality, leading to a low sample efficiency and slow convergence. Bayesian optimization, on the other hand, can work very well for moderate dimensionality (the maximum dimension of about 20 is recommended in the literature), and the sampling decisions appear to be optimum despite being sequential. In addition, Bayesian optimization is considered to be derivative free but a black box method.

Despite a good performance, there is a need to provide better theoretical understanding of Bayesian optimization methods, define stopping rules, and improve learning of surrogate models, especially for non-Gaussian processes and for systems containing non-linear transformations. Bayesian optimization can be used to efficiently find hyperparameters of machine learning models [13]. Recently, maximization of a composite function, $g\left(h\left(x\right)\right)$, has been considered in the literature, where the inner function, $h$, is expensive to evaluate and a black box, whereas the outer function, $g$, is easy to evaluate.

4. Bayesian Monte Carlo simulations

Monte Carlo simulations are a general class of random sampling methods. They can be used to solve problems such as (2)–(4) using numerical algorithms. Simulations are often used to solve problems involving exploration of high-dimensional parameter spaces in order to perform sensitivity analysis, performance analysis, and system optimization. Unlike analytical approaches, computer simulations are more flexible, require fewer assumptions, and allow for more complex and, thus, more realistic models to be considered. In silico experiments are much more cost-effective than laboratory experiments conducted in life sciences. Digital twins are now often built to monitor complex engineering systems. Synthetic data from simulations can be used to train machine learning models and to avoid expensive or hard-to-get real-world measurements. However, there are often problems with reproducibility and validation of simulation results, especially when the simulation procedures and models become more complex. In general, it is extremely difficult to distinguish among algorithmic errors; incorrect choice and use of algorithms; the errors in mathematical derivations; conceptual errors, for example, due to invalid assumptions; and the errors in algorithm implementations.

Vast amounts of generated simulation data are usually reduced into a few summary statistics. This results in a loss of potentially valuable information, unless the statistics are also sufficient within a given context. There is rarely a deeper analysis of simulation techniques to explain why the implemented algorithms work or do not work, under what conditions, and whether they could be improved further. Thus, there is a need to enhance traditional Monte Carlo simulations and define the best practices of how to plan and conduct computer experiments in order to maximize the information gain and extract maximum useful knowledge from these simulations [14].

Improving the usefulness of Monte Carlo simulations requires that they are interpreted as statistical estimators of quantities of interest including identifying relevant causal associations [15, 16]. For instance, the sample statistics are point estimates, whereas it may be more useful to extract Bayesian posterior or credible intervals of these statistics when performing simulations. Bayesian analysis also enables counterfactual reasoning in order to find how the performance or system behavior would change, if the simulated model was altered. In particular, if the average performance of a system, $f\left(x\right)$, is evaluated as $\int f\left(x\right)p\left(x\right)\mathit{dx}$ over a distribution of factors, $p\left(x\right)$, then the counterfactual performance of a modified system having the factor distribution $p^{\ast }\left(x\right)$ can be defined as

$$\int f\left(x\right)p^{\ast }\left(x\right)\mathit{dx}=\int f\left(x\right)\underset{w\left(x\right)}{\underbrace{\frac{p^{\ast }\left(x\right)}{p\left(x\right)}}}p\left(x\right)\mathit{dx}.$$E30

Eq. (30) is akin to IS Eq. (5), that is, the system response, $f\left(x\right)$, is scaled by weights, $w\left(x\right)$, when assessing its hypothetical performance. Such approach can yield the following systematic and general strategy for designing systems via computer simulations:

1. Build a structural equation model (SEM) of the system using available performance data;

2. Use the SEM to define points of interventions and measurements;

3. Infer changes in the performance distribution under different hypothetical modifications;

4. Select and validate the modification with the best hypothetical improvement by implementing it in a real-world system.

A basic strategy for increasing the information gain of Monte Carlo simulations is to augment the number of observation points. In Figure 1A, the simulation inputs and outputs are represented by random variables, $X$ and $Y$, respectively, and the augmented output is represented by a random variable, $Z$. The corresponding posteriors and likelihoods are then computed as

Figure 1.

The input–output causal relationships in augmented Monte Carlo simulations.

$$p\left(X|Y\right)=\frac{p\left(Y|X\right)p\left(X\right)}{p\left(Y\right)}p\left(Y|X\right)=\sum _{Z}p\left(Y|X,Z\right)p\left(Z|X\right)$$E31

$$p\left(Z|Y\right)=\frac{p\left(Y|Z\right)p\left(Z\right)}{p\left(Y\right)}p\left(Y|Z\right)=\sum _{X}p\left(Y|X,Z\right)p\left(X|Z\right)$$E32

$$p\left(X|Z\right)=\frac{p\left(Z|X\right)p\left(X\right)}{p\left(Z\right)}p\left(Y,Z,X\right)=p\left(Y|Z,X\right)p\left(Z|X\right)p\left(X\right).$$E33

The causal relationships can be represented by probabilistic graphs, which are referred to as structural causal models (SCMs). These models are similar to SEM as well as allow capturing non-linear relationships. The SCMs are effectively Bayesian networks with directed edges indicating causal effects (since these effects are generally asymmetric) rather than statistical dependencies. The graph cycles are not allowed in order to avoid a variable to be a cause of itself. The endogenous noises are not shown explicitly.

The fundamental idea of studying causal effects is by accommodating interventions. Such a framework has been devised by Pearl [8], and it is now known as Do-calculus. The intervention denoted as $\mathrm{do}\left(X\right)$ enforces a certain deterministic or a random value from some distribution on the variable $X$. If this causes a change in the posterior, $p\left(Y|\mathrm{do}\left(X\right)\right)$, then $X$ and $Y$ are causally related. However, in general, observing $p\left(Y|X\right)$ is not sufficient to determine causality. The basic rules of Pearl’s Do-calculus are:

1. Observations can be inserted or deleted in conditional probabilities;

2. Actions and observations can be exchanged in conditional probabilities;

3. Actions can be inserted or deleted in conditional probabilities.

Moreover, if the causal effect is identifiable, then the causal effect statement can be transformed into a probability expression containing only the observed variables. Unknown causal dependencies can be replaced with conditional probabilities. The significance of these rules is that they allow performing causal inferences using traditional methods of statistical inference for Bayesian networks. This can be used for removing confounding bias, recovering from a selection bias, defining surrogate experiments, and extrapolating and transferring causal knowledge to other similar SCMs.

To illustrate the basic causal components of SCMs, consider three random variables, $A$, $B$, and $C$. The causal chain $A\to B\to C$ and the causal fork $A\leftarrow B\to C$ imply that $A$ and $C$ are generally dependent; however, they are independent, conditioned on $B$. On the other hand, the causal collider $A\to B\leftarrow C$ implies that $A$ and $C$ are generally independent, but they are dependent, conditioned on $B$. Thus, conditioning on $B$ in the case of causal chain or fork, blocks (or D-separates) the path between $A$ and $C$, whereas conditioning on $B$ opens up such a path in the case of causal collider.

These rules can be used to discuss causality in different SCM representations of augmented Monte Carlo simulation. In particular, the SCM in Figure 1B contains direct causal paths, $X\to Z$, $Z\to Y$, and $X\to Y$. There is also a backdoor path between $Z$ and $Y$, that is, $Z\leftarrow X\to Y$. Thus, $X$ is a common cause, or a confounder, so conditioning on $X$ blocks the backdoor path and enables causal inference.

The SCM in Figure 1C contains unmeasured or unobserved (e.g., exogenous, outside the model) variable, $U$. This is an example of so-called confounding by indication. There are two confounded associations, that is, $X\to Z\to Y$ and $U\to X\to Z\to Y$. Moreover, although conditioning on $U$ is not possible, conditioning on $X$ removes any unmeasured confounding. In the case of SCM in Figure 1D, conditioning on $X$ is sufficient to block the backdoor path.

The SCM in Figure 1E contains two unmeasured variables, $U_1$ and $U_2$. There is no bias without any conditioning. However, fixing $X$ as $\mathrm{do}\left(X\right)$ will induce a selection bias by opening the backdoor path, $Z\leftarrow U_2\to X\leftarrow U_1\to Y$, between $Z$ and $Y$. On the other hand, conditioning on $X$ will create direct as well as backdoor association between $Z$ and $Y$. The causal relationships in the remaining SCMs in Figures 1F–H are trivial.

The augmented Monte Carlo simulation in Figure 1A may require more sophisticated processing of the observed inputs and outputs as suggested in Figure 2 with the aim of improving the information gain and achieving explainability [17]. In particular, the summary statistics should be calculated between variables $X$ and $Y$ as in traditional Monte Carlo simulations and then also between $X$ and $Z$ and $Z$ and $Y$ in augmented Monte Carlo simulations. This should also improve detectability of events and anomalies. More importantly, it may be possible to define formal rules to automate building SEM and SCM models [18]. The models can be then fitted to observations and analyzed to discover causal relationships. The schematic of SEM or SCM is depicted in Figure 3, where U’s represent the inferred latent variables. The associations, $A_i$, interconnect the identified events, $E_j$, for example, using Bayesian rules of statistical dependencies and causal inferences.

Figure 2.

A structure of explainable Monte Carlo simulations.

Figure 3.

A schematic SCM for the explainable Monte Carlo simulation in Figure 2.

5. Conclusions

This chapter discussed the most important Bayesian frameworks and their applications. These frameworks play a crucial role in statistical and causal inferences. Learning from observations is one of the most fundamental concepts. Learning should be unbiased, systematic, replicable, generalizable, sample and other resources efficient and also provide a sufficient information gain. Since learning problems are often mathematically intractable, Bayesian reasoning can guide development of corresponding numerical algorithms. The algorithms can be considered to be concise representations of data, which can be generated by these algorithms. Bayesian frameworks allow building much more sophisticated models and carry out more complex computer simulations.

Bayesian experiment design is especially important due to a widespread adoption of computer simulations in designing and analyzing various real-world systems. Research in machine learning is slowly moving from statistical correlations-based models to models exploiting causal relationships. As discussed in this chapter, the causal relationships in SCM (or SEM) can be analyzed by statistical inference methods, which were developed for Bayesian probabilistic models. Bayesian machine learning assigns a distribution to data labels. It improves the performance of semi-supervised learning and enables automated data labeling and data label corrections. Bayesian optimization is particularly useful for optimizing real-world complex systems, which are very expensive to evaluate. Both Bayesian optimization and Bayesian machine learning are very active areas of research with many exciting open research problems. Bayesian Monte Carlo simulations aim at providing a systematic, semi-automated, explainable simulation framework. They augment the set of observations and leverage SEM or SCM in order to move from traditional descriptive simulations to predictive and even prescriptive simulations. Such research is still little explored in the existing literature.

Acknowledgments

This research was funded by a research grant from Zhejiang University.

Abbreviations

ABC	approximate Bayesian computation
BBVI	black box variational inference
BLUE	best linear unbiased estimators
CAVI	coordinated ascent variational inference
ES	entropy search
ELBO	evidence lower bound
EM	expectation maximization
EI	expected improvement
GP	Gaussian process
IS	importance sampling
KG	knowledge gradient
KL	Kullback–Leibler
	Monte Carlo
MAP	maximum a posterior estimation,
MCMC	Markov chain Monte Carlo
MLE	maximum likelihood estimation
MMSE	minimum mean square estimation
PES	predictive entropy search
RBM	restricted Boltzmann machines
SMC	sequential Monte Carlo
SVI	stochastic variational inference
SCM	structural causal model
SEM	structural equation model