**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 realworld 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 x*ð Þ, is evaluated as <sup>Ð</sup> *f x*ð Þ*p x*ð Þ*dx* over a distribution of factors, *p x*ð Þ, then the counterfactual performance of a modified system having the factor distribution *<sup>p</sup>*<sup>∗</sup> ð Þ *<sup>x</sup>* can be defined as

$$\int f(\mathbf{x}) p^\*\left(\mathbf{x}\right) d\mathbf{x} = \int f(\mathbf{x}) \underbrace{\frac{p^\*\left(\mathbf{x}\right)}{p(\mathbf{x})} p(\mathbf{x})}\_{w(\mathbf{x})} d\mathbf{x}.\tag{30}$$

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


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

$$p(X|Y) = \frac{p(Y|X)p(X)}{p(Y)}p(Y|X) = \sum\_{Z} p(Y|X,Z)p(Z|X) \tag{31}$$

$$p(Z|Y) = \frac{p(Y|Z)p(Z)}{p(Y)}p(Y|Z) = \sum\_{X} p(Y|X,Z)p(X|Z) \tag{32}$$

$$p(X|Z) = \frac{p(Z|X)p(X)}{p(Z)}p(Y,Z,X) = p(Y|Z,X)p(Z|X)p(X). \tag{33}$$

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 Docalculus. The intervention denoted as doð Þ *X* enforces a certain deterministic or a random value from some distribution on the variable *X*. If this causes a change in the posterior, *p Y*ð Þ jdoð Þ *X* , then *X* and *Y* are causally related. However, in general,

**Figure 1.** *The input–output causal relationships in augmented Monte Carlo simulations.* observing *p Y*ð Þ j*X* is not sufficient to determine causality. The basic rules of Pearl's Do-calculus are:


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* ! *B* ! *C* and the causal fork *A B* ! *C* imply that *A* and *C* are generally dependent; however, they are independent, conditioned on *B*. On the other hand, the causal collider *A* ! *B 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* ! *Z*, *Z* ! *Y*, and *X* ! *Y*. There is also a backdoor path between *Z* and *Y*, that is, *Z X* ! *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* ! *Z* ! *Y* and *U* ! *X* ! *Z* ! *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*<sup>1</sup> and *U*2. There is no bias without any conditioning. However, fixing *X* as doð Þ *X* will induce a selection bias by opening the backdoor path, *Z U*<sup>2</sup> ! *X U*<sup>1</sup> ! *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

**Figure 2.** *A structure of explainable Monte Carlo simulations.*

**Figure 3.**

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

associations, *Ai*, interconnect the identified events, *Ej*, for example, using Bayesian rules of statistical dependencies and causal inferences.
