# Neural Networks Assisted Metropolis-Hastings for Bayesian Estimation of
Critical Exponent on Elliptic Black Hole Solution in 4D Using Quantum
Perturbation Theory

Armin<EMAIL_ADDRESS>Department of Mathematics and Statistics,
Memorial University of Newfoundland, St. John’s, NL, Canada.  Ehsan
<EMAIL_ADDRESS>Roberto J<EMAIL_ADDRESS>
###### Abstract

It is well-known that the critical gravitational collapse produces continuous
self-similar solutions characterized by the Choptuik critical exponent,
$\gamma$. We examine the solutions in the domains of the linear perturbation
equations, considering the numerical measurement errors. Specifically, we
study quantum perturbation theory for the four-dimensional Einstein-axion-
dilaton system of the elliptic class of $\text{SL}(2,\mathbb{R})$
transformations. We develop a novel artificial neural network-assisted
Metropolis-Hastings algorithm based on quantum perturbation theory to find the
distribution of the critical exponent in a Bayesian framework. Unlike existing
methods, this new probabilistic approach identifies the available
deterministic solution and explores the range of physically distinguishable
critical exponents that may arise due to numerical measurement errors.

##  1 Introduction

Black holes are known to be the final state of gravitational collapse. They
have a well-known characteristic property. They are entirely defined by mass,
angular momentum, and charge. Choptuik [1] shows that there seems to be
another parameter or the fourth quantity that establishes the collapse itself.
Following up on Christodoulou’s work on the spherically symmetric collapse of
scalar fields [2], Choptuik looked into the idea that a critical behaviour can
show discrete spacetime self-similarity. The amplitude of the scalar field
fluctuation shows that $p$ must be greater than the critical value
$p_{\text{crit}}$ to form a black hole. Moreover, for values of $p$ above this
threshold, the mass of the black hole $M_{\text{bh}}$ results a scaling law as

$r_{S}(p)\propto M_{\text{bh}}(p)\propto(p-p_{\text{crit}})^{\gamma}\,.$ (1)

The critical exponent in 4d for a single real scalar field is given by
$\gamma\simeq 0.37$ [1, 3, 4], while for general dimension ($d\geq 4$) [5, 6]
these relations are modified as

$r_{S}(p)\propto(p-p_{\text{crit}})^{\gamma}\,,\quad
M_{\text{bh}}(p)\sim(p-p_{\text{crit}})^{(D-3)\gamma}.$ (2)

Various numerical simulations with other fields are carried out [7, 8, 9, 10,
11, 12]. As an instance the collapse of a perfect fluid is done in [13, 14, 5,
15]. In [14] they explored $\gamma\simeq 0.36$ and therefore it was realised
in [16] that $\gamma$ might be universal for all matter field coupled to
gravity in four dimensions. It was first in [5, 15, 17] that the authors found
the critical exponent may be found by working out the perturbations of the
solutions. To achieve the perturbations, one needs to perturb any field $h$
(be it the metric or the matter content) as

$h=h_{0}+\varepsilon\,h_{-\kappa},$ (3)

where the perturbation $h_{-\kappa}$ does have the scaling
$-\kappa\in\mathbb{C}$ that labels the different modes. Within the possible
values of $\kappa$, we express the most relevant mode $\kappa^{*}$ as the
highest value of $\real(\kappa)$.

It was argued in [5, 15, 17] that the $\kappa^{*}$ is related to the critical
exponent by $\gamma=\frac{1}{\real\kappa^{*}}$. The axial symmetry was
analyzed in [18], and the critical solution for the shock waves was
investigated by [19].

For the first time, the axion-dilaton critical collapse solution coupled to
gravity in four dimensions was determined by [20] where they found the value
$\gamma\simeq 0.2641$, hence arguing the doubts relating to the universality
of $\gamma$ in four dimensions. The first motivation for studying critical
solutions in the axion-dilaton system is the AdS/CFT correspondence [21],
which relates the Choptuik exponent to the imaginary part of quasinormal
modes, and the dual conformal field theory [22]. Additional motivations
include the holographic description of black hole formation [6] and the
broader physics of black holes and their applications to holography and string
theory [23, 24].

In type IIB string theory, there is significant interest in exploring
gravitational collapse in spaces that can asymptotically approach
$AdS_{5}\times S^{5}$ where the matter content is described by the axion-
dilaton system and the self-dual 5-form field. The paper [25] recently
investigated entire families of distinguishable continuous self-similar
solutions of the Einstein-axion-dilaton system in four and five dimensions,
including all three conjugacy classes of $\mathrm{SL}(2,\mathbb{R})$. This
study builds upon works of [26, 27]. By applying advanced analytic and
numerical methods described in [28], the perturbed critical solution of a
four-dimensional elliptic critical collapse was computed and determined the
value of $\gamma$ to be around 0.2641, as previously reported by [20]. The
findings provide strong confidence in our ability to determine other critical
exponents across different dimensions and for various classes of solutions.

The study by [29] utilized regression models to estimate nonlinear critical
functions. Subsequently, [30] introduced several methods, including truncated
power basis, natural spline, and penalized B-spline regression models, to
estimate the nonlinear functions relevant to black hole physics, specifically
in the axion-dilaton case. Recently, in [31], artificial neural networks were
applied to find black hole solutions within the parabolic class in higher
dimensions. Additionally, in [32], the Hamiltonian Monte Carlo method was
proposed to analyze the complexity of elliptic black holes within a Bayesian
framework. Lastly, [33] employed the sequential Monte Carlo approach to
investigate the multimodal posterior distributions of critical functions in
hyperbolic equations. This probabilistic approach helps identify existing
solutions in the literature and finds all possible solutions that might arise
due to measurement errors.

Unlike all existing methods in the literature, we treat the critical exponent
in this work as a random variable for the first time, thereby better capturing
the uncertainty and numerical measurement errors associated with the perturbed
equations of motion. To achieve this, we have developed a novel statistical
method that integrates the power and flexibility of Markov Chain Monte Carlo
(MCMC) and artificial neural networks, estimating the posterior probability
distribution of the critical exponent within a Bayesian framework. Our
approach introduces a generic formalism based on artificial neural networks
and the Metropolis-Hastings algorithm [34, 35] to estimate critical exponent,
addressing their inherent complexity. Using quantum perturbation theory, we
apply the Bayesian estimation method to determine the critical exponent for
the elliptic black hole solution in 4d, estimating the Choptuik exponent
within the domain of the equations. We investigate the range of possible
values for the critical exponent in the Einstein-axion-dilaton system’s 4d
elliptic black hole solution. In this novel iterative approach, each iteration
combines the Metropolis-Hastings algorithm with artificial neural networks to
simultaneously solve the nonlinear perturbed equations of motion and identify
the stochastically most likely values of the critical exponent. We have
investigated the Bayesian estimation of the critical exponent based on various
perturbed equations independently and by considering all the perturbed
equations simultaneously. This was done to examine the impact of various
perturbed equations on stochastic accept-or-reject transitions. Unlike
traditional methods, this new probabilistic approach provides the established
solution and explores the range of physically distinguishable critical
exponents that may arise due to numerical measurement errors.

This paper is structured as follows: Section 2 explores the Einstein-axion-
dilaton system. Sections 3 details the quantum perturbation analysis and the
perturbed equations of motion. Section 4 presents the statistical methods
employed, including polynomial regression, artificial neural networks, and
Markov Chain Monte Carlo. Section 5 provides the numerical studies that
examine the posterior distribution and Bayesian estimation of the critical
exponent.

##  2 The Einstein-axion-dilaton System

The Einstein-axion-dilaton system coupled to gravity in $d$ dimensions is
defined by the following action:

$S=\int
d^{d}x\sqrt{-g}\left(R-\frac{1}{2}\frac{\partial_{a}\tau\partial^{a}\bar{\tau}}{(\imaginary\tau)^{2}}\right)\;,$
(4)

which is known by the effective action of type II string theory [36, 37]. The
axion-dilaton is expressed by $\tau\equiv a+ie^{-\phi}$. This action enjoys
the $\mathrm{SL}(2,\mathbb{R})$ symmetry

$\tau\rightarrow M\tau\equiv\frac{\alpha\tau+\beta}{\gamma\tau+\delta}\;,$ (5)

where $\alpha$, $\beta$, $\gamma$, $\delta$ are real parameters satisfying
$\alpha\delta-\beta\gamma=1$. If the quantum effects are considered then the
$\mathrm{SL}(2,\mathbb{R})$ symmetry gets exchanged to
$\mathrm{SL}(2,\mathbb{Z})$ and this S-duality was revealed to be the non-
perturbative symmetry of IIB string theory [38, 39, 40]. Taking the above
action, one derives the equations of motion as

$R_{ab}=\tilde{T}_{ab}\equiv\frac{1}{4(\imaginary\tau)^{2}}(\partial_{a}\tau\partial_{b}\bar{\tau}+\partial_{a}\bar{\tau}\partial_{b}\tau)\,,$
(6)
$\nabla^{a}\nabla_{a}\tau+\frac{i\nabla^{a}\tau\nabla_{a}\tau}{\imaginary\tau}=0\,.$
(7)

We take the spherical symmetry for both background and perturbations and the
general form of the metric in $d$ dimensions is given by

$ds^{2}=(1+u(t,r))(-b(t,r)^{2}dt^{2}+dr^{2})+r^{2}d\Omega_{d-2}^{2}\,,$ (8)
$\tau=\tau(t,r).$ (9)

The angular part of the metric in $d$ dimensions is $d\Omega_{d-2}^{2}$. One
can find out a scale-invariant solution by requiring the fact that under a
scale transformation $(t,r)\rightarrow(\Lambda t,\Lambda r)$ so that the line
element can be changed as $ds^{2}\rightarrow\Lambda^{2}ds^{2}$. Hence all the
functions must be scale invariant, i.e. $u(t,r)=u(z)$, $b(t,r)=b(z)$,
$z\equiv-r/t$. The effective action (4) is invariant under the
$\mathrm{SL}(2,\mathbb{R})$ transformation (5), hence $\tau$ also does need to
be invariant up to an $\mathrm{SL}(2,\mathbb{R})$ transformation, which means
that if the following holds

$\tau(\Lambda t,\Lambda r)=M(\Lambda)\tau(t,r),$ (10)

we then call a system of $(g,\tau)$, which respects the above properties of a
continuous self-similar (CSS) solution. It is worth highlighting that various
cases get related to other classes of
$\evaluated{\derivative{M}{\Lambda}}_{\Lambda=1}$ [25]. Hence, $\tau$ can take
three different forms [26], called elliptic, hyperbolic and parabolic cases.
In this paper, we deal with the elliptic ansatz and its form for the axion-
dilaton case can be cast as follows

$\displaystyle\tau(t,r)=i\frac{1-(-t)^{i\omega}f(z)}{1+(-t)^{i\omega}f(z)}\,,$
(11)

where $\omega$ is a real parameter to be known and the function $f(z)$ needs
to satisfy $\absolutevalue{f(z)}<1$ for the elliptic case. By replacing the
CSS ansätze in the equations of motion, we would be able to get a differential
system of equations for $u(z)$, $b(z)$, $f(z)$. Making use of spherical
symmetry, one can show that $u(z),u^{\prime}(z)$ can be expressed in terms of
$b(z)$ and $f(z)$ so that finally, we are left with some ordinary differential
equations (ODEs) as follows

$\displaystyle b^{\prime}(z)$ $\displaystyle=B(b(z),f(z),f^{\prime}(z))\,,$
(12) $\displaystyle f^{\prime\prime}(z)$
$\displaystyle=F(b(z),f(z),f^{\prime}(z))\,.$ (13)

These equations in elliptic 4d can be written in a closed form as

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
b^{\prime}+{z(b^{2}-z^{2})\over
b(-1+|f|^{2})^{2}}f^{\prime}\bar{f}^{\prime}-{i\omega(b^{2}-z^{2})\over
b(-1+|f|^{2})^{2}}(f\bar{f}^{\prime}-\bar{f}f^{\prime})-{\omega^{2}z|f|^{2}\over
b(-1+|f|^{2})^{2}},$ $\displaystyle 0$ $\displaystyle=$ $\displaystyle
f^{\prime\prime}-{z(b^{2}+z^{2})\over b^{2}(-1+|f|^{2})^{2}}f^{\prime
2}\bar{f}^{\prime}+{2\over(1-|f|^{2})}\left(1-{i\omega(b^{2}+z^{2})\over
2b^{2}(1-|f|^{2})}\right)\bar{f}f^{\prime 2}$ (14)
$\displaystyle+{i\omega(b^{2}+2z^{2})\over
b^{2}(-1+|f|^{2})^{2}}ff^{\prime}\bar{f}^{\prime}+{2\over z}\left(1+{i\omega
z^{2}(1+|f|^{2})\over(b^{2}-z^{2})(1-|f|^{2})}\right.$
$\displaystyle+\left.{\omega^{2}z^{4}|f|^{2}\over
b^{2}(b^{2}-z^{2})(1-|f|^{2})^{2}}\right)f^{\prime}+{\omega^{2}z\over
b^{2}(-1+|f|^{2})^{2}}f^{2}\bar{f}^{\prime}+$
$\displaystyle{2i\omega\over(b^{2}-z^{2})}\left(\frac{1}{2}-{i\omega(1+|f|^{2})\over
2(1-|f|^{2})}\right.-\left.{\omega^{2}z^{2}|f|^{2}\over
2b^{2}(-1+|f|^{2})^{2}}\right)f.$

The above equations have five singularities [26] located at $z=\pm 0$,
$z=\infty$ and $z=z_{\pm}$ where the last singularities are known by
$b(z_{\pm})=\pm z_{\pm}$. $z=z_{+}$ is just a mere coordinate singularity [20,
26], thus $\tau$ should have been regular across it and this constraint
actually translates to having the finite value for $f^{\prime\prime}(z)$ as
$z\rightarrow z_{+}$.

One can actually realize that the vanishing of the divergent part of
$f^{\prime\prime}(z)$ provides us with an equation that is a complex valued
constraint at $z_{+}$ which can be indicated by
$G(b(z_{+}),f(z_{+}),f^{\prime}(z_{+}))=0$ where the final form of $G$ for the
elliptic class is given as [25]

$\displaystyle G(f(z_{+}),f^{\prime}(z_{+}))=$
$\displaystyle\,2z\bar{f}(z_{+})\left(-2\omega^{2}\right)f^{\prime}(z_{+})$
(15)
$\displaystyle+f(z_{+})\bar{f}(z_{+})\left(2z_{+}\bar{f}(z_{+})(-2+2i\omega+2)f^{\prime}(z_{+})+2i\omega\left(2+\omega^{2}\right)\right)$
$\displaystyle-\frac{2z_{+}(2+2i\omega-2)f^{\prime}(z_{+})}{f(z_{+})}$
$\displaystyle+2\omega(\omega-i)f(z_{+})^{2}\bar{f}(z_{+})^{2}-2\omega(\omega+i)\,.$

The initial conditions are known by the smoothness of the solution. By
applying polar coordinate as $f(z)=f_{m}(z)e^{if_{a}(z)}$, we can see that all
equations are invariant under a global redefinition of the phase of $f(z)$, so
this means that $f_{a}(0)=0$. Having set the regularity condition at $z=0$ and
making use of the residual symmetries in the equations of motions (14), we
reveal the initial boundary conditions of the equations of motion as follows

$\displaystyle b(0)=1,f_{m}(0)=x_{0},\quad\quad\quad\quad
f_{m}^{\prime}(0)=f_{a}^{\prime}(0)=f_{a}(0)=0\;,$ (16)

where $x_{0}$ is a real parameter and $(0<x_{0}<1)$. Therefore, we have two
constraints (the vanishing of the real and imaginary parts of $G$) and two
parameters $(\omega,x_{0})$. The discrete solution in four dimensions was
found in [28]. Note that these solutions are explicitly found by root-finding
procedure as well as numerically integrating all the equations of motion. For
instance, the solution in 4d elliptic case was discovered to be [41, 26] as

$\omega=1.176,\quad\absolutevalue{f(0)}=0.892,\quad z_{+}=2.605$ (17)

where the black hole solutions in higher dimensions for the axion-dilaton
system are obtained in [42].

##  3 Quantum Perturbation Analysis

In this section, we would like to derive the quantum perturbation equations
for black holes in the elliptic class in four dimensions.

This method can be implemented as an extensive method that holds for all
dimensions and all matter content. Some of the steps are taken from [43] and
[44] while we could remove $u(t,r)$ and its derivatives from all the
equations, making use of some algebraic computations444Some similar
perturbations of spherically symmetric solutions for Horava Gravity were
carried out in [45].. One can perturb the exact solutions $h_{0}$ (where $h$
denotes either $b$, or $f$) that are explored in Section 2 as

$h(t,r)=h_{0}(z)+\varepsilon\,h_{1}(t,r)$ (18)

where $\varepsilon$ is a small number. By expanding equations in powers of
$\varepsilon$, the zeroth order part results in the background equations which
have been studied in Section 2 and the linearized equations for the
perturbations $h_{1}(t,r)$ are explored due to the linear terms in
$\varepsilon$. One can consider the perturbations of the form

$\displaystyle h(z,t)=h_{0}(z)+\varepsilon(-t)^{-\kappa}h_{1}(z).$ (19)

The four-dimensional axion-dilaton system is known to be stable, and we
consider $\real\kappa>0$. We can explore the spectrum of $\kappa$ by solving
the equations for $h_{1}(z)$. Indeed, one can point out that the general
solution to the first-order equations will be gained with the linear
combination of these modes. We would like to find out the mode $\kappa^{*}$
with the biggest real part (by assuming a growing mode, which would be
$t\rightarrow 0$, i.e. $\real\kappa>0$), which will be related to critical
Choptuik exponent as [5, 15, 17]

$\displaystyle\gamma=\frac{1}{\real\kappa^{*}}.$ (20)

We consider only real modes $\kappa^{*}$ where the values $\kappa=0$ and
$\kappa=1$ are indeed gauge modes with respect to phase of $f$ or $U(1)$ re-
definitions of $f$ and time translations as well [28] and these modes have
been eliminated from the calculations.

###  3.1 Perturbation Equations in 4D For the Elliptic Class

The derivation of the perturbations has been explained in [28]. However, we
want to analyze the quantum perturbation to explore the entire range of
perturbed equations independently and by considering all the perturbed
equations simultaneously, as well as their impact on stochastic accept-reject
transitions in finding the posterior distributions and Bayesian estimate of
the critical exponent. We here highlight briefly all perturbation equations.
The perturbation ansatz (19), for $b$, $\tau$ functions are

$b(t,r)=b_{0}(z)+\varepsilon\,(-t)^{-\kappa}b_{1}(z)\ ,$ (21)
$\tau(t,r)=i\frac{1-(-t)^{i\omega}f(t,r)}{1+(-t)^{iw}f(t,r)}\ ,$ (22)
$f(t,r)\equiv f_{0}(z)+\varepsilon(-t)^{-\kappa}f_{1}(z).$ (23)

Using the ansätze (21) for the metric functions, we can calculate the Ricci
tensor for the metric as a function of $\varepsilon$ to find the zeroth-order
and first-order parts from the limiting behaviours as

$R^{(0)}_{ab}=\lim_{\varepsilon\rightarrow 0}R_{ab}(\varepsilon)\,\quad
R^{(1)}_{ab}=\lim_{\varepsilon\rightarrow
0}\derivative{R_{ab}(\varepsilon)}{\varepsilon}\,.$ (24)

The same method is applied to the right-hand side $\tilde{T}_{ab}$ of the
field equations, which results in

$\tilde{T}^{(0)}_{ab}=\lim_{\varepsilon\rightarrow
0}\tilde{T}_{ab}(\varepsilon)\,\quad\tilde{T}^{(1)}_{ab}=\lim_{\varepsilon\rightarrow
0}\derivative{\tilde{T}_{ab}(\varepsilon)}{\varepsilon}.$ (25)

Indeed, the Einstein Field Equations (EFEs) are held order by order so that

$R^{(0)}_{ab}=\tilde{T}^{(0)}_{ab}\,,\quad R^{(1)}_{ab}=\tilde{T}^{(1)}_{ab}.$
(26)

As explained, one can remove $u$ and its first derivative in terms of
$b_{0},f$ and their first derivatives. We now combine the field equations so
that we would be able to eliminate the second-derivative terms in $b(t,r)$.
This procedure is known as Hamiltonian constraint, which can be shown by

$C(\varepsilon)\equiv
R_{tt}+b^{2}\,R_{rr}-\tilde{T}_{tt}-b^{2}\,\tilde{T}_{rr}=0\,.$ (27)

One finds the lowest-order contribution as a first-order equation relating
$b_{0}^{\prime}$ to $b_{0}$, $f_{0}$, $f_{0}^{\prime}$ as follows

$\displaystyle
b_{0}^{\prime}=\frac{((z^{2}-b_{0}^{2})f_{0}^{\prime}(z\bar{f}_{0}^{\prime}+i\omega\bar{f}_{0})+\omega
f_{0}(\omega
z\bar{f}_{0}-i(z^{2}-b_{0}^{2})\bar{f}_{0}^{\prime}))}{b_{0}(1-f_{0}\bar{f}_{0})^{2}}\,.$
(28)

In a similar way, the first correction is given by

$\evaluated{\derivative{C(\varepsilon)}{\varepsilon}}_{\varepsilon=0}=0\ ,$
(29)

which is a first-order equation relating $b_{1}^{\prime}$ to $b_{0}$,
$b_{0}^{\prime}$ $f_{0}$, $f_{0}^{\prime}$, $f_{0}^{\prime\prime}$, and the
other perturbations $b_{1}$, $f_{1}$, $f_{1}^{\prime}$, which are indeed
linear in all perturbations. For the elliptic class in 4d one obtains the
linearized equations for $b_{1}^{\prime}$ as follows

$\displaystyle(L_{1})b_{1}^{\prime}$
$\displaystyle=rt((t-2t)b_{0}+rb_{0}^{\prime})\bigg{(}-\frac{2f_{1}^{\prime}\bar{f}_{0}^{\prime}r^{2}}{t^{4}s_{0}^{2}}-\frac{2b_{1}b_{0}^{\prime}}{rt}+\frac{\kappa
2b_{0}b_{1}b_{0}^{\prime}}{r\left((t-2t)b_{0}+rb_{0}^{\prime}\right)}$ (30)
$\displaystyle\quad-\frac{4i\omega
b_{0}b_{1}\bar{f}_{0}f_{0}^{\prime}}{rts_{0}^{2}}+\frac{2i\omega
b_{0}^{2}\bar{f}_{1}f_{0}^{\prime}}{rts_{0}^{3}}+\frac{2\kappa\bar{f}_{1}f_{0}^{\prime}r}{t^{3}s_{0}^{2}}+\frac{2i\omega\bar{f}_{1}f_{0}^{\prime}r}{t^{3}s_{0}^{2}}$
$\displaystyle\quad-\frac{2\kappa
b_{0}^{2}\bar{f}_{1}f_{0}^{\prime}}{rts_{0}^{2}}+\frac{2i\omega\bar{f}_{0}f_{1}^{\prime}r}{t^{3}s_{0}^{2}}-\frac{2i\omega
b_{0}^{2}\bar{f}_{0}f_{1}^{\prime}}{rts_{0}^{2}}+\frac{4b_{0}b_{1}f_{0}^{\prime}\bar{f}_{0}^{\prime}}{t^{2}s_{0}^{2}}$
$\displaystyle\quad+\frac{2b_{0}^{2}f_{1}^{\prime}\bar{f}_{0}^{\prime}}{t^{2}s_{0}^{2}}+\frac{1}{rt^{4}s_{0}^{3}}2f_{1}\Big{(}t\omega\bar{f}_{0}^{2}(rt(-i\kappa+\omega)f_{0}-2i(r^{2}-t^{2}b_{0}^{2})f_{0}^{\prime})$
$\displaystyle\quad\quad+t(\kappa-i\omega)(-r^{2}+t^{2}b_{0}^{2})\bar{f}_{0}^{\prime}+\bar{f}_{0}\big{(}rt^{2}\omega(i\kappa+\omega)+t(\kappa+i\omega)$
$\displaystyle\quad\quad\times(r^{2}-t^{2}b_{0}^{2})f_{0}\bar{f}_{0}^{\prime}+2r(r^{2}-t^{2}b_{0}^{2})f_{0}^{\prime}\bar{f}_{0}^{\prime}\big{)}\Big{)}-\frac{2f_{0}^{\prime}\bar{f}_{1}^{\prime}r^{2}}{t^{4}s_{0}^{2}}+\frac{2b_{0}^{2}f_{0}^{\prime}\bar{f}_{1}^{\prime}}{t^{2}s_{0}^{2}}$
$\displaystyle\quad-\frac{2if_{0}(\bar{f}_{1}(r\omega
t^{2}(\kappa+i\omega)+\omega
t(2r^{2}-t^{2}b_{0}^{2})\bar{f}_{0}f_{0}^{\prime}+2ir(r^{2}-t^{2}b_{0}^{2})f_{0}^{\prime}\bar{f}_{0}^{\prime}))}{rt^{4}s_{0}^{3}}$
$\displaystyle\quad-\frac{2if_{0}t\omega(-\bar{f}_{1}^{\prime}r^{2}+2t^{2}b_{0}b_{1}\bar{f}_{0}^{\prime}+t^{2}b_{0}^{2}\bar{f}_{1}^{\prime})}{rt^{4}s_{0}^{3}}+\frac{4i\omega
f_{0}^{2}\left(r^{2}-t^{2}b_{0}^{2}\right)\bar{f}_{1}\bar{f}_{0}^{\prime}}{rt^{3}s_{0}^{3}}$
$\displaystyle+\frac{2i\omega
f_{0}^{2}\left(\bar{f}_{0}\left(-\bar{f}_{1}^{\prime}r^{2}+t(\kappa-i\omega)\bar{f}_{1}r+2t^{2}b_{0}b_{1}\bar{f}_{0}^{\prime}+t^{2}b_{0}^{2}\bar{f}_{1}^{\prime}\right)\right)}{rt^{3}s_{0}^{3}}\bigg{)}\,,$

where

$L_{1}=2b_{0}((\kappa-1)tb_{0}+rb_{0}^{\prime}),\quad
s_{0}=(f_{0}\bar{f}_{0}-1).$ (31)

It can be seen that (30) is invariant under dilations
$(t,r)\rightarrow(e^{\lambda}r,e^{\lambda}t)$, so that, by changing
coordinates to $(t,z)$, one finds that all the factors of $t$ will cancel off
and the final result just depends on $z$. Hence, we are left with just real
and linear ordinary differential equations.

If we apply the perturbation ansätze in the $\tau$ equation of motion (7) and
extract the zeroth and first-order parts and remove all $u_{0}^{\prime}$,
$u_{1}$, $u_{1}^{\prime}$, then the resulting zeroth-order part will be an
equation including $b_{0}$, $b_{0}^{\prime}$ $f_{0}$, $f_{0}^{\prime}$,
$f_{0}^{\prime\prime}$, while substituting $b_{0}^{\prime}$ according to (12)
and solving for $f_{0}^{\prime\prime}$, we obtain a second-order equation for
$f_{0}$, which is indeed the explicit form of (13). The first-order part
includes $b_{1}$, $b_{1}^{\prime}$, $f_{1}$, $f_{1}^{\prime}$, and
$f_{1}^{\prime\prime}$ linearly, while the equations depend on the zeroth-
order functions as well as their derivatives non-linearly as background
functions. Hence, we are left with the perturbations for
$f_{1}^{\prime\prime}$ given by

$\displaystyle(L_{3})f_{1}^{\prime\prime}=$ $\displaystyle
r^{2}\bigg{(}it\omega
f_{0}^{2}(\bar{f}_{1}b_{0}^{\prime}+\bar{f}_{0}b_{1}^{\prime})+f_{1}b_{0}^{\prime}(\kappa
t-it\omega-t(\kappa-2i\omega)f_{0}\bar{f}_{0}+rf_{0}^{\prime}\bar{f}_{0})$
(32)
$\displaystyle\quad-r(b_{1}^{\prime}f_{0}^{\prime}+b_{0}^{\prime}f_{1}^{\prime})+f_{0}\big{(}b_{1}^{\prime}(-it\omega+r\bar{f}_{0}f_{0}^{\prime})+rb_{0}^{\prime}(\bar{f}_{1}f_{0}^{\prime}+\bar{f}_{0}f_{1}^{\prime})\big{)}\bigg{)}$
$\displaystyle+rt^{2}b_{0}^{2}\big{(}f_{1}\bar{f}_{0}b_{0}^{\prime}f_{0}^{\prime}-b_{1}^{\prime}f_{0}^{\prime}-b_{0}^{\prime}f_{1}^{\prime}+f_{0}(\bar{f}_{1}b_{0}^{\prime}f_{0}^{\prime}+\bar{f}_{0}(b_{1}^{\prime}f_{0}^{\prime}+b_{0}^{\prime}f_{1}^{\prime}))\big{)}$
$\displaystyle-t^{2}b_{0}^{3}\bigg{(}2r\bar{f}_{1}f_{0}^{\prime
2}-2tf_{1}^{\prime}+4r\bar{f}_{0}f_{0}^{\prime}f_{1}^{\prime}+f_{1}\bar{f}_{0}(2tf_{0}^{\prime}-rf_{0}^{\prime\prime})$
$\displaystyle\quad+f_{0}\big{(}2t\bar{f}_{0}f_{1}^{\prime}+\bar{f}_{1}(2tf_{0}^{\prime}-rf_{0}^{\prime\prime})\big{)}\bigg{)}$
$\displaystyle-
rb_{0}\bigg{(}t^{2}\omega(\omega-i)f_{0}^{2}\bar{f}_{1}+2r\big{(}-r\bar{f}_{1}f_{0}^{\prime
2}+(t(1+\kappa-i\omega)-2r\bar{f}_{0}f_{0}^{\prime})f_{1}^{\prime}\big{)}$
$\displaystyle\quad\quad+f_{1}\Big{(}t^{2}(-\kappa^{2}+\kappa(-1+2i\omega)+\omega(i+\omega))+t^{2}(\kappa+\kappa^{2}+2i\kappa\omega$
$\displaystyle\quad\quad\quad+2\omega(-i+\omega)\big{)}\bar{f}_{0}f_{0}+r\bar{f}_{0}(2t(-1+2\kappa-i\omega)f_{0}^{\prime}+rf_{0}^{\prime\prime})\Big{)}$
$\displaystyle\quad\quad+rf_{0}(-2t(1+\kappa+i\omega)\bar{f}_{0}f_{1}^{\prime}+\bar{f}_{1}(-2it(-i+\omega)f_{0}^{\prime}+rf_{0}^{\prime\prime}))\bigg{)}$
$\displaystyle-
b_{1}\bigg{(}rt^{2}\omega(-i+i\kappa+\omega)f_{0}^{2}\bar{f}_{0}-t\big{(}r^{2}(-2+\kappa+2i\omega)+6t^{2}b_{0}^{2}$
$\displaystyle\quad-2rtb_{0}b_{0}^{\prime}\big{)}f_{0}^{\prime}-2r(r^{2}-3t^{2}b_{0}^{2})\bar{f}_{0}f_{0}^{\prime
2}-r(r^{2}-3t^{2}b_{0}^{2})f_{0}^{\prime\prime}$
$\displaystyle\quad+f_{0}\Big{(}rt^{2}\omega(i-i\kappa+\omega)+\bar{f}_{0}\big{(}t(r^{2}(-2+\kappa-2i\omega)+6t^{2}b_{0}^{2}$
$\displaystyle\quad-2rtb_{0}b_{0}^{\prime})f_{0}^{\prime}+r(r^{2}-3t^{2}b_{0}^{2})f_{0}^{\prime\prime}\big{)}\Big{)}\bigg{)}\,,$

where

$L_{3}=rb_{0}(r^{2}-t^{2}b_{0}^{2})s_{0}\,.$ (33)

The above equations are scale-invariant. Hence, they turn to an ordinary
differential equation after making use of a change of coordinates to $(z,t)$.
Thus, the system of ordinary linear differential equations is given by

$\displaystyle b_{1}^{\prime}$
$\displaystyle=B_{1}(b_{1},f_{1},f_{1}^{\prime})\,,$ (34) $\displaystyle
f_{1}^{\prime\prime}$ $\displaystyle=F_{1}(b_{1},f_{1},f_{1}^{\prime})\,.$
(35)

Note that $b_{1}$ and $f_{1}$ are now linear functions that still depend non-
linearly on the unperturbed solution. We also have a quadratic dependence on
$\kappa$ as well. Note that these equations also do include the same
singularities as appeared in the unperturbed system of equations which means
that they are also singular for $z=0$ and $b^{2}(z)=z^{2}$. The modes are now
explored by finding the $\kappa$ values that are related to smooth solutions
of the perturbed equations (34), (35) where they need to satisfy the proper
boundary conditions as follows.

We now focus on the boundary conditions needed to solve (34) and (35). First
of all, at $z=0$, we re-scale the time coordinate so that $b_{1}(0)=0$. Also,
using the regularity condition for the axion-dilaton at $z=0$, we find that
$f_{1}^{\prime}(0)=0$ so that the freedom in $f$ is reduced to $f_{1}(0)$
which is an unknown complex parameter. On the other hand, at $z_{+}$ (we
recall that $z_{+}$ is defined by the equation $b(z_{+})=z_{+}$) all equations
and perturbations are regular so that all the second derivatives
$\partial_{r}^{2}f(t,r)$, $\partial_{r}\partial_{t}f(t,r)$,
$\partial_{t}^{2}f(t,r)$ should be finite as $z\rightarrow z_{+}$. Hence,
$f_{0}^{\prime\prime}(z)$ and $f_{1}^{\prime\prime}(z)$ are also finite as
$z\rightarrow z_{+}$. We introduce $\beta=b_{0}(z)-z$ and then expand
$f_{0}^{\prime\prime},f_{1}^{\prime\prime}$ near the singularity, as

$\displaystyle f_{0}^{\prime\prime}(\beta)$
$\displaystyle=\frac{1}{\beta}G(h_{0})+\mathcal{O}(1)\,,$ (36) $\displaystyle
f_{1}^{\prime\prime}(\beta)$
$\displaystyle=\frac{1}{\beta^{2}}\bar{G}(h_{0})+\frac{1}{\beta}H(h_{0},h_{1}|\kappa)+\mathcal{O}(1)\,,$
(37)

where we have defined the following equations:

$\displaystyle h_{0}=(b_{0}(z_{+}),f_{0}(z_{+}),f_{0}^{\prime}(z_{+})),\quad
h_{1}=(b_{1}(z_{+}),f_{1}(z_{+}),f_{1}^{\prime}(z_{+}))\;.$ (38)

The vanishing unperturbed complex constraint is given by $G(h_{0})=0$ at
$z_{+}$, and we checked that it implies $\bar{G}(h_{0})=0$ at $z_{+}$ as well,
which means that

$G(h_{0})=0\,\Rightarrow\,\bar{G}(h_{0})=0\,.$ (39)

Hence, we are left just with the complex-valued constraint

$H(h_{0},h_{1}|\kappa)=0\,,$ (40)

that is linear in $h_{1}$.

We now try to solve this constraint for $f_{1}^{\prime}(z_{+})$ in terms of
$f_{1}(z_{+})$, $b_{1}(z_{+})$, $\kappa$ and $h_{0}$. Hence, this condition
reduces the number of free parameters in the boundary conditions at $z_{+}$ to
a real number $b_{1}(z_{+})$ and a complex $f_{1}(z_{+})$. In conclusion, we
have six real unknowns, which are $\kappa$ and the five-component vector as

$X=(\real f_{1}(0),\,\imaginary f_{1}(0),\,\real f_{1}(z_{+}),\,\imaginary
f_{1}(z_{+}),\,b_{1}(z_{+}))\;,$ (41)

and the system of linear ODE’s (34) and (35) where the total real order is
indeed five. From the numerical procedure of [28] and given a set of boundary
conditions $X$, we integrate from $z=0$ to an intermediate point
$z_{\text{mid}}$, and similarly, we integrate backwards from $z_{+}$ to
$z_{\text{mid}}$. Finally, we collect the values of all functions
$(b_{1},\real f_{1},\imaginary f_{1},\real f_{1}^{\prime},\imaginary
f_{1}^{\prime})$ at $z_{\text{mid}}$ and encode the difference between the two
integrations in a difference function $D(\kappa;X)$. By definition,
$D(\kappa;X)$ is linear in $X$ thus it has a representation in matrix form

$D(\kappa;X)=A(\kappa)X\;,$ (42)

where $A(\kappa)$ is a $5\times 5$ real matrix depending on $\kappa$. Thus, we
only need to find the zeroes of $D(\kappa;X)$ and this can be achieved by
evaluating $\det A(\kappa)=0$. We carry out the root search for the
determinant as a function of $\kappa$ where the root with the biggest value
will be related to the Choptuik exponent through (20). Note that the perturbed
equations of motion are singular whenever the factor
$\left(\kappa-1-z\frac{b_{0}^{\prime}}{b_{0}}\right)$ in the denominator
vanishes, so that the numerical procedure fails at a particular point. We can
get an estimate for the values of $\kappa$ giving rise to this singular
behaviour, and we can also find the entire region that leads to numerical
failure. However, this apparent problem does not affect our evaluation of the
critical exponent because, in most cases, the most relevant mode $\kappa^{*}$
lies outside that particular failure region.

##  4 Statistical Approach

In this section, we detail the novel approach we have developed to estimate
the distribution of the critical exponent $\kappa$ in the elliptic 4d case of
the equations of motion with perturbations. In figure 1 we show the main steps
of the proposed pipeline.

Figure 1: Bayesian approach with artificial neural networks assisted
Metropolis-Hastings for the estimation of the distribution of $\kappa$. We
start with the estimation of the closed form for the unperturbed equations of
motion following polynomial regression. Then, a Bayesian artificial neural
networks assisted Metropolis-Hastings is followed, sampling $\kappa$ from a
prior distribution, and using an ANNs solver for the solution of the perturbed
equations of motion in each of the Metropolis-Hastings iterations.

Technically, we offer a new Bayesian approach with artificial neural networks
(ANNs) assisted Metropolis-Hastings for the estimation of the distribution of
the mentioned $\kappa$ that is related to the critical exponent $\gamma$. To
do so, we must treat $\kappa$ a random variable, so we use the Bayesian
strategy to find the posterior distribution of $\kappa$ based on the perturbed
DE system. We emphasize again that this is first time in the literature where
the critical exponent is treated as a random variable. As can be seen in
figure 1, our approach starts with the use a polynomial regression technique
to estimate the closed form of the unperturbed critical collapse functions
under analysis. To properly feed the polynomial regression technique, we
follow an ANN-based approach [32] to numerically find the unperturbed critical
collapse functions in the entire domain of the DE system. Our approach uses
the ANN estimates of the unperturbed critical collapse functions so that the
polynomial regression technique can be used for the estimation of the final
closed form. As our objective is the Bayesian analysis of the perturbed
equations, we use the estimated closed-form to incorporate the perturbations
related to $\kappa$. The last step consists of the implementation of ANN-
assisted Metropolis-Hastings to estimate the distribution of of $\kappa$ which
is related to the critical exponent. We detail each step in the following
subsections, providing all the details for the statistical methods employed.

###  4.1 Polynomial regression for the estimation of the closed form of the
unperturbed equations

Polynomial regression models are common tools in data science to model
multivariate data and find the relationship between variables. Let $({\bf
x},{\bf Z})$ be the input multivariate data, where ${\bf x}$ denotes the
response vector of sample size $n$ and ${\bf Z}$ denotes the design matrix of
size $(n\times p)$ with $p$ independent variables $({\bf z}_{1},\ldots,{\bf
z}_{p})^{\top}$. Given the training data of size $n$, the multivariate
regression

$\displaystyle x_{i}={\bf z}_{i}^{\top}{\bf\beta}+\epsilon_{i},$ (43)

where $\epsilon_{i}$ are independent and identically distributed (iid)
variables from the standard Gaussian distribution. The regression model (43)
essentially converts predicting the response differential equation variable
$x=g({\bf z})$ into estimating the unknown coefficients of the model. The
least squares (LS) method is a common approach to estimating the regression
model coefficients (43) minimizing the $l_{2}$ norm between the predicted
responses and observed responses as

$\displaystyle{\widehat{\beta}}=\arg\min_{\beta}||{\bf x}-{\bf
Z}{\bf\beta}||^{2}_{2}.$ (44)

One can easily show that the LS estimate of unknown coefficients can be
obtained by $\widehat{\bf\beta}_{LS}=({\bf Z}^{\top}{\bf Z})^{-1}{\bf Z}{\bf
x}$ as the solution to (44). Hence, the differential equation response $x(z)$
at every space-time point $z_{new}$ is predicted by ${\widehat{x}}_{new}={\bf
z}_{new}{\widehat{\beta}}_{LS}$ [46].

The polynomial regression model extends the properties of multivariate
regression by employing the higher-order terms of the independent variables to
better predict the nonlinear pattern of the response function $x=g(z)$. The
polynomial regression of order $J$ is given by

$\displaystyle x_{i}=\sum_{j=l}^{J}z_{i}^{j}\beta_{j}+\epsilon_{i},$ (45)

where $\beta=(\beta_{1},\ldots,\beta_{J})$ are the model coefficients that
should be estimated. One can easily represent (45) in matrix form as (43)
where the columns of the design matrix ${\bf Z}$ induce the explanatory
variable $z_{i},i=1,\ldots,n$ raised to various polynomial powers
$j=1,\ldots,J$. From the least squares method (44), one can predict the non-
linear responses of the differential equations at every point $z_{new}$ in the
domain of the equations of motion by ${\widehat{x}}_{new}={\bf
z}_{new}{\widehat{\beta}}_{LS}$.

Technically, we follow the described polynomial regression technique to first
obtain a closed form of the unperturbed equations of motions. As it has been
detailed, the polynomial regression needs pairs of corresponding inputs and
outputs. For generating them, for our problem, we follow an ANNs-based
approach [32] to numerically find the unperturbed critical collapse functions
in the entire domain of the DE system. In other words, we employ the technique
in [32] to generate the training data of the polynomial regressor.

###  4.2 ANNs assisted Metropolis-Hastings for the estimation of the
distribution of the critical exponent

We do have a closed form of the unperturbed equations of motion thanks to the
polynomial regression detailed in previous section. However, our objective
consists of finding the distribution of the parameter of $\kappa$ in the
perturbed equations of motion.

With the estimations of the closed form for the DE variables
$b_{0}(z),|f_{0}(z)|$ and $\arg(f_{0}(z))$ of the unperturbed equations of
motion, we proceed to incorporate these estimates and update the perturbed
equations of motion (34) and (35). We consider these perturbed equations of
motion, as our underlying DE system to estimate the distribution of the
critical exponent.

Consider a system of $H$ differential equations where ${\bf
x}(t)=\left(x_{1}(t),\ldots,x_{H}(t)\right)$ represents the multivariate
solutions evaluated at space-time $t$ to the system of differential equations
(DEs)

$\displaystyle\frac{d}{dt}x_{h}(t)=g_{h}({\bf x}(t)|\theta),$ (46)

where the DE variable $x_{h}(t)$ corresponds to $h$ DE and $\theta$ denotes
the collection of the unknown parameters of the DE. In our particular case
$\theta=\kappa$. Due to the high non-linearity of the black hole equations of
motion, the exact solution to the DE equations can be observed. Rather than
the exact DE solution, one may observe a perturbed version of the DE variable
with numerical measurement errors through numerical experiments. To take this
uncertainty into statistical methods, let ${\bf y}_{h}=(y_{h1},\ldots,y_{hn})$
denote the observed version of the DE variable $x_{h}(t)$ at $n$ space-time
points. Hence, let $y_{h,i}$ follow a Gaussian distribution with mean
$x_{h}(t_{i}|\theta)$ and variance $\sigma_{h}^{2}$. Given the Gaussian
distribution of the observed data, the uncertainty in the DE data can be
modelled by the likelihood function of $y_{i}$ as

$\displaystyle L(\theta|{\bf
y})=\prod_{h=1}^{H}\prod_{i=1}^{n}(1/\sigma_{i}^{2})\exp\left\\{-{(y_{hi}-x_{h}(t_{i}|\theta)}/{2\sigma_{h}^{2}}\right\\},$
(47)

The likelihood function translates the problem of finding the solution to a
system of DEs into a Gaussian process with unknown parameters. One can
estimate the unknown parameters of the DE system given the observed DE
variables by maximizing the likelihood function.

According to the non-linear nature of the DE variables, the likelihood
function (47) appears nonlinear with multiple optimums, indicating the
sensitivity of the maximum likelihood estimate of the DE variables.
Technically, in this work, we propose a Bayesian approach incorporating the
numerical measurement errors in estimating the critical exponent parameter in
the equations of motion (34),(35). This framework requires prior knowledge of
model parameters. This prior knowledge is incorporated into equations through
a prior probability density function. Let $\pi(\theta;{\bf\alpha})$ denote
this prior distribution, with ${\bf\alpha}$ representing the vector of hyper-
parameters. In the experiments, see Section 5, we employ different prior
distributions.

So, we are able to first generate $\kappa$ candidates from a proposal
distribution. Then, we propose to apply fully connected ANNs using the
$\kappa$ candidate and find the solution to the perturbed equations of motion
corresponding to the $\kappa$ candidate. This means that our model will
generate for each $\kappa$ candidate a solution using ANNs, assisting in an
iterative fashion the deployed Bayesian approach. We show now in detail how
this ANNs step is done.

Standard ANNs encompass multi-layer perceptrons, which transfer the data
between the layers through linear and non-linear functions [47]. In this
research, we focus on fully connected ANNs to deal with the non-linearity of
the elliptic perturbed equations of motion in the Einstein-axion-dilaton
system in 4d. Let $\mathcal{N}^{L}({\bf x},t,\psi)$ denote a neural networks
with $L$ layers which map the input dimensions $R^{in}$ to $R^{out}$. Let
${\bf W}^{l}$ and ${\bf b}^{l}$ denote the weight matrix and bias vector,
which regress the neurons in layer $l$ on $l-1$, respectively. Accordingly,
the response vector before activation in layer $l$ is given by

$\displaystyle\begin{array}[]{ll}\text{input layer:}&N^{0}({\bf
x},t,\psi)={\bf x},\\\ \text{hidden layers:}&N^{l}({\bf x},t,\psi)=\eta({\bf
W}^{l}N^{l-1}({\bf x},t,\psi)+{\bf b}^{l}),\\\ \forall 1\leq l\leq L-1&\\\
\text{output layer:}&N^{L}({\bf x},t,\psi)=\eta({\bf W}^{L}N^{L-1}({\bf
x},t,\psi)+{\bf b}^{L}),\\\ \end{array}$

where the output of layer $l$ is obtained after applying the activation
function, that is $\eta(\cdot)$.

The universal approximation theorem guarantees that the neural networks can
approximate any function, making the neural networks a state-of-the-art method
for solving nonlinear equations [47]. One must define a loss function
$\mathcal{L}$ to assess the performance of the ANNs and a back-propagation
algorithm to adjust the weight matrices and bias parameters of the ANNs based
on the gradient of the loss functions. The back-propagation algorithm
calculates the gradients of the ANNs with respect to the response variable to
ensure that the estimates of the ANN’s coefficients meet the requirements of
the loss function. This process translates solving the differential equations
of motion into an optimization problem, estimating the differential variables
based on minimizing the squared residues of the ANNs. Accordingly, for this
purpose, we introduce the $l_{2}$ loss function as

$\mathcal{L}\left(\mathcal{N}^{L}({\bf
x},t,\psi)\right)=\left(\mathcal{N}^{L}({\bf x},t,\psi)-g\right)^{2},$

where $\mathcal{N}^{L}({\bf x},t,\psi)$ represents the trained neural networks
approximating the underlying function.

Overall, at this point, our method disposes of a solution of the system of DEs
pertinent to the perturbed elliptic class of black hole equations of motion,
for a sampled particular value of $\kappa$. However, our objective is to find
the probability distribution of the critical exponent, so that it can be
considered, for the first time in the literature, as a random variable.
Therefore, for each of the solutions associated to sampled values of $\kappa$
we must compute the likelihood function (47) and find its acceptance
probability

Using information contained in the likelihood function (47) and the prior
distribution, one can find the posterior distribution of the unknown
parameters $\pi({\bf\theta}|{\bf y})$. The posterior distribution enables us
to capture the pattern of the parameters of the DE system using all the
uncertainty involved in calculating the DE variables. The posterior
distribution of $\theta$ is calculated by

$\displaystyle\pi({\bf\theta}|{\bf y})=\frac{L({\bf\theta}|{\bf
y})\pi({\bf\theta};{\bf\alpha})}{\int_{\theta}L({\bf\theta}|{\bf
y})\pi({\bf\theta};{\bf\alpha})}.$ (48)

Nevertheless, the posterior distribution (48) lacks a closed-form solution
because of the high non-linearity of the DE variables and the complex marginal
distributions of the DE variables. For this reason, one has to use iterative
methods such as Markov Chain Monte Carlo to find numerically the posterior
distribution of the parameters of the DE system [48].

Markov Chain Monte Carlo (MCMC) is a widely acceptable numerical technique
within the Bayesian framework to estimate the posterior distribution of the DE
system. The power of the MCMC approach relies on the fact that it transfers
the estimation task into the sampling process from the posterior distribution
of interest. Extensive research has been conducted in the literature about the
properties of the MCMC approach, including Metropolis-Hastings, Gibbs
sampling, and Sequential Monte Carlo [49] and [35]. In its nutshell, MCMC
employs various probabilistic techniques to generate a Markov chain of samples
approximating the target posterior distribution. When the chain is selected
long enough, the Markov chain will eventually reach a stationary state that
accurately fluctuates around the true form of the posterior distribution.

For our implementation, we follow the Metropolis-Hastings (MH) algorithm [34],
which is an established MCMC approach. The MH method calculates the
probability of the transition between the current state of the chain
$\theta^{(t)}$ and a candidate state $\theta^{*}$ simulated stochastically
from an independent proposal distribution $q(\theta)$. The MH method applies a
probabilistic step to accept or reject the proposed candidate. Let
$\theta^{(t)}$ represent the $t$-th state of the MCMC chain. Then we accept
the proposed state $\theta^{*}$ as the next state of the chain with
probability

$\displaystyle\min\left\\{1,\frac{\pi({\bf\theta}^{*}|{\bf
y})q({\bf\theta}^{(t)})}{\pi({\bf\theta}^{(t)}|{\bf
y})q({\bf\theta}^{*})}\right\\}.$ (49)

If the candidate ${\bf\theta}^{*}$ is accepted, then
${\bf\theta}^{(t+1)}={\bf\theta}^{*}$; otherwise,
${\bf\theta}^{(t+1)}={\bf\theta}^{(t)}$. The entire MCMC method is replicated
$N$ times to generate a sequence of $\\{{\bf\theta}^{(t)};t=1,\ldots,N\\}$
samples from the target posterior distribution $\pi({\bf\theta}|{\bf y})$.

Overall, we have detailed a novel model that allows for a Bayesian estimation
of the critical exponent on the elliptic black hole solution in 4d. In the
next section, we show all the numerical studies performed, and for the first
time in the literature, we are able to provide a probability distribution of
$\kappa$.

##  5 Numerical Studies

In this section, we plan to find the distribution of the critical exponent
$\kappa$ in the elliptic 4d equations of motion in the Bayesian framework. To
do so, the critical exponent is treated as a random variable, so we use the
Bayesian strategy to find the posterior distribution of $\kappa$ based on the
perturbed DE system. In this numerical study, we investigate the equations of
motion as the input DE system.

Hatefi et al. (2023) in [32] proposed ANNs to solve the unperturbed DE
variable. We applied the properties of polynomial regression and found the
closed-form polynomial regression estimates for the unperturbed critical
collapse functions. According to the fact that the goal of this research is to
estimate the critical exponent and that the critical exponent is only observed
in perturbed DE variables, we first followed [32] and used the ANNs based on
$\omega=1.176$ to numerically find the unperturbed critical collapse functions
in the entire domain of the DE system. To find the regression estimate of the
unperturbed critical collapse functions, we generated equally spaced points
${z_{1},\ldots,z_{n}}$ of size $n=1000$ from the domain of the DE system. We
obtained the ANN estimates of the unperturbed critical collapse functions
evaluated at each of the 1000 points. Once we obtained the ANN estimates, we
used 750 observations as training and 250 as testing data to estimate the
closed-form polynomial regression of the unperturbed critical collapse
functions, as described in Subsection 4.1. Hence, the closed-form polynomial
estimates of the unperturbed critical collapse functions are derived by

$\displaystyle\widehat{b_{0}(z)}$
$\displaystyle=1.005-0.187z+0.480z^{2}-0.004z^{3},$ (50)
$\displaystyle\widehat{|f_{0}(z)|}$
$\displaystyle=0.919-0.122z-0.028z^{2}-0.002z^{3},$ (51)
$\displaystyle\widehat{\arg(f_{0}(z))}$
$\displaystyle=-0.011+0.041z+0.047z^{2}-0.012z^{3}.$ (52)

Figure 2: The polynomial estimates of the unperturbed critical collapse
functions $b_{0}(z),\arg(f_{0}(z))$ and $|f_{0}(z)|$ of the equations of
motion based on polynomial orders of $l=1$ (black), $l=2$ (red) and $l=3$
(blue). The dots show the true functional form evaluated at the validation
points.

Figure 2 shows the training and test performance of the polynomial estimates
for the critical collapse functions for polynomial orders $l=1,2,3$. One can
easily see that the polynomial of order $l=3$ can capture very well the non-
linearity of the unperturbed critical collapse functions. Once we found the
unperturbed DE variables $b_{0}(z),|f_{0}(z)|$ and $\arg(f_{0}(z))$, we
incorporated these estimates and updated the perturbed equations of motion. In
the next step, we consider this perturbed equations of motion, as our
underlying DE system to estimate the critical exponent.

We proposed an artificial neural network-assisted Metropolis-Hastings to
obtain the posterior distribution of $\kappa$. We first generated $\kappa$
candidates from a proposal distribution to do so. Here, we considered two
proposal distributions, including the Uniform distribution between (3,4) and
the Gaussian distribution with a mean of 3.5 and variance of 1 to deal with
the singularity issue, and the range included the established solution
available in the literature [28].

In the next step, we apply fully connected neural networks using the $\kappa$
candidate and find the solution to the perturbed equations of motion
corresponding to the $\kappa$ candidate. We implemented the fully connected
ANNs using the Python package neurodiffeq [50]. The neural networks used
contain 4 hidden layers, each with 16 neurons. The training was done for 50
epochs. We then computed the likelihood function (47) and found the acceptance
probability of the Metropolis-Hasting algorithm under non-informative uniform
distribution. Due to the fact that the range of DE variables is dramatically
different, we assign $\sigma_{b_{1}}=7$ and $\sigma_{g}=2$ in the likelihood
function of $\kappa$ given the perturbed DE variables.

It is shown in the literature that the perturbed $b_{1}(z)$ DE variable is
linear with respect to $\kappa$; however, the perturbed $f_{1m}(z)=g_{m}(z)$
contains the higher orders of $\kappa$ [28]. To investigate the effects of the
different DE variables on the MCMC chains and the Bayesian distribution of
$\kappa$, we evaluated the performance of the stochastic accept-reject method
based on three scenarios. These scenarios include the accept-reject step based
on only the observed perturbed DE variable of $b_{1}(z)$, accept-reject based
on only the observed perturbed DE variable $g_{m}(z)$ and accept-reject based
on both observed perturbed DE variables $g_{m}(z)$ and $b_{1}(z)$. Finally,
the entire neural network-assisted Metropolis-Hasting approach was carried out
for 2000 chains, where in each iteration, the perturbed equations of motion
were computed using artificial neural networks. We discard the first 10% of
the history of the MCMC chain to wash out the effect of the initialization
step on the performance of the MCMC chain.

Figures 3-8 show the posterior distributions and trace plots of the $\kappa$
parameter based on 1800 MCMC chains (after discarding the burn-in period)
using our artificial neural networks-assisted Metropolis-Hastings approach.
The results are based on three accept-reject approaches and two proposal
distributions, including uniform and Gaussian distributions. From Figures 3-8,
we observe that the distribution’s posterior mean and posterior mode, as two
Bayesian estimates of the $\kappa$ range between 3.2 and 3.8, which verifies
the findings in the literature that $\kappa^{*}\approx 3.7858$ for the
equations of motion in 4d elliptic class [28]. From the trace plot, we easily
observed that the MCMC chain easily searches the domain of the parameter of
interest and supports the solution between 3.2 and 3.8.

As one needs the ANNs in each iteration of the MCMC chain to evaluate the
accept-reject steps, we show the results of the loss function of ANNs in
Figures 9-14 (in the Appendix). To show the convergence of the ANNs, we
computed the difference between training and validation loss functions in the
last epoch of the ANNs over 1800 MCMC samples. In each figure, we also show
the entire trajectory of the training and validation loss functions over all
the epochs for a randomly selected MCMC sample. From Figures 9-14, we clearly
observe the loss differences are almost very small of magnitude $10^{-6}$
fluctuating around zero. This highlights the convergence of the ANNs in
estimating the critical exponent in a Bayesian framework.

The unique solution for the elliptic case in four dimensions was achieved in
[28]. Indeed, the behaviour of $detA(k)$ near the last crossing of the
horizontal axis was estimated to be

$\kappa^{*}_{4E}\approx 3.7858\,,$ (53)

which gives rise to a Choptuik exponent as $\gamma_{4E}\approx 0.2641$. In
this paper, we have used various statistical methods to actually explore the
entire range of the critical exponent rather than exploring the last crossing.
Surprisingly, our results

$3.2<\kappa^{*}_{4E}<3.8$ (54)

is in perfect agreement with the literature, see [41]. Hence, we can actually
take this very non-trivial match as evidence that shows all our numerical sets
up are reliable methods that can be extended to any dimension and any matter
content.

Figure 3: The posterior distribution and trace plot of the $\kappa$ based on
MCMC samples from the Metropolis-Hastings algorithm using only the perturbed
DE variable $b_{1}(z)$ to accept or reject the transitions under the Uniform
proposal distribution. Figure 4: The posterior distribution and trace plot of
the $\kappa$ based on MCMC samples from the Metropolis-Hastings algorithm
using only the perturbed DE variable $f_{1m}(z)=g_{m}(z)$ to accept or reject
the transitions under the Uniform proposal distribution. Figure 5: The
posterior distribution and trace plot of the $\kappa$ based on MCMC samples
from the Metropolis-Hastings algorithm using both the perturbed DE variables
$b_{1}(z)$ and $g_{m}(z)$ to accept or reject the transitions under the
Uniform proposal distribution.

##  6 Conclusion

In this paper, we advance the current literature on critical exponents for
black hole solutions by examining all solutions within the domains of the
linear perturbation equations of motion. Specifically, we investigate the
quantum perturbation theory for a four-dimensional Einstein-axion-dilaton
system under an elliptic class of $\text{SL}(2,\mathbb{R})$ transformations.
Utilizing quantum perturbation theory, we propose artificial neural network-
assisted Metropolis-Hastings algorithms for Bayesian estimation of the
critical exponent of the elliptic black hole solution in 4d.

Unlike conventional methods, we investigated the range of possible values for
critical exponents using quantum perturbation theory. We conducted a thorough
analysis not only of the perturbed differential equation (DE) variables
$b_{1}(z)$ and $f_{1m}(z)=g_{m}(z)$ but also of the DE variables for both
$b_{1}(z)$ and $g_{m}(z)$ equations simultaneously. This comprehensive
approach allowed us to assess their effects on stochastic accept-reject
transitions within the posterior distributions of the critical exponent. Our
innovative probabilistic method stands out from existing techniques by
offering the established solution while simultaneously exploring the entire
spectrum of physically distinguishable critical exponents, accounting for
potential numerical measurement errors. This advancement provides a more
robust and inclusive understanding of critical exponents, highlighting the
limitations of previous approaches and paving the way for more accurate
predictions in complex systems.

Our new methods introduce innovative approaches for exploring the potential
range of allowed values for the critical exponent, applicable to all different
conjugacy classes of the $\mathrm{SL}(2,\mathbb{R})$ transformation. Moreover,
it is important to emphasize that these methods can be extended to various
types of matter content. This is a direction we are keen to pursue in future
research.

It is interesting to note that, besides the critical exponents’ dependence on
matter content, dimension, and various solutions of self-critical collapse
[25], these exponents can span an entire range of values rather than being
confined to a single localized value. Therefore, we conclude that the
conjecture regarding the universality of the Choptuik exponent does not hold.
However, some universal behaviours might be embedded in the combinations of
critical exponents and other parameters of the given theory that our current
analysis has not accounted for. Our new findings clearly indicate that the
standard expectations of statistical mechanics do not transfer to the context
of critical gravitational collapse.

##  7 Acknowledgments

E. Hatefi would like to thank Pablo Diaz, E. Hirschmann, Philip Siegmann, L.
Alvarez-Gaume, and A. Sagnotti for their useful discussions and support. E.
Hatefi is supported by the María Zambrano Grant of the Ministry of
Universities of Spain. Armin Hatefi acknowledges the support from the Natural
Sciences and Engineering Research Council of Canada (NSERC).

## References

  * [1] M.W. Choptuik, Universality and Scaling in Gravitational Collapse of a Massless Scalar Field, Phys. Rev. Lett. 70, 9 (1993).
  * [2] D. Christodoulou, “The Problem of a Self-gravitating Scalar Field,” Commun. Math. Phys. 105 (1986) 337; “Global Existence of Generalized Solutions of the Spherically Symmetric Einstein Scalar Equations in the Large,” Commun. Math. Phys. 106 (1986) 587;“The Structure and Uniqueness of Generalized Solutions of the Spherically Symmetric Einstein Scalar Equations,” Commun. Math. Phys. 109 (1987) 591.
  * [3] R. S. Hamade and J. M. Stewart, “The Spherically symmetric collapse of a massless scalar field,” Class. Quant. Grav. 13 (1996) 497 [arXiv:gr-qc/9506044].
  * [4] C. Gundlach, “Critical phenomena in gravitational collapse,” Phys. Rept. 376 (2003) 339 [gr-qc/0210101].
  * [5] T. Koike, T. Hara, and S. Adachi, “Critical Behavior in Gravitational Collapse of Radiation Fluid: A Renormalization Group (Linear Perturbation) Analysis,” Phys. Rev. Lett. 74 (1995) 5170 [gr-qc/9503007].
  * [6] L. Alvarez-Gaume, C. Gomez and M. A. Vazquez-Mozo, “Scaling Phenomena in Gravity from QCD,” Phys. Lett. B 649 (2007) 478 [hep-th/0611312].
  * [7] M. Birukou, V. Husain, G. Kunstatter, E. Vaz and M. Olivier, “Scalar field collapse in any dimension,” Phys. Rev. D 65 (2002) 104036 [gr-qc/0201026].
  * [8] V. Husain, G. Kunstatter, B. Preston and M. Birukou, “Anti-de Sitter gravitational collapse,” Class. Quant. Grav. 20 (2003) L23 [gr-qc/0210011].
  * [9] E. Sorkin and Y. Oren, “On Choptuik’s scaling in higher dimensions,” Phys. Rev. D 71, 124005 (2005) [arXiv:hep-th/0502034].
  * [10] J. Bland, B. Preston, M. Becker, G. Kunstatter and V. Husain, “Dimension-dependence of the critical exponent in spherically symmetric gravitational collapse,” Class. Quant. Grav. 22 (2005) 5355 [gr-qc/0507088].
  * [11] E. W. Hirschmann and D. M. Eardley, “Universal scaling and echoing in gravitational collapse of a complex scalar field,” Phys. Rev. D 51 (1995) 4198 [gr-qc/9412066].
  * [12] J. V. Rocha and M. Tomašević, “Self-similarity in Einstein-Maxwell-dilaton theories and critical collapse,” Phys. Rev. D 98 (2018) no.10, 104063 [arXiv:1810.04907 [gr-qc]].
  * [13] L. Alvarez-Gaume, C. Gomez, A. Sabio Vera, A. Tavanfar and M. A. Vazquez-Mozo,“Critical gravitational collapse: towards a holographic understanding of the Regge region,” Nucl. Phys. B 806 (2009) 327 [arXiv:0804.1464 [hep-th]].
  * [14] C. R. Evans and J. S. Coleman, “Observation of critical phenomena and selfsimilarity in the gravitational collapse of radiation fluid,” Phys. Rev. Lett. 72 (1994) 1782 [gr-qc/9402041].
  * [15] D. Maison, “Non-Universality of Critical Behaviour in Spherically Symmetric Gravitational Collapse,” Phys. Lett. B 366 (1996) 82 [gr-qc/9504008].
  * [16] A. Strominger and L. Thorlacius, “Universality and scaling at the onset of quantum black hole formation,” Phys. Rev. Lett. 72 (1994) 1584 [hep-th/9312017].
  * [17] E. W. Hirschmann and D. M. Eardley, ‘Critical exponents and stability at the black hole threshold for a complex scalar field,” Phys. Rev. D 52 (1995) 5850 [gr-qc/9506078].
  * [18] A. M. Abrahams and C. R. Evans, “Critical behavior and scaling in vacuum axisymmetric gravitational collapse,” Phys. Rev. Lett. 70 (1993) 2980.
  * [19] L. Alvarez-Gaume, C. Gomez, A. Sabio Vera, A. Tavanfar and M. A. Vazquez-Mozo, “Critical formation of trapped surfaces in the collision of gravitational shock waves,” JHEP 0902, 009 (2009) [arXiv:0811.3969 [hep-th]].
  * [20] E. W. Hirschmann and D. M. Eardley, “Criticality and bifurcation in the gravitational collapse of a selfcoupled scalar field,” Phys. Rev. D 56 (1997) 4696 [gr-qc/9511052].
  * [21] J. M. Maldacena, “The Large N limit of superconformal field theories and supergravity,”Int. J. Theor. Phys. 38 (1999), 1113-1133, Adv. Theor. Math. Phys. 2, arXiv:hep-th/9711200, E. Witten, “Anti-de Sitter space and holography,” Adv. Theor. Math. Phys. 2 (1998), 253-291,hep-th/9802150, S. Gubser, I. R. Klebanov and A. M. Polyakov, “Gauge theory correlators from noncritical string theory,” Phys. Lett. B 428 (1998), 105-114, hep-th/9802109.
  * [22] D. Birmingham, “Choptuik scaling and quasinormal modes in the AdS / CFT correspondence,” Phys. Rev. D 64 (2001), 064024 [arXiv:hep-th/0101194 [hep-th]].
  * [23] E. Hatefi, A. Nurmagambetov and I. Park, “ADM reduction of IIB on $\mathcal{H}^{p,q}$ to dS braneworld,” JHEP 04 (2013), 170, arXiv:1210.3825 , “$N^{3}$ entropy of $M5$ branes from dielectric effect,” NPB866 (2013), 58-71, arXiv:1204.2711, S. de Alwis, R. Gupta, E. Hatefi and F. Quevedo, “Stability, Tunneling and Flux Changing de Sitter Transitions in the Large Volume String Scenario,” JHEP 11 (2013), 179, arXiv:1308.1222; E. Hatefi and I. Y. Park, Phys. Rev. D 85 (2012), 125039 [arXiv:1203.5553 [hep-th]]; E. Hatefi, JHEP 05 (2010), 080,[arXiv:1003.0314 [hep-th]].
  * [24] E. Hatefi, Phys. Rev. D 86 (2012), 046003, [arXiv:1203.1329 [hep-th]]; E. Hatefi and I. Y. Park, Nucl. Phys. B 864 (2012), 640-663 [arXiv:1205.5079 [hep-th]]; E. Hatefi, JCAP 09 (2013), 011, [arXiv:1211.5538 [hep-th]]; E. Hatefi, JHEP 07 (2013), 002, [arXiv:1304.3711 [hep-th]]; E. Hatefi, JHEP 04 (2013), 070, [arXiv:1211.2413 [hep-th]]; E. Hatefi, JHEP 11 (2013), 204, [arXiv:1307.3520 [hep-th]].
  * [25] R. Antonelli and E. Hatefi, “On self-similar axion-dilaton configurations,” , JHEP 03 (2020), 074 [arXiv:1912.00078 [hep-th]].
  * [26] L. Álvarez-Gaumé and E. Hatefi, “Critical Collapse in the Axion-Dilaton System in Diverse Dimensions,” Class. Quant. Grav. 29 (2012) 025006 [arXiv:1108.0078 [gr-qc]].
  * [27] L. Álvarez-Gaumé and E. Hatefi, “More On Critical Collapse of Axion-Dilaton System in Dimension Four,” JCAP 1310 (2013) 037 [arXiv:1307.1378 [gr-qc]].
  * [28] R. Antonelli and E. Hatefi, “On Critical Exponents for Self-Similar Collapse,” JHEP 03 (2020), 180 [arXiv:1912.06103 [hep-th]].
  * [29] E. Hatefi and A. Hatefi, “Estimation of Critical Collapse Solutions to Black Holes with Nonlinear Statistical Models,” Mathematics 10 (2022) no.23, 4537 [arXiv:2110.07153 [gr-qc]].
  * [30] E. Hatefi and A. Hatefi, “Nonlinear statistical spline smoothers for critical spherical black hole solutions in 4-dimension,” Annals Phys. 446 (2022), 169112 [arXiv:2201.00949 [gr-qc]].
  * [31] E. Hatefi, A. Hatefi and R. J. López-Sastre, “Analysis of black hole solutions in parabolic class using neural networks,” Eur. Phys. J. C 83 (2023) no.07, 623 [arXiv:2302.04619 [gr-qc]].
  * [32] A. Hatefi, E. Hatefi and R. J. López-Sastre, “Modeling the complexity of elliptic black hole solution in 4D using Hamiltonian Monte Carlo with stacked neural networks,” JHEP 10 (2023), 034 [arXiv:2307.14515 [gr-qc]].
  * [33] A. Hatefi and E. Hatefi, ‘Sequential Monte Carlo with cross-validated neural networks for complexity of hyperbolic black hole solutions in 4D,” Eur. Phys. J. C 83 (2023) no.11, 1083 [arXiv:2308.07907 [hep-th]].
  * [34] Robert, C. P., Casella, G.,“ Monte Carlo statistical methods“ (Vol. 2). New York: Springer, (1999).
  * [35] Murphy, K. P., “Probabilistic machine learning: Advanced topics“. MIT Press, 2023.
  * [36] A. Sen, “Strong - weak coupling duality in four-dimensional string theory,” Int. J. Mod. Phys. A 9 (1994) 3707 [hep-th/9402002].
  * [37] J. H. Schwarz, “Evidence for nonperturbative string symmetries,” Lett. Math. Phys. 34 (1995) 309 [hep-th/9411178].
  * [38] M.B. Green, J.H. Schwarz and E. Witten, 1987 Superstring Theory Vols I,II, Cambridge University Press,
  * [39] J. Polchinski, 1998 String Theory, Vols I,II, Cambridge University Press
  * [40] A. Font, L. E. Ibanez, D. Lust and F. Quevedo, “Strong - weak coupling duality and nonperturbative effects in string theory,” Phys. Lett. B 249 (1990) 35.
  * [41] D. M. Eardley, E. W. Hirschmann and J. H. Horne, “S duality at the black hole threshold in gravitational collapse,” Phys. Rev. D 52 (1995) 5397 [arXiv:gr-qc/9505041].
  * [42] E. Hatefi and E. Vanzan, “On Higher Dimensional Self-Similar Axion-Dilaton Solutions,” Eur. Phys. J. C 80 (2020), 10 [arXiv:2005.11646 [hep-th]].
  * [43] R. S. Hamade, J. H. Horne and J. M. Stewart, “Continuous Self-Similarity and $S$-Duality,” Class. Quant. Grav. 13 (1996) 2241 [arXiv:gr-qc/9511024].
  * [44] E. Hatefi and A. Kuntz, “On perturbation theory and critical exponents for self-similar systems,” Eur. Phys. J. C 81 (2021) no.15 [arXiv:2010.11603 [hep-th]].
  * [45] A. Ghodsi and E. Hatefi, “Extremal rotating solutions in Horava Gravity,” Phys. Rev. D 81 (2010) 044016 [arXiv:0906.1237 [hep-th]].
  * [46] F. Harrell Jr, “ Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis,“ Springer (2015).
  * [47] Goodfellow, I., Bengio, Y., & Courville, A. “Deep learning”, MIT press. (2016)
  * [48] Girolami, M., “Bayesian inference for differential equations. “ Theoretical Computer Science, 408(1), 4-16, (2008).
  * [49] Bishop, C. M., and Nasrabadi, N. M., “Pattern recognition and machine learning“ (Vol. 4, No. 4, p. 738). New York: Springer, (2006).
  * [50] F. Chen, D. Sondak, P. Protopapas, P.Mattheakis, M. Liu, S.Agarwal, D. & Di Giovanni, M., “NeuroDiffEq: A Python package for solving differential equations with neural networks,” Journal Of Open Source Software. 5, 1931 (2020)

## Appendix

Figure 6: The posterior distribution and trace plot of the $\kappa$ based on
MCMC samples from the Metropolis-Hastings algorithm using only the perturbed
DE variable $b_{1}(z)$ to accept or reject the transitions under the Gaussian
proposal distribution. Figure 7: The posterior distribution and trace plot of
the $\kappa$ based on MCMC samples from the Metropolis-Hastings algorithm
using only the perturbed DE variable $f_{1m}(z)=g_{m}(z)$ to accept or reject
the transitions under the Gaussian proposal distribution. Figure 8: The
posterior distribution and trace plot of the $\kappa$ based on MCMC samples
from the Metropolis-Hastings algorithm using both the perturbed DE variables
$b_{1}(z)$ and $g_{m}(z)$ to accept or reject the transitions under the
Gaussian proposal distribution. Figure 9:  The difference between the train
and validation losses in the last epochs over MCMC samples (A) and the
histogram of the loss differences (B). The train and validation losses for one
random MCMC sampling (C). The MCMC samples are from the Metropolis-Hastings
algorithm, using only the perturbed DE variable $b_{1}(z)$ to accept or reject
the transitions under the Uniform proposal distribution. Figure 10:  The
difference between the train and validation losses in the last epochs over
MCMC samples (A) and the histogram of the loss differences (B). The train and
validation losses for one random MCMC sampling (C). The MCMC samples are from
the Metropolis-Hastings algorithm, using only the perturbed DE variable
$g_{m}(z)$ to accept or reject the transitions under the Uniform proposal
distribution. Figure 11:  The difference between the train and validation
losses in the last epochs over 1800 MCMC samples (A) and the histogram of the
loss differences (B). The train and validation losses for one random MCMC
sampling (C). The MCMC samples are from the Metropolis-Hastings algorithm,
using both perturbed DE variables $b_{1}(z)$ and $g_{m}(z)$ to accept or
reject the transitions under the Uniform proposal distribution. Figure 12:
The difference between the train and validation losses in the last epochs over
MCMC samples (A) and the histogram of the loss differences (B). The train and
validation losses for one random MCMC sampling (C). The MCMC samples are from
the Metropolis-Hastings algorithm, using only perturbed DE variables
$b_{1}(z)$ to accept or reject the transitions under the Gaussian proposal
distribution. Figure 13: The difference between the train and validation
losses in the last epochs over MCMC samples (A) and the histogram of the loss
differences (B). The train and validation losses for one random MCMC sampling
(C). The MCMC samples are from the Metropolis-Hastings algorithm, using only
the perturbed DE variable $g_{m}(z)$ to accept or reject the transitions under
the Gaussian proposal distribution. Figure 14:  The difference between the
train and validation losses in the last epochs over MCMC samples (A) and the
histogram of the loss differences (B). The train and validation losses for one
random MCMC sampling (C). The MCMC samples are from the Metropolis-Hastings
algorithm, using both perturbed DE variables $b_{1}(z)$ and $g_{m}(z)$ to
accept or reject the transitions under the Gaussian proposal distribution.