This paper addresses the $\epsilon$-close parameter tuning problem for Bayesian
networks (BNs): find a minimal $\epsilon$-close amendment of probability entries
in a given set of (rows in) conditional probability tables that make a
given quantitative constraint on the BN valid. Based on the
state-of-the-art “region verification” techniques for parametric Markov
chains, we propose an algorithm whose capabilities go
beyond any existing techniques. Our experiments show that $\epsilon$-close tuning
of large BN benchmarks with up to eight parameters is feasible. In
particular, by allowing (i) varied parameters in multiple CPTs and (ii)
inter-CPT parameter dependencies, we treat subclasses of parametric BNs
that have received scant attention so far.

§ INTRODUCTION

Bayesian networks.
Bayesian networks (BNs) [21] are probabilistic graphical models that enable succinct knowledge representation and facilitate probabilistic reasoning [12].

Parametric Bayesian networks (pBNs) [5]
extend BNs by allowing polynomials in conditional probability tables (CPTs) rather than constants.

* Expert's knowledge to build and modify the network <cit.>

* The cost of the modification and the feasibility of the modification: it is favorable that the amendment is small and the parameters acquire new probabilities close to their original values.

* Uncertainty <cit.> and fluctuation of parameters <cit.>

Parameter synthesis on Markov models. Parameter synthesis is to find the right values for the unknown parameters with respect to a given constraint. Various synthesis techniques have been developed for parametric Markov chains (pMCs) ranging over e.g., the gradient-based methods [14], convex optimization [10, 11], and region verification [23].
Recently, Salmani and Katoen [[26]] have proposed a translation from pBNs to pMCs that facilitates using pMC algorithms to analyze pBNs. Proceeding from this study, we tackle a different problem [18] for Bayesian networks.

Minimal-change parameter tuning. Given a Bayesian network $B$, a hypothesis $H$, evidence $E$, $\lambda \in [0, 1]$, and a constraint of the form $\Pr(H | E) \leq \lambda$ (or $\geq  \lambda$),
[backgroundcolor=blue!5, nobreak=true, skipabove=4pt, skipbelow=0pt]
what is a minimal change—with respect to a given measure of distance—in the probability values of (a subset of) CPT rows, such that the constraint holds?

We illustrate the problem with an example of testing COVID-19 that is adopted from several medical studies [2, 20, 13]. The outcome of PCR tests only depends on whether the person is infected, while the antigen tests are less likely to correctly identify COVID-19 if the infection is asymptomatic (or pre-symptomatic). Figure <ref> depicts a Bayesian network that models such probabilistic dependencies.
In the original network, the probability of no COVID-19 given that both tests are positive, is $0.011089$.
Assume that in an application domain, the result of such a query is now required not to exceed $0.009$: the aim is to make this constraint hold while imposing the least change in the original network with respect to a distance measure.
This is an instance of the minimal-change parameter tuning problem. We consider the $\epsilon-$bounded variant of the problem: for a subset of modifiable CPT rows, are there new values within the distance of $\epsilon$ from the original probability values that make the constraint hold?

Main contributions.  Based on existing region verification and region partitioning techniques for pMCs,

* we propose a practical algorithm for $\epsilon$-bounded tuning.
More precisely, we find instantiations that (i) satisfy the constraint (if the constraint is satisfiable) and (ii) are $\epsilon$-close and (iii) lean towards the minimum distance instantiation depending on a coverage factor $0 \leq \eta \leq 1$.

We propose two region expansion schemes to realize $\epsilon$-closeness of the results both for Euclidean distance and for CD distance [7].

* Contrary to the existing techniques that restrict to hyperplane solution spaces, we handle pBNs with multiple parameters in multiple distributions and multiple CPTs.

Our experiments on our prototypical implementation [<https://github.com/baharslmn/pbn-epsilon-tuning>] indicate that $\epsilon$-bounded tuning of up to 8 parameters for large networks with 100 variables is feasible.

Paper organization. Section 2 includes the basic notations and Sec. 3 the preliminaries on parametric Bayesian networks. Section 4 introduces parametric Markov chains and the region verification techniques thereof. Section 5 details our main contributions and Sec. 6 our experimental results. Section 7 concludes the paper with an overview of the related studies that were not mentioned before.

Contributions. Our contributions are as follows.

* We proposed an algorithm for minimal change parameter tuning on Bayesian networks that scales to the simultaneous variation of up to $8$ parameters in large BN benchmarks with thousands “hyper-parameters".

* We empirically indicated the applicability of our approach to subclasses of the problem that yield sensitivity functions that are not necessarily multi-linear in terms of the varied parameters. Our completes the conventional setting for Bayesian networks that mostly impose restrictions on the inter-CPT parameter dependencies.

* We empirically evaluated our approach against the existing tuning tools for Bayesian networks.

* We performed a detailed study on the efficiency of our algorithm for several benchmarks with multiple settings.

Antigen: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8683002/> cited, <https://www.cochranelibrary.com/cdsr/doi/10.1002/14651858.CD013705.pub2/full> cited

PCR: <https://www.health.govt.nz/covid-19-novel-coronavirus/covid-19-health-advice-public/covid-19-testing/covid-19-test-results-and-their-accuracy>
Symptoms: <https://www.ijidonline.com/article/S1201-9712(20)30139-9/pdf> cited

§ BACKGROUND

Let $V$ be a set of $m$ random variables $v_1, \ldots, v_m$ and $D_{v_i}$ the domain of variable $v_i$.
For $A \subseteq V$,

$\mbox{\it Eval}(A)$ denotes the set of joint configurations for the variables in $A$.

Let $X \,{=}\, \{x_1, \ldots, x_n\}$ be a set of $n$ real-valued parameters.
A parameter instantiation is a function $u\,{:}\, X \,{\rightarrow}\,\mathbb{R}$ that maps each parameter to a value.

All parameters are bounded; i.e., $lb_{x_i} \,{\leq}\, u(x_i) \,{\leq}\, ub_{x_i}$ for $x_i$.
Let $I_i \,{=}\, [lb_{x_i}, ub_{x_i}]$.

The parameter space $\mathcal{U} \subseteq \mathbb{R}^{n}$ of $X$ is the set of all possible values of $X$, i.e., the hyper-rectangle spanned by the intervals $I_i$ for all $i$.

Polynomials $f$ over $X$ are functions $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$ where $f(u)$ is obtained by replacing each occurrence of $x_i$ in $f$ by $u(x_i)$; e.g., for $f = 2x_1^2{+}x_2$, $u(x_1) = 3$ and $u(x_2) = 2$, $f(u) = 20$.
Let $f[u]$ denote such substitution.

Parametric distributions. Let $\Distr(D)$ denote the set of probability distributions over $D$. Let $\mathbb{Q}[X]$ be the set of multivariate polynomials with rational coefficients over $X$. A parametric probability distribution is the function $\mu: D \rightarrow \mathbb{Q}[X]$ with $\sum_{d \in D} \mu(d) \,{=}\, 1$.  Let $\mbox{\it pDistr}(D)$ denote the set of parametric probability distributions over $D$ with parameters in $X$. Instantiation $u$ is well-formed for $\mu \,{\in}\, \mbox{\it pDistr}(D)$ iff $0 \leq \mu(d)[u] \leq 1$ and $\Sigma_{d \in D} \, \mu(d)[u] = 1$.

Co-variation scheme.
A co-variation scheme $cov: \Distr(D) \times D \to \pDistr(D)$ maps a probability distribution $\mu$ to a parametric distribution $\mu'$ based on a given  $d \in D$.

Distance measure.
The function $d: \mathcal{U} \times \mathcal{U}\rightarrow \mathbb{R}^{\geq 0}$ is a distance measure if for all $u, u', u'' \in \mathcal{U}$, it satisfies: (I) positiveness: $d(u,u') \geq 0$, (II) symmetry: $d(u,u') = d(u',u)$, and (III) triangle inequality: $d(u,u'') \leq d(u,u') + d(u',u'')$.

§ PARAMETRIC MARKOV CHAINS

Parametric Markov chains are an extension of Markov chains (MCs) that allow parametric polynomials as transition labels:

A parametric Markov chain (pMC) $\mathcal{M}$ is a tuple $(S, s_{I}, X, \mathcal{P})$ where $S$ is a finite set of states, $s_I \in S$ is the initial state, $X$ is a finite set of real-valued parameters, and $\mathcal{P} : S  \rightarrow \pDistr(S)$ is a transition function over the states.

Figures <ref> and <ref> are examples of pMCs.

The parametric transition probability of going from state $s$ to $t$ is given by $\mathcal{P}(s, t)$.
Intuitively, instead of a concrete probability, a transition is equipped with a polynomial over the parameters $X$.
A pMC with $X = \emptyset$ and $\mathcal{P} : S \rightarrow \mbox{\it Distr}(S)$ is a Markov chain (MC).
Applying the instantiation $u: X \rightarrow \mathbb{R}$ to the pMC $\mathcal{M}$ yields $\mathcal{M}[u]$ by replacing each transition function $f \in \mathbb{Q}[X]$ in $\mathcal{M}$ by $f[u]$.
An instantiation $u$ is well-defined (for $\mathcal{M}$) if the transition function induced by $u$ yields probability distributions, i.e., if $\mathcal{P}(s, .)[u] \in \mbox{\it Distr}(S)$, for each $s \in S$.
For the well-defined instantiation $u$,  $\mathcal{M}[u]$ is an MC.  The well-defined instantiation $u$ is graph-preserving for $\M$ if the topology is preserved, i.e., $\mathcal{P}(s,s') \neq 0$ implies $\mathcal{P}(s,s')[u] \neq 0$ for all states $s$ and $s'$. A region $R$ is well-defined (graph-preserving) if $u$ is well-defined (graph-preserving) for all $u \in R$. We assume in the sequel that all instantiations are graph-preserving.

For the MC $M = (S, s_I, \mathcal{P})$ and the state $s \in S$, let $succ(s)$ denote the set of all successors of $s$ in $M$, i.e., $s' \in S$ with $\mathcal{P}(s,s') > 0$. Let $\mathbb{Q}_{\neq 0}[X]$ denote the set of multivariate polynomials with rational coefficients over $X$ excluding the constant polynomial $0$.  Let $S_{\modif} \subseteq S$ be the set of modifiable states in $M$. In what follows, let $s_m \in S_{\modif}$.

The pMC $\M_p = (\boldsymbol{S}, \boldsymbol{s_I}, X, \boldsymbol{\mathcal{P}})$ is a parametrization of $\M$ over the parameters set $X$ with respect to $S_{\modif}$  iff $\boldsymbol{S} = S$, $\boldsymbol{s_I} = s_I$, and for $s, s' \in S$,

\begin{align*}
\mathcal{P}(s_m, s') =  f  \in \mathbb{Q}[X] \text{ for } s' \in

\end{align*}

\begin{align*}
&\mathcal{P}(s, s') = 0 \iff \boldsymbol{\mathcal{P}}(s,s') = 0
\quad \text{ and } \\
& \mathcal{P}(s, s') \neq 0 \implies \boldsymbol{\mathcal{P}}(s,s') = f  \in \mathbb{Q}_{\neq 0}[X].
\end{align*}

The parametrization thus preserves the topology of the MC in the obtained pMC and does not introduce new edges between the states. The pMC $\M_p$ is a graph-preserving parametrization of MC $\M$ iff it is a parametrization of $\M$

\begin{align*}
& \text{ and } \quad \mathcal{P}(s,s') = 1 \iff \boldsymbol{\mathcal{P}}(s,s') = 1.
\end{align*}

The pMC $\M = ({S}, {s_I}, X, \boldsymbol{P})$ is a linear proportional parametrization of $M$ over the parameters set $X$ iff $\M$ is a parametrization of $M$ over $X$ and

\begin{align*}
& \boldsymbol{\mathcal{P}}(s_m, s') {=} x \in X \text{ for some }s' {\in} succ(s_m), \\
& \boldsymbol{\mathcal{P}}(s_m, s'') = (1 - x) \cdot \frac{\mathcal{P}(s_m, s'')}{1 - \mathcal{P}(s_m,s')} \text{ and for } s'' {\neq} s' \in succ(s_m), \\
& \text { and } \boldsymbol{\mathcal{P}}(s, \cdot) = \mathcal{P}(s, \cdot) \text{ for other states}.
\end{align*}

Ratio-preserving. The proportional parametrization preserves the ratio between the original outgoing transition probabilities of each state of $\M$ in the obtained pMC $\M_p$, i.e.,

\begin{equation*}
\text{for }s \in S \text{ and } s, s' \in succ(s), \frac{\mathcal{P}(s,s')}{\mathcal{P}(s,s'')} = \frac{\boldsymbol{\mathcal{P}}(s,s')}{\boldsymbol{\mathcal{P}}(s,s'')}.
\end{equation*}

Well-defined. Each instantiation $u$ with $0 \leq u(x) < 1$ is well-defined for the pMC $\M$ that is obtained by a proportional parametrization
and each instantiation $u$ with $0 < u(x) < 1$ is graph-preserving.
Note that for the state $s \in S$ with only two outgoing successors $s', s''$, the parametric transition probabilities induced by the proportional parametrization scheme are simply $x$ and $1-x$ for some $x \in X$.

Reachability constraints.
A reachability constraint $\varphi = (T, \sim, \lambda)$ is defined over pMC $\mathcal{M} = (S, s_I, X, \mathcal{P})$ with targets $T \subseteq S$, threshold $0 \leq \lambda \leq 1$, and $\sim \ \in \{\leq, \geq\}$.
For pMC $\mathcal{M}$ and instantiation $u: X \rightarrow \mathbb{R}$:

\begin{align*}
\mathcal{M}[u] \models \varphi \quad \text{if and only if} \quad \Pr_{\mathcal{M}[u]}(\lozenge T) \sim \lambda,
\end{align*}

where $\Pr_{\mathcal{M}[u]}(\lozenge T)$ denotes the probability that MC $\mathcal{M}[u]$ reaches a state in $T$.

Consider the pMC $\mathcal{M}_{\B, E}$ in Fig. <ref>. Let .165!$T = \{ s_{10}\}$. Then, .7!$\Pr_{\M_{\B,E}}(\lozenge s_{10}) = \dfrac{361}{34900 \cdot p \cdot q+8758 \cdot q+361}$.

\begin{align*}
\Pr_{\M_{E}}(\lozenge s_{10}) = \dfrac{361}{34900 \cdot p \cdot q+8758 \cdot q+361}.
\end{align*}

By the pMC semantics,

\begin{align*}
\centering
\begin{split}
& p_{s_0} = \covidyes \cdot p_{s_1} + \covidno \cdot p_{s_2} \, \\
&  p_{s_1} = \sympyescovidyes \cdot p_{s_3} + \sympnocovidyes \cdot p_{s_4} \,  \\
& p_{s_2} = \sympyescovidno \cdot p_{s_5} + \sympnocovidno\cdot p_{s_6} \, \\
&  p_{s_3} = \antigenposYYPAR \cdot p_{s_7} + ({\antigennegYYPAR}) \cdot p_{s_0} \, \\
&  p_{s_4} = \antigenposYN \cdot p_{s_7} +\antigennegYN\cdot p_{s_0} \, \\
\end{split}
\begin{split}
&  p_{s_5} = \antigenposNY \cdot p_{s_8} +\antigennegNY \cdot p_{s_0} \, \\
&  p_{s_6} = \antigenposNN \cdot p_{s_8} +\antigennegNN \cdot p_{s_0} \, \\
&  p_{s_7} = \PCRposCOVIDYPAR \cdot p_{s_9} +(\PCRnegCOVIDYPAR) \cdot p_{s_0} \, \\
&  p_{s_8} = \PCRposCOVIDN \cdot p_{s_{10}} +\PCRnegCOVIDN \cdot p_{s_0} \, \\
& p_{s_{9}} = 0, \, \, p_{s_{10}} = 1    \, \\
\end{split}
\end{align*}

\begin{align*}
\begin{split}
& p_{s_0} = \covidyes \cdot p_{s_1} + \covidno \cdot p_{s_2} \, \\
&  p_{s_1} = \sympyescovidyes \cdot p_{s_3} + \sympnocovidyes \cdot p_{s_4} \,  \\
& p_{s_2} = \sympyescovidno \cdot p_{s_5} + \sympnocovidno\cdot p_{s_6} \, \\
&  p_{s_3} = \antigenposYY \cdot p_{s_7} + \antigennegYY \cdot p_{s_0} \, \\
&  p_{s_4} = \antigenposYN \cdot p_{s_7} +\antigennegYN\cdot p_{s_0} \, \\
\end{split}
\begin{split}
&  p_{s_5} = \antigenposNY \cdot p_{s_8} +\antigennegNY \cdot p_{s_0} \, \\
&  p_{s_6} = \antigenposNN \cdot p_{s_8} +\antigennegNN \cdot p_{s_0} \, \\
&  p_{s_7} = \PCRposCOVIDY \cdot p_{s_9} +\PCRnegCOVIDY \cdot p_{s_0} \, \\
&  p_{s_8} = \PCRposCOVIDN \cdot p_{s_{10}} +\PCRnegCOVIDN\cdot p_{s_0} \, \\
& p_{s_{9}} = 0, \, \, p_{s_{10}} = 1    \, \\
\end{split}
\end{align*}

Note that for the pMC $\M'_E$ of the pBN $\B'$, (see Ex. <ref>, Fig. <ref>)

\begin{align*}
\Pr_{\M'_{E}}(\lozenge s_{10}) = \dfrac{3800 \cdot (t \cdot (9 \cdot t + s))}{34200 \cdot t^2+349 \cdot p \cdot q+151 \cdot r \cdot q+3800 \cdot s \cdot t}.
\end{align*}

Let $R \subseteq \mathbb{R}^{n}$ with $n = |X|$ and let

\begin{align*}
\mathcal{M}, R \models \varphi \quad \text{if and only if} \quad \forall u \in R. \ \mathcal{M}[u] \models \varphi.
\end{align*}

We now define the parameter synthesis problems for pMC $\M$ and reachability constraint $\varphi$ that are relevant to our setting.

Partition region $R$ into $R_{+}$, $R_{-}$, and $R_{?}$ such that:

\begin{align*}
& R_{+} \, \subseteq \,
\underbrace{\{u \in R \mid \mathcal{M}[u] \models \varphi\},}_{\mbox{\footnotesize satisfying instantiations}} \quad  R_{-} \, \subseteq \,
\underbrace{\{u \in R \mid \mathcal{M}[u] \models \neg \varphi\},}_{\mbox{\footnotesize refuting instantiations}}
\end{align*}

and $R_{?} = R \setminus (R_{+} \cup R_{-})$ with $||R_{?}|| \leq (1{-}\eta){\cdot}||R||$ for some given coverage factor $0 \leq \eta \leq 1$.

The sub-region $R_{?}$ denotes the fragment of $R$ that is inconclusive for $\varphi$.
This fragment should cover at most fraction $1{-}\eta$ of $R$'s volume.
Region partitioning is exact if $R_{?} = \emptyset$.

For the region $R$ and the specification $\phi$, the problem is to check whether:

\begin{align*}
& \underbrace{\mathcal{M}, R \models \varphi}_{R \text{ is accepting}} \ \text{ or } \   \underbrace{\mathcal{M}, R \models \neg \varphi}_{R \text{ is rejecting}} \  \text{ or } \
\underbrace{\mathcal{M}, R \not \models \varphi \,\wedge\, \mathcal{M}, R\not \models \neg \varphi}_{R \text{ is inconclusive}}.
\end{align*}

For the Markov chain $M$, the parametric Markov chain $\M$ (of $M$), the original instantiation $u_0$, and the reachability specification $\varphi$, the minimal change parameter tuning problem is defined analogously to Def. <ref>.

4.1. Parameter Lifting. The parameter lifting algorithm (PLA) [23] is an abstraction technique that reduces region verification to a probabilistic model checking problem for which efficient algorithms exist [16].
It first removes parameter dependencies by making all parameters unique.
Then it considers for each parameter $x_i$ only its bounds $lb_{x_i}$ and $ub_{x_i}$ within the region $R$.
This yields a non-parametric Markov Decision Process (MDP).

A Markov decision process (MDP) is a tuple $\mathcal{M} = (S,s_I,\mbox{\it Act},\mathcal{P})$ with a finite set $S$ of states, an initial state $s_I \in S$, a finite set $\mbox{\it Act}$ of actions, and a (partial) transition probability function $\mathcal{P} : S \times \mbox{\it Act} \rightarrow \mbox{\it Distr}(S)$.

While resolving the non-determinism in the obtained MDP, a trade-off between maximizing or minimizing $x_i$ may occur. This occurs in particular when parameter $x_i$ repeats in the outgoing transitions of multiple states, see e.g., parameter $t$ in Fig. <ref>.
PLA handles the issue by a relaxation step that introduces intermediate parameters; see e.g.,  Fig. <ref>, where parameter $t$ in the outgoing edges of state $C{=}no, A{=}pos$ is replaced by $t'$.
The relaxation step yields an over-approximation of the region verification problem.

Consider the pMC in Fig. <ref>(a) and  region $t \in [\antigenposNNMin, \antigenposNNMax]$.
Fig. <ref> shows the pMC after relaxation, e.g., parameter $t$ in the outgoing transitions of state $C{=}no,A{=}pos$ is replaced by $t'$.
Fig. <ref> shows the MDP obtained by substituting the parameters with their extremal values.
Consider e.g., state $C=no,S=no$.
Its outgoing dotted (blue) edges have probability $\antigenposNNMin$ and $1{-}\antigenposNNMin$ obtained from $lb(t)$. The probabilities $\antigenposNNMax$ and $1{-}\antigenposNNMax$ of the dashed (purple) edges stem from $ub(t)$.

After substitution, region verification reduces to a simple model-checking query that over-approximates the original higher-order polynomial, e.g., Example <ref>.
The region accepts $\varphi$ for the pMC if the resulting MDP satisfies $\varphi$ for all schedulers.
If all the schedulers satisfy $\neg \varphi$, the region is rejecting. If the region is neither accepting nor rejecting (i.e., inconclusive), it is partitioned into smaller subregions that are more likely to be either accepting or rejecting.

Partitioning ends when at least $\eta\%$ of the region is conclusive.

4.2. Parametric Markov Chain for pBNs.
To enable the use of PLA to parameter tuning of pBNs, we map pBNs onto pMCs as proposed in [26, 27]; see the next example.

Variables states. For the variable $v_i \in V$, let $D^*_{v_i} = D_{v_i} \cup \{*\}$. A variables state is a function $s \colon D^*_{v_1} \times \ldots \times D^*_{v_m}$ that maps each variable $v_i \in V$ either to some $d_i$ or to the don't care value. Intuitively speaking, $v_i = *$ if the value of variable $v \in V$ is either not yet determined or not needed any longer. Let $\mbox{\it States}(V)$ denote the set of all variables states.
Let $s \in \mbox{\it States}(V)$.

* We write $s \models (v_i = d_i)$ iff $s(v_i) = d_i$ with $d_i \in D_{v_i}$.

* We have $s \models (v_i = s(v_i) \wedge \ldots \wedge v_m = s(v_m))$.

* For $V' \subseteq V$ and $s' \in Eval(V')$, we write

\begin{align*}
s \models s' \text{ iff } \bigwedge_{v \in V'} s(v_i) = s'(v_i).
\end{align*}

* For the logical formula $\alpha$ over the atomic propositions $(v_i = d_i)$, we write

\begin{align*}
s \models \alpha \text{ iff }  (v_i = s(v_i) \wedge \ldots \wedge v_m = s(v_m)) \implies \alpha.
\end{align*}

Given the pBN $\mathcal{B} = (V,W,X,\Theta)$, the topological order $\varrho$ on $V$, and the state $s \in \mbox{\it States}(V)$, let $s[v_i{=}d_i] = s^*$, where

* $s^* \in \mbox{\it States}(V)$,

* $s^*(v_i) = d_i$,

* $s^*(v) = s(v)$ for $v \neq v_i$ satisfying

\begin{align*}
\varrho(v) < \varrho(v_i) \text{ and } \exists c \in \mbox{\it children}(v) \,.\, \varrho(c) > \varrho(v_i), \text{ and}
\end{align*}

* $s^{*}(v) = *$ for all other variables.

Let $\mathcal{B}=(V,W,X, \Theta)$ be a pBN and $\varrho = (v_1, \dots, v_m)$ be a topological order over $V$. The pMC of $\mathcal{B}$ is $\mathcal{M}^{\varrho}_{\mathcal{B}}=(S, s_0, X, P)$, where:

* $s_0 =(v_1=*, \ldots, v_m=*)$ is the initial state,

* $S \subseteq \mbox{\it States}(V)$ is the set of states defined as follows:

* $S = \bigcup\limits_{0 \leq i \leq m} S_i$

* $S_0 = \{s_0\}$

* $S_{i{+}1} = \{ s_{i}[v_{i{+}1}{=}d_{i{+}1}] \,\mid\, \text{for } s_{i} \in S_{i} \text{ and } d_{i{+}1} \in D_{v_{i{+}1}}\} \text{ for } 0 \leq i < m$

* and $P \colon S \times S \to [0,1]$ is the transition probability functions defined by the following rules.

* Let $0 \leq i < m$, $\overline{par} \in \mbox{\it Eval}(\mbox{\it parents}(v_{i{+}1}))$, and $s_{i{+}1} = s_{i}[v_{i{+}1}{=}d_{i{+}1}] \in S_{i{+}1}$. For each $s_{i} \in S_{i}$ and $d_{i{+}1} \in D_{v_{i{+}1}}$,

\begin{align*}
P(s_{i}, s_{i{+}1}) = \Theta(\overline{par})(d_{i{+}1}) \text{ iff } s_{i} \models \overline{par}.
\end{align*}

* $P(s_m, s_m) = 1$ for each $s_m \in S_m$.

* $P(s_i, s_j) = 0$, otherwise.

Consider the acyclic pMC in Fig. <ref> for the pBN in Fig. <ref> for the topological ordering $\varrho: C < S < A < P$.

Initially, all variables are don't care.
The initial state can evolve into $\mbox{\it C}{=}yes$ and $\mbox{\it C}{=}no$.
Its transition probabilities $\covidyes$ and $\covidno$ come from the CPT entries of .
Generally, at “level” $i$ of the pMC, the outgoing transitions are determined by the CPT of the $i{+}1$-st variable in $\varrho$.

The following proposition [27] relates inference in pBNs to reachability objectives in pMCs.

For the pBN $\B$, the topological order $\varrho$, and the pMC $\M_{\B}^{\varrho}$,

\begin{align*}
&\Pr_{\B}(E) \ = \ 1 - \Pr_{\M^{\varrho}_{\B}}(\lozenge \, \neg E)
\quad \text{and} \\
\Pr_{\B} (H \mid E) \ = \ \dfrac{1 - \Pr_{\M^\varrho_\B}(\lozenge\,(\neg H \lor \neg E))}{1 - \Pr_{\M^{\varrho}_{\B}}(\lozenge \, \neg E)},
\end{align*}

where the latter equality requires $\Pr_{\B}(\neg E) < 1$.

pBN constraints relate to reachability constraints in its pMC.

Consider the pBN $\B_1$ from Fig. <ref> and the query .55!$\Pr_{\B}(\text{C} = \text{no} \,|\,\text{A} = \text{pos} \wedge \text{P} = \text{pos} )$ from Example <ref>. Let $\M^{\varrho}_{\B}$ be the pMC in Fig. <ref>. The query coincides with

\begin{align*}
& \dfrac{1 - \Pr_{\M^\varrho_\B}(\lozenge\, \text{C}=\text{yes} \lor \text{A} = \text{neg} \lor \text{P} = \text{neg} )}{1 - \Pr_{\M^\varrho_\B}(\lozenge\, \text{A} = \text{neg} \lor \text{P} = \text{neg})} \\
%& = \dfrac{1- (0.999658 -0.0038 \cdot q)}{1- (0.9913379 - 0.033155 \cdot p + 0.0038 \cdot q) } \\
& = \dfrac{361}{34900 \cdot p \cdot q+8758 \cdot q+361} = f_{\B, \phi}; \text{ see Ex. \ref{ex:constraint}.}
%& = \dfrac{1- 0.999693}{1 - 0.9674473} ~\approx~ 0.011089.
\end{align*}

The iterative procedure of region-based parameter tuning for $n=2$ and $\gamma = \nicefrac{1}{2}$.

§ REGION-BASED MINIMAL CHANGE TUNING

Our approach to the parameter tuning problem for a pBN and constraint $\phi$ consists of two phases.
The $\epsilon$-close partitioning

exploits parameter lifting on the pMC for reachability constraint $\varphi$ of $\phi$.
We start with a rectangular region enclosing the original instantiation $u_0$.
This region is iteratively expanded if it rejects $\varphi$.
Inconclusive regions are partitioned.
Any iteration that finds some accepting subregion ends the first phase.
The second phase (Alg. <ref>)
extracts a minimal-change instantiation from the set of accepting sub-regions.

5.1. Obtaining Epsilon-Close Subregions. Alg.   <ref> is the main thread of our approach.
Its inputs are the pBN $\B=(V,W,X,\Theta)$, the initial instantiation $u_0: X \to \mathbb{R}$, and the constraint $\phi$.
It starts with obtaining the pMC  $\mathcal{M}_{\B}$ of the pBN and reachability constraint $\varphi$ of $\phi$ (lines 1–2). The output is either (i) $u_0$ if $\mathcal{M}[u_0] \models \varphi$, or (ii) the instantiation $u$ with a minimal distance from $u_0$ or (iii) no instantiation if $\varphi$ is infeasible[That is, in the $\eta = 100\%$ of the checked parameter space, no satisfying instantiation has been found.].  The hyper-parameters of the algorithm are the coverage factor $0 < \eta < 1$, the region expansion factor $0 < \gamma < 1$, and the maximum number of iterations $K \in \mathbb{N}$.
The hyper-parameters $\gamma$ and $K$ steer the iterative procedure. They initially determine the distance bound $\epsilon$, see line 3.
The bound $\epsilon$ determines region $R_{\epsilon}$, i.e., how far the region bounds deviate from $u_0$.

PLA then verifies $R_{\epsilon}$. If $R_{\epsilon}$ is rejecting, $\epsilon$ is extended by the factor $\gamma^{-1}$ (l. 9).
Figure <ref> visualizes the procedure for $\gamma{=}\nicefrac{1}{2}$ and $n{=}2$.
At iteration $i$, (a) PLA is invoked on $R_{\epsilon}$ and either (b) $R_+ \neq \emptyset$ or (c) $R_+ = \emptyset$. For case (b), the iterative procedure ends and Alg. <ref> is called to obtain an instantiation in $R_{+}$ that is closest to $u_0$.
For case (c), the region is expanded by factor $\gamma^{-1}$ and passed to PLA, see (d). Note that the loop terminates when the distance bound $\epsilon$ reaches its upper bound $\hat{d_0}$ (l. 2). We refer to Sec. <ref> Corollary 1 and 3 for computing $\hat{d_0}$.

Region expansion schemes. Region $R_{\epsilon}$ is determined by factor $\epsilon$, see line 16. The methods and detail this for each of our distance measures.
For the EC distance (line 20), the parameters have an absolute deviation from $u_0$. We refer to Corollary 2, Sec. <ref>. For CD distance (lines 23 and 24), the deviation is relative to the initial value of the parameter. We refer to Corollary 4, Sec. <ref>.  Such deviation schemes ensure $\epsilon$-closeness of $R_{\epsilon}$ both for EC and CD distance, that is, all the instantiations in $R_{\epsilon}$ have at most distance $\epsilon$ from $u_0$.

For pBNs with a single parameter, already checked regions can be excluded in the next iterations.
Applying such a scheme to multiple parameters yet may not yield minimal-change instantiation.
To see this, let $u_0(x) = 0.15$ and $u_0(y) = 0.15$ and take $\epsilon_0 = 0.05$.
Assume $R_{\epsilon_0} = [0.1,0.2] \times [0.1,0.2]$ is rejecting.
Limiting the search in the next iteration to $R_{\epsilon_1} = [0.2,0.4] \times [0.2, 0.4]$ may omit the
accepting sub-regions that are at a closer distance from $u_0$, e.g. if the region $[0.1,0.2] \times [0.2, 0.4]$ includes some satisfying sub-regions. This is why we include the already-analyzed intervals starting from $u_0$ in the next iterations. However, as will be noted in our experimental results, this does not have a deteriorating effect on the performance: the unsuccessful intermediate iterations are normally very fast as they are often analyzed by a single model checking: the model checker in a single step determines that the entire region is rejecting and no partitioning is needed.

5.2. Obtaining a Minimal-Distance Instantiation. Once $R_+ \neq \emptyset\subseteq R_{\epsilon}$ is found, it remains to find the instantiation $u_+$ in $R_{+}$ that is closest to $u_0$.
Region $R_+$ includes infinitely many instantiations, yet minimal-distance instantiation is computable in $O(k{\cdot} n)$ for $n = |X|$ and the regions $R_{+} = \{R_{+,1}, \cdots, R_{+,k}\}$: (i) PLA ensures that the subregions $R_{+,1}$ to $R_{+,k}$ are rectangular and mutually disjoint. (ii) Due to the triangle equality of the distance measure, a minimum-distance instantiation can be picked from the bounded region.
Alg. <ref> simplifies finding $u_+$.

In the first step, we pick $u_{+,i}$ from $R_{+,i}$ that has minimal distance from $u_0$ (l. 3–18).
The idea is simple: for every $x_j \in X$, three cases are possible, see Fig. <ref>.
Let $lb_{x_j}$ be the lower bound for parameter $x_j$ in the region $R_{+,i}$ and $ub_{x_j}$ be its upper bound.
(1) If $u_0(x_j) \in [lb_{x_j},ub_{x_j}]$, we set $u_{+,i}(x_j) = u_0(x_j)$. This yields the least distance, i.e., $0$ in the dimension of $x_j$, see Fig. <ref> (left).
(2) $u_0(x_j) < lb_{x_j}$, see Fig. <ref> (middle).
Then $lb_{x_j}$ has the least distance from $u_0(x_j)$ in the dimension $x_j$, i.e., $|lb_{x_j}{-}u_0(x_j)| < |a{-}u_0(x_j)|$ for every value $a \in [lb_x, ub_x]$  and similarly for CD-distance, $\ln{\frac{lb_{x_j}}{u_0(x_j)}} < \ln{\frac{a}{u_0(x_j)}}$.
In this case, we set $u_{+,i}(x_j) = lb_{x_j}$ to minimize the distance in dimension $x$.
(3) By symmetry, for $u_0(x_j) > ub_{x_j}$, we set $u_{+,i}(x_j) = ub_{x_j}$; see Fig. <ref> (right).
It remains to compute the distance for each candidate in $U_+ = \{u_{+,1}, \cdots, u_{+,k}\}$ and pick $u_+$ closest to $u_0$ (l. 19).

Obtaining the minimal-change value for parameter $x \in X$.

We apply the entire algorithm on the pBN in Fig. <ref>, constraint $\phi = \Pr(\text{C}{=}\text{no} \mid \text{A}{=} \text{pos} \wedge \text{P}{=}\text{pos} ) \leq 0.009$, $u_0: p \rightarrow 0.72, q \rightarrow 0.95$, and $\eta{=}0.99$.
The pMC $\M_{\B,E}$ is given in Fig. <ref> and the reachability constraint is $\varphi: \Pr(\lozenge s_{10}) \leq 0.009$.
We take $\hat{d_0} {=} \sqrt{n}$, $\gamma {=} \nicefrac{1}{2}$, and $\epsilon_0 {=} \nicefrac{1}{16}$.

* The initial region $R_{\epsilon}$ is obtained by $\epsilon{=}\nicefrac{1}{8}$:

\begin{align*}
[0.72{-}\nicefrac{1}{8\sqrt{2}}, 0.72{+} \nicefrac{1}{8\sqrt{2}}] \times [0.95{-} \nicefrac{1}{8\sqrt{2}}, 0.95{+}\nicefrac{1}{8\sqrt{2}}].
\end{align*}

* Running PLA on $R_{\epsilon}$ gives no sub-region accepting

* Expanding the region by the factor $\gamma^{-1} = 2$ yields:

\begin{align*}
[0.72{-}\nicefrac{1}{4\sqrt{2}}, 0.72{+} \nicefrac{1}{4\sqrt{2}}] \times [0.95{-} \nicefrac{1}{4\sqrt{2}}, 0.95{+}\nicefrac{1}{4\sqrt{2}}].
\end{align*}

* Running PLA gives $12$ regions accepting $\varphi$, e.g., $R_{+,1}: [0.92075, 0.960375] \times [0.97475,1]$ and
$R_{+,2}: [0.960375, 0.97028175] \times [0.9431875, 0.9495]$.

* Running Alg. <ref> on $R_{+,1}$ through $R_{+,12}$ gives minimal distance candidates $u_{+,1}$ through $u_{+,12}$; e.g., $u_{+,2}(p){=}lb_p{=}0.960375$ as $u_0(p) < lb_p$ and $u_{+,2}(q){=}ub_q{=}0.9495$ as $ub_q < u_0(q)$; see Fig. <ref>.

* Running Alg. <ref> using the distance measure of choice—here: EC-distance—for the candidates $U_{+}$ returns $u_{+}$ closest to $u_0$: $u_+(p) = 0.92075$ and $u_+(q) = 0.97475$; distance $0.040913125$.

Let us shortly discuss the hyper-parameters: $\eta$, $\gamma$, and $K$. The hyper-parameter $\eta$ is the coverage factor for the region partitioning. Our region partitioning algorithm, i.e., parameter lifting (line 17 Alg. <ref>) works based on this factor: the procedure stops when $R_?$ is at most $(1-\eta) \cdot ||R||$ of the entire region $R$, e.g., for $\eta = 0.99$, it continues until at least $99\%$ of $R$ is either rejecting or accepting and at most $1\%$ is unknown. This relates $\eta$ to the approximation bounds of our algorithm, see Problem statement <ref>. Intuitively speaking, the coverage factor $\eta$ means that if the algorithm provides a minimal-distance value $D(u,u_0) = d$, then with the probability $1-\eta$ there may exist a smaller value distance than $d$ that works too but has not been found. One can thus consider $\eta$ as the confidence factor for the minimality of the results. The hyper-parameter $\gamma$ specifies the factor by which the region is expanded at each iteration, see Line 9, Alg. <ref>. The hyper-parameters $\gamma$ and $K$ specify the size of the initial region, see Lines 3 and 16 of Alg. <ref>. Our experiments (not detailed in the paper) with $\gamma=\{0.2, 0.5, 0.8\}$ and $K=\{2, \cdots, 12\}$ reveal that $\gamma=\nicefrac{1}{2}$ and $K = 6$ gave the best balance. Large values for $\gamma \, (0.8)$ and small $K$ lead to unnecessary computations in the initial iteration for the simple cases i.e., when small perturbations of the parameters make the constraint satisfied. Small values for $\gamma\, (0.2)$ lead to large regions in the next iterations due to the expansion by $\gamma^{-1}$.

§.§ Tailoring the Algorithm for CD-Distance

While the Euclidean distance takes the absolute change of each parameter into account, other distance measures consider the relative change in order to bind the effect that is imposed on the joint probability distribution when altering the parameters.
Our region-based parameter tuning approach is adaptable for other distance measures: (I) the last line of Algorithm <ref> changes to computing distance measures of choice. Second and more importantly, (II) the region expansion acquires a non-uniform scheme. continue and introduce epsilon/gamma as a vector (dimension $n = |X|$ depending on the original parameter values

§.§ Accuracy Bounds

§ EXPERIMENTAL EVALUATION

We empirically evaluated our approach using a prototypical realization on top of the probabilistic model checker Storm [15] (version 1.7.0).
As baseline, we used the latest version (10.4) of Bayesserver[<https://www.bayesserver.com/>], a commercial BN analysis tool that features sensitivity analysis and parameter tuning for pBNs.
It supports pBNs with a single parameter only.
We parametrized benchmarks from  repository and defined different constraints.
We (i) parametrized the CPTs of the parents (and grandparents) of the evidence nodes, and (ii) used the SamIam tool [<http://reasoning.cs.ucla.edu/samiam/>]

to pick the CPT entries $\Theta_{\modif}$ most relevant to the constraint.
To get well-formed BNs, we used the linear proportional co-variation, see Def. <ref>.
In all experiments, we picked evidences from the last nodes in the topological order; they have a long dependency path to the BN roots.
This selection reflects the worst-case in [25].

Previous studies [25] indicate that inference on BNs using probabilistic model checking is, in practice, computationally more expensive when the query nodes are in the last nodes of the topological order. To fairly reflect the behavior of our approach in the evaluation results, we accordingly picked the evidence nodes from the last nodes in the topological order with a long dependency path to the BN roots.

We took $\gamma{=}\nicefrac{1}{2}$ and $K{=}6$ for our experiments, see Alg.<ref>.  We conducted all our experiments on a 2.3 GHz Intel Core i5 processor with 16 GB RAM.

RQ1: Comparison to Bayesserver.
We used small, medium, and large BNs from the  repository.
For each BN, we parameterized one distribution in a single CPT to align with the restrictions from Bayesserver.
Figure <ref> indicates the results, comparing the tuning times (in sec) of Bayesserver (x-axis) to those of our implementation (y-axis).
The latter includes the time for region refinements by PLA.
For each pair (pBN, constraint), we did experiments for coverage factors $\eta = 1{-}10^{-i}  \text{ for } i \in \{1, \cdots, 6\}$.
The influence of $\eta$ on the exactness of our results is addressed under RQ2.

Storm outperforms Bayesserver for several benchmarks (, , , and multiple instances of ), whereas Bayesserver is faster by about an order of magnitude for the other pBNs, such as . Explanation: Bayesserver exploits specific methods for one-way sensitivity analysis and relies on the linear form of the sensitivity function. These techniques are very efficient yet not applicable to pBNs with multiple parameters in multiple CPTs. For our experiments on those subclasses, see RQ4. Such subclasses are not supported by Bayesserver and—to the best of our knowledge—not any existing BN tool. This applies e.g., also to Bayesfusion which considers only the change of a single parameter at a time, and SamIam which is limited to the single parameter and single-CPT.

RQ2: Sensitivity to the coverage factor $\eta$.
For each pBN and constraint, we decreased (the refinement factor) $1{-}\eta$ in a step-wise manner by a factor $10^{-1}$.
To quantify the tightness of our results, we measured how our approximately-close instantiation, denoted $u_{\storm}$, differs from the absolute minimum-distance from Bayesserver, denoted $u_{\bayesserver}$.
Figure <ref> (log-log scale) plots the tightness of the results $|u_{\storm} - u_{\bayesserver}|$ (y-axis) against the refinement factor (x-axis).
Figures <ref> and <ref> indicate the number of iterations and the tuning time (log scale, seconds) for each refinement factor.

(I) Mostly, $1{-}\eta$ bounds the difference between our approximately-close solution and Bayesserver's solution.

For e.g., $\eta = 1 {-} 10^{-4}$, the difference is at most $10^{-4}$.
(II) On increasing the coverage, the difference to the true minimal distance rapidly decreases.
((III) The computation time moderately increases on increasing the coverage, but the number of iterations was mostly unaffected. Explanation: (I, II) Recall that the value of $1{-}\eta$ bounds the size of the unknown regions; see Def. <ref>. This indicates why $|u_{\storm} - u_{\bayesserver}|$ relates to $1{-}\eta$. (III) At a higher coverage factor, the region partitioning is more fine-granular possibly yielding more accepting regions to analyze. Therefore the computation becomes more expensive. The timing is, however, not correlated to the number of iterations. This is because the $\epsilon$-close iterations before the last iteration are often completed by a single region verification and are very fast. Similar observations have been made for $n{>}1$ parameters; see RQ4.

RQ3: Sensitivity to the threshold $\lambda$.
We varied the constraint's threshold ($\lambda$) by steps of $\nicefrac{1}{20}$ for the benchmarks  , and  with $n{=}1$ and $\eta{=}99.99999\%$.
Figures <ref>, <ref>, and <ref> display the outcomes with the x-axis indicating the threshold and the y-axis indicating the tuning time (in seconds), the distance (log-scale), and the number of iterations.

Findings: By strengthening the threshold, the possibly satisfying regions get further away from $u_0$. Thus the distance, the number of iterations, and sometimes the tuning time grow. Similar findings are valid for $n>1$ parameters; see RQ4, Table <ref>. Explanation: Region refinement starts with small regions in the close vicinity of the original values of the parameters. Therefore, for the constraints close to the original probability $\Pr_{\B[u_0], \phi}$, the number of iterations is low, the distance is naturally small, and the minimal-change tuning is completed faster without the need to analyze larger regions.

RQ4: Scaling the number of parameters.
We took the  and  benchmarks and parameterized them in multiple ways.
Their pMCs have $9,947$ and $1,118$ states and $17,948$ and $3,147$ transitions respectively.
The set of parameters for each pBN is including (and doubles the number of) parameters in the previous pBN.

Table <ref>(left) and (right) list the results for  and .
We list for each pBN, the number of affected CPTs, the number of parameters, the threshold $\lambda$, and the coverage $\eta$.
E.g.,  with $par{=}8$ has $8$ parameters occurring in 4 CPTs.
The columns EC, iter, and $t(s)$ report the EC-distance, the number of iterations, and the total time in seconds (incl. model building time, time for region refinement, and tuning time) respectively.
TO and MO indicate time-out (30 minutes) and memory-out ($>$16 GB).

(I) Approximately-close parameter tuning is feasible for pBNs with up to $8$ parameters.
This is significantly higher than the state of the art—one parameter.
As the number of sub-regions by PLA grows exponentially, treating more parameters is practically infeasible.
(II) More parameters often ease finding satisfying instantiations. E.g., the threshold $\geq 0.25$ is unsatisfiable for  with $n{=}2$, but is satisfied with $n{=}4$.
(III) The results for multiple parameter pBNs confirm the findings for RQ2 and RQ3; see the rows for each pBN instance with (i) varied coverage and (ii) varied threshold.
(IV) The unsatisfiability of a constraint can be computed fast (with $100\%$ confidence), regardless of the number of parameters, see e.g.,  with $8$ parameters for the constraint $\leq 0.001$.
For infeasibility, a single verification suffices; no partitioning is needed.

RQ5: Handling pBNs with parameter dependencies. Parameter lifting algorithm (Section 4.1) enables handling models with parameter dependencies, see e.g., Example <ref>: the parameters are allowed in multiple local distributions and the pBN sensitivity function is of higher degree. We extended our experiments to such cases for  and : we parameterized the entries $\theta_1 \cdots \theta_k \in \Theta_{\modif}$ over the same parameter $x \in X$ when the original values of $\theta_1 \cdots \theta_k$ were the same in the original BN. The eleventh row in Table <ref> (left) and the eighth and eleventh rows in Table <ref> (right) correspond to such cases. The term $8^*(2)(2)$ e.g., denotes that out of the $8$ parameters, 2 parameters repeatedly occurred in two distinct distributions. The term $8^*(3)$ denotes that out of the $8$ parameters, one was occurring in 3 distinct distributions.

Findings: (I) Our method is applicable to pBNs with parameter dependencies where the sensitivity function is of a higher degree. (II) For the same coverage factor, the pBNs with parameter dependency are more expensive to analyze. See e.g., the two rows for $\texttt{win95pts}$ with $8$ parameters, threshold $\geq 0.115$, and the coverage $80\%$. This is due to more complex sensitivity functions that give a higher number of sub-regions to verify. (III) The pBNs with parameter dependency yielded notably smaller distances.

§ EPILOGUE

Related work.
Kwisthout and van der Gaag [[18]]
studied the theoretical complexity of tuning problems. Renooij [[24]] studied the properties of the co-variation schemes for BN tuning. She shows that the linear proportional scheme optimizes the CD distance for single-CPT pBNs. Similar are the studies by Bolt and van der Gaag [[3, 4]] that consider tuning heuristics for distance optimizations. Peng and Ding [[22]] propose an iterative proportional fitting procedure (IPFP) to minimize the KL-divergence distance for a set of constraints. The method does not scale to large networks. Santos et al. [[28]] exploits linear programming for BN parameter tuning. Yak- aboski and Santos [[30]] consider a new distance measure for parameter learning and tuning of BNs. Leonelli [[19]] considers nonlinear sensitivity functions, yet only for single parameter pBNs and Ballester-Ripoll and Leonelli [[1]] efficiently compute the derivatives of sensitivity functions to select the most relevant parameters to a query, yet they limit to single parameter variation.

Related work, unstructured:
Model repair tunes the transition probabilities in an MC to satisfy a constraint $\varphi$ while optimizing a cost function.
<cit.> reduce it to a non-linear optimization problem, whereas <cit.> exploit abstraction.

For Bayesian networks:

[24] studies the co-variation schemes for tuning, in particular the characteristics of the proportional scheme and proves that it optimizes the CD distance for single CPT.
Similar are: [3] <cit.> [4] that consider new heuristics for parametrization.

Relevant to us: [19] sensitivity analysis beyond linearity. Considers non-linear sensitivity function for a single parameter

<cit.>: an iterative proportional fitting procedure (IPFP), minimizes the KL-divergence distance for tuning with multiple constraints. does not scale.

[18] theoratical complexity of different variants of BN tuning problems

Most Recent work (that I've found): [1] that computes the sensitivity values to quantify the affect of a single parameter) and to find the most relevant parameters for a given query. Does not go beyond sensitivity analysis and is not applied to the tuning of multiple parameters

[6] proposes approaches for sensitivity analysis and tuning of pBNs with parameters in a single CPT and extends the theoretical findings to more general pBNs. They show that pBNs with more parameters often lead to smaller distances. (confirmed in our experiments)

To our knowledge, no prior work practically scales to large networks for minimal tuning of multiple parameters in multiple CPTs, and no prior work on applying region verification techniques for the minimal change tuning.

* Darwiche papers <cit.>, distance measures [7], single CPT [6]
<cit.> These are mostly cited through the text

* Sensitivity functions, related papers of Van der Gaag and L.C., Renooij These can be cited simply as the ground basis papers for parameter tuning and sensitivity analysis [9] <cit.> [18] [24]

* Our own previous work (be carefull about the double-blind review)

* tuning heuristics, Bolt

These are tuning heuristics paper that we can write a sentence in the conclusion, that remain for future work[3] <cit.> [4]

[28] exploits linear programming for BN parameter tuning [30] consideres a new distance measure for parameter learning and tuning of BNs

<cit.> a survey on how to deal with uncertainties in CPT entries. Could be cited in the intro.

For the rest, I should look more to select which is more important to include

* Credal networks <cit.> <cit.> <cit.>, Interval Bayesian networks (?),

* Robustness analysis <cit.> <cit.>, <cit.>, <cit.>, <cit.>

This seems relevant to me. I wrote a bit about it, above
[1] 2022 paper, but still limited to single parameter fluctuations in the experiments and limited to sensitivity values.

Conclusion.A novel algorithm for parameter tuning in Bayesian networks is presented and experimentally evaluated.
Whereas existing algorithms come with severe restrictions—single parameters and/or linear functions—our approach is applicable to multiple (in practice about 8) parameters, large BNs (up to 100 variables), and polynomial functions.

Future work includes considering balanced tuning heuristic [4] and using monotonicity of parameters [29].

Future work and ideas:

* Exploiting the monotonicity of the parameters to optimize our approach to make the region verification and region splitting more efficient in practice by having smaller regions

* Investigating new heuristics for taking the region candidates and new region extension schemes, e.g., non-uniform or threshold-dependent extension schemes.

* To investigate the optimal candidates for varied parameters and their inter-dependency

* To obtain accuracy bounds for the region-based tuning method formally

* Addressing other distance measures such as KL-divergence and specific distance measures for n-way parameter tuning with parameter dependencies. I think KL would be a problem, as it is not a "distance measure" and does not satisfy triangle inequality which is needed for

§ ACKNOWLEDGEMENT

This research was funded by the ERC AdG Projekt
FRAPPANT (Grant Nr. 787914). We kindly thank Alexandra Ivanova for her implementation efforts and Tim Quatmann for the fruitful discussions.

[1]
Rafael Ballester-Ripoll and Manuele Leonelli.
You only derive once (YODO): automatic differentiation for
efficient sensitivity analysis in Bayesian networks.
In PGM, volume 186 of Proceedings of Machine Learning
Research, pages 169–180. PMLR, 2022.

[2]
Pablo Barreiro, Jesús San-Román, Maria del Mar Carretero, and
Francisco Javier Candel.
Infection and infectivity: Utility of rapid antigen tests for the
diagnosis of covid-19.
Revista Española de Quimioterapia, 34(Suppl1):46, 2021.

[3]
Janneke H. Bolt and Linda C. van der Gaag.
Balanced tuning of multi-dimensional Bayesian network classifiers.
In ECSQARU, volume 9161 of Lecture Notes in Computer
Science, pages 210–220. Springer, 2015.

[4]
Janneke H. Bolt and Linda C. van der Gaag.
Balanced sensitivity functions for tuning multi-dimensional
Bayesian network classifiers.
Int. J. Approx. Reason., 80:361–376, 2017.

[5]
Enrique F. Castillo, José Manuel Gutiérrez, and Ali S. Hadi.
Sensitivity analysis in discrete Bayesian networks.
IEEE Trans. Syst. Man Cybern. Part A, 27(4):412–423, 1997.

[6]
Hei Chan and Adnan Darwiche.
Sensitivity analysis in Bayesian networks: From single to multiple
In UAI, pages 67–75. AUAI Press, 2004.

[7]
Hei Chan and Adnan Darwiche.
A distance measure for bounding probabilistic belief change.
Int. J. Approx. Reason., 38(2):149–174, 2005.

[8]
Hei Chan.
Sensitivity analysis of probabilistic graphical models.
University of California, Los Angeles, 2005.

[9]
Veerle M. H. Coupé and Linda C. van der Gaag.
Properties of sensitivity analysis of Bayesian belief networks.
Ann. Math. Artif. Intell., 36(4):323–356, 2002.

[10]
Murat Cubuktepe, Nils Jansen, Sebastian Junges, Joost-Pieter Katoen, and Ufuk
Synthesis in pMDPs: A tale of 1001 parameters.
In ATVA, volume 11138 of Lecture Notes in Computer
Science, pages 160–176. Springer, 2018.

[11]
Murat Cubuktepe, Nils Jansen, Sebastian Junges, Joost-Pieter Katoen, and Ufuk
Convex optimization for parameter synthesis in MDPs.
IEEE Trans. Autom. Control., 67(12):6333–6348, 2022.

[12]
Adnan Darwiche.
Modeling and Reasoning with Bayesian Networks.
Cambridge University Press, 2009.

[13]
Jacqueline Dinnes, Pawana Sharma, Sarah Berhane, Susanna S van Wyk, Nicholas
Nyaaba, Julie Domen, Melissa Taylor, Jane Cunningham, Clare Davenport, Sabine
Dittrich, et al.
Rapid, point-of-care antigen tests for diagnosis of sars-cov-2
Cochrane Database of Systematic Reviews, 7, 2022.

[14]
Linus Heck, Jip Spel, Sebastian Junges, Joshua Moerman, and Joost-Pieter
Gradient-descent for randomized controllers under partial
In VMCAI, volume 13182 of Lecture Notes in Computer
Science, pages 127–150. Springer, 2022.

[15]
Christian Hensel, Sebastian Junges, Joost-Pieter Katoen, Tim Quatmann, and
Matthias Volk.
The probabilistic model checker storm.
Int. J. Softw. Tools Technol. Transf., 24(4):589–610, 2022.

[16]
Joost-Pieter Katoen.
The probabilistic model checking landscape.
In LICS, pages 31–45. ACM, 2016.

[17]
Doga Kisa, Guy Van den Broeck, Arthur Choi, and Adnan Darwiche.
Probabilistic sentential decision diagrams.
In KR. AAAI Press, 2014.

[18]
Johan Kwisthout and Linda C. van der Gaag.
The computational complexity of sensitivity analysis and parameter
In UAI, pages 349–356. AUAI Press, 2008.

[19]
Manuele Leonelli.
Sensitivity analysis beyond linearity.
Int. J. Approx. Reason., 113:106–118, 2019.

[20]
Hiroshi Nishiura, Tetsuro Kobayashi, Takeshi Miyama, Ayako Suzuki, Sung-mok
Jung, Katsuma Hayashi, Ryo Kinoshita, Yichi Yang, Baoyin Yuan, Andrei R
Akhmetzhanov, et al.
Estimation of the asymptomatic ratio of novel coronavirus infections
Int. Journal of Infectious Diseases, 94:154–155, 2020.

[21]
Judea Pearl.
Probabilistic reasoning in intelligent systems: networks of
plausible inference.
Morgan kaufmann, 1988.

[22]
Yun Peng and Zhongli Ding.
Modifying Bayesian networks by probability constraints.
In UAI, pages 459–466. AUAI Press, 2005.

[23]
Tim Quatmann, Christian Dehnert, Nils Jansen, Sebastian Junges, and
Joost-Pieter Katoen.
Parameter synthesis for Markov models: Faster than ever.
In ATVA, volume 9938 of Lecture Notes in Computer
Science, pages 50–67, 2016.

[24]
Silja Renooij.
Co-variation for sensitivity analysis in Bayesian networks:
Properties, consequences and alternatives.
Int. J. Approx. Reason., 55(4):1022–1042, 2014.

[25]
Bahare Salmani and Joost-Pieter Katoen.
Bayesian inference by symbolic model checking.
In QEST, volume 12289 of Lecture Notes in Computer
Science, pages 115–133. Springer, 2020.

[26]
Bahare Salmani and Joost-Pieter Katoen.
Fine-tuning the odds in Bayesian networks.
In ECSQARU, volume 12897 of Lecture Notes in Computer
Science, pages 268–283. Springer, 2021.

[27]
Bahare Salmani and Joost-Pieter Katoen.
Fine-tuning the odds in Bayesian networks.
CoRR, abs/2105.14371, 2021.

[28]
Eugene Santos, Qi Gu, and Eunice E. Santos.
Bayesian knowledge base tuning.
Int. J. Approx. Reason., 54(8):1000–1012, 2013.

[29]
Jip Spel, Sebastian Junges, and Joost-Pieter Katoen.
Are parametric Markov chains monotonic?
In ATVA, volume 11781 of Lecture Notes in Computer
Science, pages 479–496. Springer, 2019.

[30]
Chase Yakaboski and Eugene Santos.
Bayesian knowledge base distance-based tuning.
In WI, pages 64–72. IEEE Computer Society, 2018.