Finding an ϵ-Close Minimal Variation of Parameters in Bayesian Networks

1 Introduction

Bayesian networks. Bayesian networks (BNs) [Pearl, 1988] are probabilistic graphical models that enable succinct knowledge representation and facilitate probabilistic reasoning [Darwiche, 2009]. Parametric Bayesian networks (pBNs) [Castillo et al., 1997] extend BNs by allowing polynomials in conditional probability tables (CPTs) rather than constants.

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 [Heck et al., 2022], convex optimization [Cubuktepe et al., 2018; Cubuktepe et al., 2022], and region verification [Quatmann et al., 2016]. Recently, Salmani and Katoen [2021a] 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 [Kwisthout and van der Gaag, 2008] for Bayesian networks.

Minimal-change parameter tuning. Given a Bayesian network B, a hypothesis H, evidence E, λ ∈ [0, 1], and a constraint of the form Pr(H|E) ≤ λ (or ≥ λ),

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 COVID19 that is adopted from several medical studies [Barreiro et al., 2021; Nishiura et al., 2020; Dinnes et al., 2022]. 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 presymptomatic). Figure 1a 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 ϵ-bounded variant of the problem: for a subset of modifiable CPT rows, are there new values within the distance of ϵ 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 ϵ-bounded tuning. More precisely, we find instantiations that (i) satisfy the constraint (if the constraint is satisfiable) and (ii) are ϵclose and (iii) lean towards the minimum distance instantiation depending on a coverage factor 0 ≤ η ≤ 1.
• We propose two region expansion schemes to realize ϵcloseness of the results both for Euclidean distance and for CD distance [Chan and Darwiche, 2005].
• 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 1 indicate that ϵ-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.

1 https://github.com/baharslmn/pbn-epsilon-tuning

2 Background

Variables. Let V be a set ofmrandomvariables v1, . . . , vm and Dvi the domain of variable vi. For A ⊆ V , Eval(A) denotes the set of joint configurations for the variables in A.

Parameters. Let X={x1, . . . , xn} be a set of n realvalued parameters. A parameter instantiation is a function u :X→R that maps each parameter to a value. All parameters are bounded; i.e., lbxi ≤u(xi)≤ubxi for xi. Let Ii = [lbxi , ubxi ]. The parameter space U ⊆ Rn of X is the set of all possible values of X, i.e., the hyper-rectangle spanned by the intervals Ii for all i.

Substitution. Polynomials f over X are functions f : Rn → R where f(u) is obtained by replacing each occurrence of xi in f by u(xi); e.g., for f = 2x21+x2, u(x1) = 3 and u(x2) = 2, f (u) = 20. Let f [u] denote such substitution.

Parametric distributions. Let Distr(D) denote the set of probability distributions over D. Let Q[X] be the set of multivariate polynomials with rational coefficients over X. A parametric probability distribution is the function µ : D → Q[X] with ∑d∈D µ(d)=1. Let pDistr(D) denote the set of parametric probability distributions over D with parameters in X. Instantiation u is well-formed for µ∈pDistr(D) iff 0 ≤ µ(d)[u] ≤ 1 and Σd∈D µ(d)[u] = 1.

Distancemeasure. The function d : U×U →R≥0 is a distance measure if for all u, u′ , u′′ ∈ U, it satisfies: (I) positiveness: d(u, u′ ) ≥ 0, (II) symmetry: d(u, u′ ) = d(u′ , u), and (III) triangle inequality: d(u, u′′ ) ≤ d(u, u′ ) + d(u′ , u′′ ).

3 Parametric Bayesian Networks

A parametric Bayesian network (pBN) is a BN in which a subset of entries in the conditional probability tables (CPTs) are polynomials over the parameters in X. Let parvi denote the set of parents for the node vi in the graph G.

Definition 1. The tuple B=(G,X,Θ) is a parametric Bayesian network (pBN) with directed acyclic graph G=(V,W) over random variables V={v1, . . . , vm}, set of parameters X={x1, . . . , xn}, and parametric CPTs Θ={Θvi | vi∈V } where Θvi : Eval(parvi )→pDistr(Dvi ).

The CPT row Θv(par) is a parametric distribution over Dv given the parent evaluation par. The CPT entry Θv(par)(d), short θ(v,d,par), is the probability that v=d given par. A pBN without parameters, i.e., X = ∅, is an ordinary BN. A pBN B defines the parametric distribution function PrB over Eval(V ). For well-formed instantiation u, BN B[u] = (G,Θ[u]) is obtained by replacing the parameter xi in the parametric functions in the CPTs of B by u(xi).

Example 1. Fig. 1a shows a pBN over variables V = {C, S,A, P} (initial letters of node names) and parameters X = {p, q}. Instantiating with u0(p) = 0.72 and u0(q) = 0.95 yields the BN B[u0] as indicated using dashed boxes.

pBN constraints. Constraints PrB(H|E) ∼ λ involve a hypothesis H, an evidence E, ∼∈ {≤, ≥} and threshold 0 ≤ λ ≤ 1. For the instantiation u : X → R,

Sensitivity functions. Sensitivity functions are rational functions that relate the result of a pBN query to the parameters [Castillo et al., 1997; Coup´e and van der Gaag, 2002].

Example 2. For the pBN B1 in Fig. 1a, let φ1 be the constraint: PrB1 (C = no | A = pos ∧ P = pos) ≤ 0.009. The left-hand side reduces to the sensitivity function

Using instantiation u with u(p) = 0.92075 and u(q) = 0.97475, fB1,φ1 [u] = 0.008798 ≤ 0.009, i.e., B1[u] |= φ1.

Higher degree sensitivity functions. Contrary to the existing literature that is often restricted to multi-linear sensitivity functions, we are interested in analyzing pBNs with sensitivity functions of higher degrees.

Example 3. Consider a variant of the COVID-19 example, where the person only takes the antigen test twice rather than taking both the antigen and PCR tests to diagnose COVID19; see Fig. 2 for the parametric CPTs. Then,

Parametrization. Let Θmodif denote the CPT entries in BN B = (G, Θ) that are explicitly changed.

Definition 2 (BN parametrization). Tuple B = (G,X,Θ′ ) is a parametrization of BN B = (G,Θ) over X w.r.t. Θmodif if

The parametrization B is monotone iff for each parametric entry f , f is a monotonic polynomial. The parametrization B is valid iff B is a pBN, i.e., iff Θ′v (par)∈pDistr(Dv) for each random variable v ∈V and its parents configuration par.

A valid monotone parametrization is obtained by the linear proportional [Coup´e and van der Gaag, 2002] that is the most established scheme in the literature to parametrize BNs.

Definition 3 (Linear proportional parametrization). A linear proportional parametrization over X maps the CPTs Θ onto the parametric CPTs Θ′ w.r.t. Θmodif if

Note that we allow the repetition of parameters in multiple distributions and multiple CPTs, see e.g., Fig. 2.

Formal problem statement. Consider BN B = (G,Θ), its valid parametrization B = (G,X,Θ′ ) over X with constraint φ. Let u0 be the original value of the parameters X in B and d : U × U → R≥0 a distance measure. Let ˆ d0 denote an upperbound for d{u∈U}(u0, u). The minimum-distance parameter tuning problem is to find umin = argmin{u∈U | B[u]|=φ} d(u, u0). Its generalized variant for 0 ≤ ϵ ≤ ˆ d0 and 0 ≤ η ≤ 1 is:

Definition 4. The (ϵ, η)-parameter tuning problem is to find u∈U s.t. B[u] |=φ, d(u, u0)≤ϵ, and d(u, umin)≤(1-η) · ˆ d0.

( ˆ d0, 1)-tuning gives the minimum-distance tuning. We call η the coverage factor that determines, intuitively speaking, the minimality of the result; we discuss this in Sec. 5.4. We consider two distance measures.

Euclidean distance (EC distance). The ECdistance between u and u′ (both in U) is defined as:

Corollary 1. For n = |X|, ˆ d0 = √n is an upperbound for the EC distance of any u0 from any instantiation u∈U.

Corollary 2. Let 0 ≤ ϵ ≤ √n and u0(xi) - ϵn ≤ u(xi) ≤ u0(xi) + ϵn . Then, EC(u, u0) ≤ ϵ.

Chan-Darwiche distance (CD distance). CD is a measure to quantify the distance between the probability distributions of two BNs, defined as [Chan and Darwiche, 2005]:

where both 0/0 = 1 and ∞/∞ = 1 by definition. Whereas the EC distance can be computed in O(n), this is for CD distance NP-complete [Kwisthout and van der Gaag, 2008]. It has known closed-forms only for limited cases, e.g., single parameter and single-CPT pBNs [Chan and Darwiche, 2004]. Let Θv be the CPT that is parametrized to Θ′v by a monotone parametrization. Let the minimum (maximum) probability entry of Θv be θmin (θmax) parametrized to θ′min with x (θ′max with y). Let lbx, ubx, lby, uby ̸= 0, 1 be the upperbounds and the lowerbounds of x and y in U.

Corollary 3. An upperbound for the CD distance of u0 from any instantiation u∈U is:

as derived from the single CPT closed-form [Chan, 2005].

Corollary 4. Let 0≤ϵ≤ ˆ d0 and α=eϵ/2 . Let for each θ ∈ Θv, α-1 · θ′ [u0] ≤ θ′ [u] ≤ α · θ′ [u0] . Then, CD(u, u0) ≤ ϵ.

We use the proportional parametrization (Def. 3) that has several useful properties: e.g., it is zero-preserving and minimizes the CD distance for single-CPT pBNs [Renooij, 2014].

4 Parametric Markov Chains

Parametric Markov chains are an extension of Markov chains (MCs) that allow parametric polynomials as transition labels: Definition 5. A parametric Markov chain (pMC) M is a tuple (S, sI,X,P) where S is a finite set of states, sI ∈ S is the initial state, X is a finite set of real-valued parameters, and P : S → pDistr(S) is a transition function over the states.

Example 4. Figures 1b and 1c are examples of pMCs.

Reachability constraints. A reachability constraint ϕ = (T, ∼, λ) is defined over pMC M = (S, sI , X, P) with targets T ⊆ S, threshold 0 ≤ λ ≤ 1, and ∼ ∈ {≤,≥}. For pMC M and instantiation u : X → R:

where PrM[u](♦T) denotes the probability that MC M[u] reaches a state in T .

Example 5. Consider the pMC MB,E in Fig. 1c. Let T = {s10}. Then, PrMB,E (♦s10) = 361 34900 · p · q + 8758 · q + 361 . Let R ⊆ Rn with n = |X| and let

We now define the parameter synthesis problems for pMCM and reachability constraint ϕ that are relevant to our setting.

Definition 6 (Region partitioning). Partition region R into R+, R-, and R? such that:

and R? = R \ (R+ ∪ R-) with ||R?|| ≤ (1-η)·||R|| for some given coverage factor 0 ≤ η ≤ 1.

The sub-region R? denotes the fragment of R that is inconclusive for ϕ. This fragment should cover at most fraction 1-η of R's volume. Region partitioning is exact if R? = ∅.

Definition 7 (Region verification). For the region R and the specification φ, the problem is to check whether:

4.1 Parameter Lifting

The parameter lifting algorithm (PLA) [Quatmann et al., 2016] is an abstraction technique that reduces region verification to a probabilistic model checking problem for which efficient algorithms exist [Katoen, 2016]. It first removes parameter dependencies by making all parameters unique. Then it considers for each parameter xi only its bounds lbxi and ubxi within the region R. This yields a non-parametric Markov Decision Process (MDP).

Definition 8. A Markov decision process (MDP) is a tuple M = (S, sI, Act,P) with a finite set S of states, an initial state sI ∈ S, a finite set Act of actions, and a (partial) transition probability function P : S × Act → Distr(S).

While resolving the non-determinism in the obtained MDP, a trade-off between maximizing or minimizing xi may occur. This occurs in particular when parameter xi repeats in the outgoing transitions of multiple states, see e.g., parameter t in Fig. 3a. PLA handles the issue by a relaxation step that introduces intermediate parameters; see e.g., Fig. 3b, 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.

Example 6. Consider the pMC in Fig. 3(a) and region t ∈ [0.0075, 0.0125]. Fig. 3b shows the pMC after relaxation, e.g., parameter t in the outgoing transitions of state C=no,A=pos is replaced by t′ . Fig. 3c 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 0.0075 and 1-0.0075 obtained from lb(t). The probabilities 0.0125 and 1-0.0125 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 3. The region accepts ϕ for the pMC if the resulting MDP satisfies ϕ for all schedulers. If all the schedulers satisfy ¬ϕ, 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 η% 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 [Salmani and Katoen, 2021a; Salmani and Katoen, 2021b]; see the next example.

Example 7. Consider the acyclic pMC in Fig. 1b for the pBN in Fig. 1a for the topological ordering ϱ : C < S < A < P. Initially, all variables are don't care. The initial state can evolve into C=yes and C=no. Its transition probabilities 0.05 and 0.95 come from the CPT entries of C. Generally, at "level" i of the pMC, the outgoing transitions are determined by the CPT of the i+1-st variable in ϱ.

pBN constraints relate to reachability constraints in its pMC.

Example 8. Consider the pBN B1 from Fig. 1a and the query PrB(C = no | A = pos ∧ P = pos) from Example 2. Let MϱB be the pMC in Fig. 1b. The query coincides with

5 Region-basedMinimal Change Tuning

Our approach to the parameter tuning problem for a pBN and constraint φ consists of two phases. The ϵ-close partitioning exploits parameter lifting on the pMC for reachability constraint ϕ of φ. We start with a rectangular region enclosing the original instantiation u0. This region is iteratively expanded if it rejects ϕ. Inconclusive regions are partitioned. Any iteration that finds some accepting subregion ends the first phase. The second phase (Alg. 2) extracts a minimal-change instantiation from the set of accepting sub-regions.

5.1 Obtaining Epsilon-Close Subregions

Alg. 1 is the main thread of our approach. Its inputs are the pBNB = (V,W,X,Θ), the initial instantiation u0 : X →R, and the constraint φ. It starts with obtaining the pMCMB of the pBN and reachability constraint ϕ of φ (lines 1-2). The output is either (i) u0 ifM[u0] |= ϕ, or (ii) the instantiation u with a minimal distance from u0 or (iii) no instantiation if ϕ is infeasible, that is, in the η = 100% of the checked parameter space, no satisfying instantiation has been found. The hyperparameters of the algorithm are the coverage factor 0 < η < 1, the region expansion factor 0 < γ < 1, and the maximum number of iterations K ∈ N. The hyper-parameters γ and K steer the iterative procedure. They initially determine the distance bound ϵ, see line 3. The bound ϵ determines region Rϵ, i.e., how far the region bounds deviate from u0. PLA then verifies Rϵ. If Rϵ is rejecting, ϵ is extended by the factor γ-1 (l. 9). Figure 4 visualizes the procedure for γ=1/2 and n=2. At iteration i, (a) PLA is invoked on Rϵ and either (b) R+ ̸= ∅ or (c) R+ = ∅. For case (b), the iterative procedure

Algorithm 1: Minimal change tuning

ends and Alg. 2 is called to obtain an instantiation in R+ that is closest to u0. For case (c), the region is expanded by factor γ-1 and passed to PLA, see (d). Note that the loop terminates when the distance bound ϵ reaches its upper bound ˆ d0 (l. 2). We refer to Sec. 3 Corollary 1 and 3 for computing ˆ d0.

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

5.2 Obtaining aMinimal-Distance Instantiation

Once R+ ̸= ∅ ⊆ Rϵ is found, it remains to find the instantiation u+ in R+ that is closest to u0. Region R+ includes infinitely many instantiations, yet minimal-distance instantiation is computable in O(k·n) for n = |X| and the regions R+ = {R+,1, · · · ,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. 2 simplifies finding u+. The idea (l. 4-14) is to pick u+,i fromR+,i that has minimal distance from u0. Let lbxj be the lower bound for parameter xj in the region R+,i and ubxj be its upper bound. We set the value for each parameter based on three possible cases. (1) If u0(xj)∈ [lbxj , ubxj ], we set u+,i(xj)=u0(xj). This yields the least distance, i.e., 0 in the dimension of xj, see Fig. 5

Algorithm 2: R+-minimal distance instantiation

(left). (2) The second case is u0(xj)<lbxj , see Fig. 5 (middle). Then the lower bound lbxj has the least distance from u0(xj) in the dimension xj , i.e., |lbxj -u0(xj)| < |a-u0(xj)| for every value a ∈ [lbxj , ubxj ] and similarly for CD distance, ln lbxj u0(xj) < ln a u0(xj) . In this case, we set u+,i(xj)=lbxj to minimize the distance in dimension xj. (3) By symmetry, for u0(xj) > ubxj , we set u+,i(xj) = ubxj ; see Fig. 5 (right). It remains to compute the distance for each candidate in U+ = {u+,1, · · · , u+,k} and pick u+ closest to u0 (l. 15).

5.3 Example

We apply the entire algorithm on the pBN in Fig. 1a, constraint φ = Pr(C=no | A=pos ∧ P=pos) ≤ 0.009, u0 : p → 0.72, q → 0.95, and η=0.99. The pMCMB,E is given in Fig. 1c and the reachability constraint is ϕ : Pr(♦s10) ≤ 0.009. We take ˆ d0=√n, γ=1/2, and ϵ0=1/16.

1. The initial region Rϵ is obtained by ϵ=1/8: [0.72-1/8√2, 0.72+1/8√2] × [0.95-1/8√2, 0.95+1/8√2].
2. Running PLA on Rϵ gives no sub-region accepting ϕ.
3. Expanding the region by the factor γ-1 = 2 yields: [0.72-1/4√2, 0.72+1/4√2] × [0.95-1/4√2, 0.95+1/4√2].
4. Running PLA gives 12 regions accepting ϕ, e.g., R+,1 : [0.92075, 0.960375] × [0.97475, 1] and R+,2 : [0.960375, 0.97028175] × [0.9431875, 0.9495].
5. Running Alg. 2 on R+,1 through R+,12 gives minimal distance candidates u+,1 through u+,12; e.g., u+,2(p)=lbp=0.960375 as u0(p) < lbp and u+,2(q)=ubq=0.9495 as ubq < u0(q); see Fig. 5.
6. Running Alg. 2 using the distance measure of choicehere: EC distance-for the candidates U+ returns u+ closest to u0: u+(p) = 0.92075 and u+(q) = 0.97475; distance 0.040913125.

5.4 Discussion

Let us shortly discuss the hyper-parameters: η, γ, and K. The hyper-parameter η is the coverage factor for the region partitioning. Our region partitioning algorithm, i.e., parameter lifting (line 17 Alg. 1) works based on this factor: the procedure stops when R? is at most (1 - η) · ||R|| of the entire region R, e.g., for η = 0.99, it continues until at least 99% of R is either rejecting or accepting and at most 1% is unknown. This relates η to the approximation bounds of our algorithm, see Problem statement, Def. 4. Intuitively speaking, the coverage factor η means that if the algorithm provides a minimal-distance value D(u, u0) = d, then with the probability 1 - η there may exist a smaller value distance than d that works too but has not been found. One can thus consider η as the confidence factor for the minimality of the results.

The hyper-parameter γ specifies the factor by which the region is expanded at each iteration, see Line 9, Alg. 1. The hyper-parameters γ and K specify the size of the initial region, see Lines 3 and 16 of Alg. 1. Our experiments (not detailed in the paper) with γ = {0.2, 0.5, 0.8} and K = {2, · · · , 12} reveal that γ = 1/2 and K = 6 gave the best balance. Large values for γ (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 γ (0.2) lead to large regions in the next iterations due to the expansion by γ-1 .

6 Experimental Evaluation

We empirically evaluated our approach using a prototypical realization on top of the probabilistic model checker Storm [Hensel et al., 2022] (version 1.7.0). As baseline, we used the latest version (10.4) of Bayesserver 2 , 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 bnlearn 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 3 to pick the CPT entries Θmodif most relevant to the constraint. To get well-formed BNs, we used the linear proportional parametrization, see Def. 3. 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 [Salmani and Katoen, 2020]. We took γ=1/2 and K=6 for our experiments, see Sec. 5.4. 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 bnlearn repository. For each BN, we parameterized one distribution in a single CPT to align with the restrictions from Bayesserver. Figure 6a 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 η = 1-10-i for i ∈ {1, · · · , 6}. The influence of η on the exactness of our results is addressed under RQ2.

Findings. Storm outperforms Bayesserver for several benchmarks (cancer, sachs, win95pts, and multiple instances of hailfinder), whereas Bayesserver is faster by about an order of magnitude for the other pBNs, such as hepar2. We noted that 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.

2 https://www.bayesserver.com/

3 http://reasoning.cs.ucla.edu/samiam/

RQ2: Sensitivity to the coverage factor η. For each pBN and constraint, we decreased (the refinement factor) 1-η in a step-wise manner by a factor 10-1 . To quantify the tightness of our results, we measured how our approximatelyclose instantiation, denoted uStorm, differs from the absolute minimum-distance from Bayesserver, denoted uBayes. Figure 6b (log-log scale) plots the tightness of the results |uStorm-uBayes| (y-axis) against the refinement factor (x-axis). Figures 6c and 6d indicate the number of iterations and the tuning time (log scale, seconds) for each refinement factor.

Findings. (I) Mostly, 1-η bounds the difference between our approximately-close solution and Bayesserver's solution. For e.g., η = 1-10-4 , the difference is at most 10-4 . Recall from Def. 6 that the value of 1-η bounds the size of the unknown regions in the region partitioning procedure. This explains why |uStorm - uBayes| is bounded by 1-η. (II) On increasing the coverage, the difference to the minimum distance rapidly decreases. (III) The computation time moderately increases on increasing the coverage, but the number of iterations was mostly unaffected. 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. (IV) The timing is not correlated to the number of iterations. This is because the ϵ-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 λ. We varied the constraint's threshold (λ) by steps of 1/20 for the benchmarks alarm hepar2, and hailfinder with n=1 and η=99.99999%. Figures 7a, 7b, and 7c 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. For small perturbations the original BNs (B[u0]) may already satisfy the constraints. In general, for the constraints close to the original probability PrB[u0],φ, the number of iterations is low, the distance is naturally small, and the minimal-change tuning is completed fast. This is because region refinement starts with small regions in the close vicinity of the original values of the parameters without the need to analyze larger regions. By strengthening the threshold, the possibly satisfying regions get further away from u0. Thus the distance, the number of iterations, and sometimes the tuning time grow. For very large perturbations the constraint may become infeasible with the confidence factor η = 100%. Similar findings are valid for n>1 parameters; see RQ4.

1

RQ4: Scaling the number of parameters. We took the win95pts and alarm 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 1(left) and (right) list the results for win95pts and alarm. We list for each pBN, the number of affected CPTs, the number of parameters, the threshold λ, and the coverage η. E.g., win95pts 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).

Findings. (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 ≥ 0.25 is unsatisfiable for win95pts 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., alarm with 8 parameters for the constraint ≤ 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 3: 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 win95pts and alarm: we parameterized the entries θ1 · · · θk ∈ Θmodif over the same parameter x ∈ X when the original values of θ1 · · · θk were the same in the original BN. The eleventh row in Table 1 (left) and the eighth and eleventh rows in Table 1 (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 win95pts with 8 parameters, threshold ≥ 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.

7 Epilogue

Related work. Kwisthout and van der Gaag [2008] studied the theoretical complexity of tuning problems. Renooij [2014] 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 [2015; 2017] that consider tuning heuristics for distance optimizations. Peng and Ding [2005] 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. [2013] exploits linear programming for BN parameter tuning. Yak- aboski and Santos [2018] consider a new distance measure for parameter learning and tuning of BNs. Leonelli [2019] considers nonlinear sensitivity functions, yet only for single parameter pBNs and BallesterRipoll and Leonelli [2022] efficiently compute the derivatives of sensitivity functions to select the most relevant parameters to a query, yet they limit to single parameter variation.

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 [Bolt and van der Gaag, 2017] and using monotonicity of parameters [Spel et al., 2019].

Acknowledgements

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.