Paper deep dive
Interpretable Causal Discovery via Causal-Effect Constraints
Cixuan Zhang, Guy Van den Broeck, Benjie Wang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/14/2026, 4:46:05 AM
Summary
The paper introduces a method for conditional causal discovery that addresses the computational challenge of inferring causal structures under rare-event constraints, such as large causal effects. By casting the problem as Bayesian inference and adapting adaptive multilevel splitting (AMS) techniques, the authors propose a framework that drives a particle population toward constrained regions in the joint graph-parameter space. This approach enables the estimation of posterior tail probabilities and provides pathway-level summaries to aid scientific exploration, validated on synthetic data and the Sachs protein dataset.
Entities (6)
Relation Signals (5)
Adaptive Multilevel Splitting → usedin → Conditional Causal Discovery
confidence 95% · Our approach integrates adaptive multilevel splitting (AMS) with an MCMC kernel over the joint graph–parameter space to address conditional causal discovery.
Sachs Protein Dataset → usedforvalidation → Conditional Causal Discovery
confidence 92% · we show in a case study on the Sachs protein dataset how our method can be used to aid scientific exploration
Causal-Effect Constraint → defines → Conditional Posterior
confidence 90% · we target the posterior over causal graphs and parameters conditional on an event such as a causal-effect constraint.
Linear Gaussian Structural Equation Model → implements → Causal Bayesian Networks
confidence 88% · In this paper, the generic mechanism parameters θ are instantiated by a linear-Gaussian structural equation model.
Bayesian causal discovery → handles → Epistemic Uncertainty
confidence 85% · We take a Bayesian approach to this problem in which we represent uncertainty over graph structure and parameters.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Causal discovery aims to uncover the underlying causal relationships given data generated from a system. The goal, however, is not merely to predict causal edges given data, but also to be able to interpret and explain either observed or hypothesized phenomena, such as a particularly large causal effect. We consider this task of conditional causal discovery and cast it as a Bayesian inference problem, in which we target the posterior over causal graphs and parameters conditional on an event such as a causal-effect constraint. Unfortunately, this poses a computational challenge: existing approaches to Bayesian causal discovery struggle when the event has small posterior mass. To address this, we adapt rare-event estimation techniques to perform inference the joint graph-parameter space. Our method gradually drives a particle population toward the constrained region while maintaining samples that approximate the conditional posterior. Empirical evaluation on synthetic graphs validates the accuracy of our approach at small and large scales, and we show in a case study on the Sachs protein dataset how our method can be used to aid scientific exploration by providing pathway-level summaries.
Tags
Links
- Source: https://arxiv.org/abs/2608.12640v1
- Canonical: https://arxiv.org/abs/2608.12640v1
Trouble viewing inline? Open PDF directly →
Full Text
79,962 characters extracted from source content.
Expand or collapse full text
Interpretable Causal Discovery via Causal-Effect Constraints Cixuan Zhang Affiliation: Computer Science Dept. Yale University New Haven, Connecticut, USA Guy Van den Broeck Affiliation: Computer Science Dept. University of California, Los Angeles Los Angeles, California, USA Benjie Wang Affiliation: Computer Science Dept. University of California, Los Angeles Los Angeles, California, USA Abstract Causal discovery aims to uncover the underlying causal relationships given data generated from a system. The goal, however, is not merely to predict causal edges given data, but also to be able to interpret and explain either observed or hypothesized phenomena, such as a particularly large causal effect. We consider this task of conditional causal discovery and cast it as a Bayesian inference problem, in which we target the posterior over causal graphs and parameters conditional on an event such as a causal-effect constraint. Unfortunately, this poses a computational challenge: existing approaches to Bayesian causal discovery struggle when the event has small posterior mass. To address this, we adapt rare-event estimation techniques to perform inference the joint graph-parameter space. Our method gradually drives a particle population toward the constrained region while maintaining samples that approximate the conditional posterior. Empirical evaluation on synthetic graphs validates the accuracy of our approach at small and large scales, and we show in a case study on the Sachs protein dataset how our method can be used to aid scientific exploration by providing pathway-level summaries. 00footnotetext: Our code is available at https://github.com/ZCX031116/MLS-Framework. 1 Introduction Causal discovery is often motivated not only by prediction, but by the need to obtain interpretable and actionable explanations of how a system works. In scientific domains, a practitioner may not simply ask which graph is most probable given the data; rather, they may ask which causal mechanisms could explain a particular domain-relevant phenomenon. For example, in the well-known Sachs study of protein interactions [31], one may want to understand which directed pathways could support an unusually large effect from one protein to another, and whether such pathways suggest plausible interventions or follow-up experiments. Posterior samplesp(G,B∣)p(G,B ): many plausible DAGsiijjiijjCausal-effect constraintℰij+(t)=(G,B):CEij(G,B)≥tE_ij^+(t)=\(G,B):CE_ij(G,B)≥ t\ttrare posterior eventConditional posteriorp(G,B∣,ℰij(t))p(G,B ,E_ij(t))iijjhigh-frequency pathwaysexplain the eventscore byCEijCE_ijMLS+ MCMC Figure 1: Overview of conditional causal discovery. Starting from the unconstrained Bayesian posterior over causal graphs and edge weights, we condition on a user-specified extreme causal-effect event and obtain a constrained posterior whose samples can be summarized at the pathway level. In this work, we formulate this type of “what-if” analysis as a conditional form of causal discovery. Given observational data, we seek causal structures that are both statistically plausible and consistent with a user-specified constraint, such as a large causal effect from node i to node j. This presupposes that there could be many graphs that explain the data well, and as such requires a treatment of epistemic uncertainty. We take a Bayesian approach to this problem in which we represent uncertainty over graph structure and parameters. Rather than committing to a single estimated graph, we consider the posterior distribution over directed acyclic graphs (G) and parameters (θ). Our target is the posterior conditioned on an event determined by the graph and parameters, such as the event that the signed causal effect from i to j exceeds a threshold t. The goal is to characterize the causal structures and directed pathways that remain plausible under the causal-effect constraint, while also estimating the posterior probability of the event itself. Conditional causal discovery poses a computational challenge that is not well handled by current tools. In particular, the event of interest may be rare under the unconstrained posterior, especially when the causal effect threshold t is large, or when we impose a combination of constraints. Standard posterior samplers based on MCMC may produce few or no samples satisfying the constraint, and simple rejection-based conditioning becomes inefficient. Even when constraint-satisfying samples are obtained, poor mixing can lead to unreliable statistical summaries of e.g. different causal paths. We therefore cast conditional causal discovery under strong (e.g. extreme-effect) constraints as a rare-event posterior inference problem. Our approach integrates adaptive multilevel splitting (AMS) [6] with an MCMC kernel over the joint graph–parameter space. Starting from posterior samples obtained from any Bayesian causal discovery method, our method progressively tightens the effect threshold while maintaining a representative particle population, producing both tail-probability estimates and conditional posterior samples. These samples can then be used to summarize which edges and directed pathways are most characteristic of the extreme-effect regime. Our contributions are as follows: • We formulate conditional causal discovery as posterior inference under extreme-effect constraints, targeting both constrained posterior samples and posterior tail probabilities. • We provide a practical rare-event inference procedure that can estimate small posterior tail probabilities while producing representative graph–parameter samples from the corresponding constrained posterior. • We validate the method on linear-Gaussian benchmarks with (d∈4,8,16,32d∈4,8,16,32) and present a Sachs dataset case study showing how extreme-effect conditioning enables pathway-level interpretation and hypothesis generation. 2 Related Work Causal discovery under background knowledge or structural constraints has been widely studied. Such constraints often encode qualitative statements about graph structure, including required or forbidden edges, paths, or ancestral relations [27, 4, 1]. Prior work has considered both identifiability and algorithmic procedures for learning graphs subject to these constraints [7, 8]. A related line of work uses interventional data as an additional source of information beyond pure observations [19, 5]. Instead, we are primarily motivated by interpretability rather than encoding fixed knowledge. Our technical approach differs in that the constraint is encoded flexibly as a quantitative score function, and in particular we encode constraints corresponding to the value of a causal-effect functional rather than a local structural statement. Moreover, rather than seeking exact identification of a single graph, we target a conditional Bayesian posterior that retains uncertainty over both graphical and parametric structure. Another line of work combines causal discovery with extreme value theory to infer causal directions when causal mechanisms are most visible in the tails of the observed distribution [17, 28, 3]. In these approaches, extremes are primarily a feature of the data-generating distribution and are used for identifiability or structure recovery. Our motivation is different: in our setting, the extreme event is a user-specified constraint on a causal-effect functional, and the inferential target is the posterior distribution conditioned on that constraint. Thus, our goal is not to identify causal directions from heavy-tailed observations, but to characterize which posterior graphs, parameters, and pathways explain an unusually large or small causal effect. For posterior inference over causal graphs, a common approach is to use Markov chain Monte Carlo (MCMC) sampling over DAGs or higher-level representations such as orders [12, 23, 36, 15]. More recent work has developed efficient DAG-space MCMC samplers with locally informed and adaptive proposals, including PARNI-DAG [24], which constructs adaptive random neighborhoods guided by posterior information and can exploit a pre-tuned skeleton to improve scalability. These graph posterior samples are useful not only for representing structural uncertainty, but also for downstream causal inference; recent benchmarking work evaluates Bayesian causal discovery methods through downstream treatment-effect estimation [11]. Alternative approaches approximate the posterior via variational inference over graphs [2, 25, 9, 37, 10, 30, 34]. These methods aim either to generate graph samples reflecting structural uncertainty or to support Bayesian model averaging for causal inference [35]. Our computational approach depends on the rare-event simulation literature, particularly adaptive multilevel splitting [6]. Splitting methods estimate small probabilities by introducing intermediate thresholds and repeatedly propagating a particle population toward rarer events. We adapt this idea to Bayesian causal discovery by defining levels through a causal-effect score and using MCMC moves over the joint graph–parameter space. This produces both an estimate of the posterior probability of the extreme-effect event and samples from the corresponding conditional posterior. 3 Preliminaries Causal Bayesian Networks A Bayesian network (BN) (G,θ)(G,θ) is a probabilistic model p()p( X) over d variables =X1,…,Xd X=X_1,…,X_d, specified by a directed acyclic graph (DAG) G and mechanism parameters θ. The graph encodes conditional independencies, while θi _i parameterizes the conditional distribution of XiX_i given its parents. The joint distribution factorizes as p(∣G,θ)=∏i=1dp(Xi∣paG(Xi),θi),p( X G,θ)= _i=1^dp (X_i _G(X_i), _i ), where paG(Xi)pa_G(X_i) denotes the parents of XiX_i in G. In this paper, the generic mechanism parameters θ are instantiated by a linear-Gaussian structural equation model. Specifically, the observed variables satisfy =B+ϵ X= XB+ ε, where B∈ℝd×dB ^d× d is the weighted adjacency matrix and ϵ∼(,Σ) ε ( b, ), with ∈ℝd b ^d and diagonal Σ∈ℝ≥0d×d _≥ 0^d× d. For a given DAG G, we impose Bij=0B_ij=0 whenever i is not a parent of j in G. Causal Bayesian networks [33, 29] add a causal interpretation to the directed edges in G: they describe how the joint distribution changes under interventions. In the linear-Gaussian SEM, the signed total causal effect from XiX_i to XjX_j is the derivative of the post-intervention mean of XjX_j with respect to an intervention on XiX_i. It has the closed form CEij(G,B)≜CE(i→j∣G,B)=[(I−B)−1]ijCE_ij(G,B) (i→ j G,B)= [(I-B)^-1 ]_ij (1) following the standard total-effect formula for linear structural equation models [32]. Since G is acyclic, B is nilpotent after a topological ordering, so (I−B)−1=I+B+B2+⋯+Bd−1(I-B)^-1=I+B+B^2+·s+B^d-1. Thus, the (i,j)(i,j) entry of (I−B)−1(I-B)^-1 aggregates the products of edge weights along all directed paths from XiX_i to XjX_j. This path-sum interpretation is useful below because our conditional posterior summaries focus on which edges and directed pathways explain unusually large or small total effects. Throughout the rest of the paper, we write a graph–weight state as Z=(G,B)Z=(G,B), and abbreviate the signed causal effect as CEij(Z)=CEij(G,B)CE_ij(Z)=CE_ij(G,B). For a target ordered pair (i,j)(i,j) and threshold t>0t>0, we use the phrase “extreme causal effect” to mean that the total causal effect falls in a user-specified posterior tail region. In particular, we define the right- and left-tail events ℰij+(t) _ij^+(t) =Z:CEij(Z)≥t, = \Z:CE_ij(Z)≥ t \, (2) ℰij−(t) _ij^-(t) =Z:CEij(Z)≤−t. = \Z:CE_ij(Z)≤-t \. When the sign is clear from context, we write ℰij(t)E_ij(t) for either tail event. These events play two roles in our inference problem. First, we estimate their posterior probability, such as ℙ(ℰij±(t)∣)P(E_ij^±(t) ). Second, we use them to define the constrained posterior p(G,B∣,ℰij±(t))p (G,B ,E_ij^±(t) ), from which we draw graph–weight samples for pathway-level summaries. Bayesian Causal Discovery Causal discovery [22, 16] is the problem of inferring the DAG G responsible for generating an observed dataset D. We make the common assumption of causal sufficiency, meaning that there are no latent confounders. Even under this assumption, a single DAG may not be reliably identifiable from finite observational data due to sampling uncertainty and Markov equivalence. Bayesian causal discovery therefore represents uncertainty through a posterior distribution over graphs and parameters rather than committing to a single structure. We place a user-specified prior p(G)p(G) over DAGs and use the BGe marginal likelihood p(∣G)p(D G) for linear-Gaussian models [13, 14]. The method itself does not require a particular graph prior; in the experiments, we use a sparse Erdős–R’enyi DAG prior, with the exact sparsity settings reported in Appendix A.1. Given G, the posterior over edge weights factorizes by node: p(B∣G,)=∏j=1dp(BpaG(Xj),j∣G,),p(B G,D)= _j=1^dp (B_pa_G(X_j),j G,D ), where coefficients outside the parent set of XjX_j are fixed to zero and each nonzero incoming-coefficient block follows a multivariate t-distribution [36]. This node-wise factorization is used later by the MCMC mutation kernel: when a graph proposal changes a node’s parent set, only the affected incoming-coefficient blocks need to be refreshed from their conditional posterior. The joint posterior over graph–weight states is p(G,B∣)∝p(G),p(∣G),p(B∣G,).p(G,B ) p(G),p(D G),p(B G,D). (3) We write π(Z)=π(G,B)=p(G,B∣)π(Z)=π(G,B)=p(G,B ) for this unconstrained posterior in the method section for compactness. 4 Conditional Causal Discovery In this section, we formulate conditional causal discovery as posterior inference under user-specified causal-effect constraints, and then describe how we sample from the resulting conditional posterior. The main challenge is that the constraint may have low probability under the posterior. If ordinary posterior samples almost never satisfy the event, then rejection sampling gives both unstable probability estimates and too few constrained samples for pathway-level interpretation. Our solution is to convert each query into a scalar score h(Z)h(Z), where larger values indicate greater progress toward the desired event, and then use adaptive multilevel splitting to reach the event through a sequence of easier conditional problems. 4.1 Conditional causal discovery as a score-level problem Let Z=(G,B)Z=(G,B) denote a graph–weight state and let π(Z)=p(G,B∣)π(Z)=p(G,B ) be the unconstrained posterior from Eq. equation 3. Throughout this section, we use π()π(A) to denote the posterior mass of an event A. A conditional query is specified by a score function h:→ℝh:Z and a final level λ⋆ _ . The score-level event, its posterior probability, and the corresponding constrained posterior are λ⋆ _ _ =Z:h(Z)≥λ⋆, =\Z:h(Z)≥ _ \, (4) p⋆ p_ =π(λ⋆), =π(A_ _ ), π⋆(Z) _ (Z) =π(Z∣λ⋆). =π\! (Z _ _ ). This formulation separates the scientific query from the computational procedure: once the query has been written through h and λ⋆ _ , the same sampler can be applied to single effects, left-tail effects, or multiple simultaneous constraints. For a single ordered pair (i,j)(i,j) and threshold t>0t>0, the right-tail event ℰij+(t)E_ij^+(t) is represented by choosing h(Z)=CEij(Z)h(Z)=CE_ij(Z) and λ⋆=t _ =t. The left-tail event ℰij−(t)E_ij^-(t) is represented by choosing h(Z)=−CEij(Z)h(Z)=-CE_ij(Z) and λ⋆=t _ =t. Thus both tails are written as the same score-level event λ⋆=Z:h(Z)≥λ⋆A_ _ =\Z:h(Z)≥ _ \, which avoids requiring separate algorithms for positive and negative extreme effects. We also support conjunctions of multiple causal-effect constraints. Let =(ic,jc,⋈c,tc,κc)c=1CC=\(i_c,j_c, _c,t_c, _c)\_c=1^C be a collection of C constraints, where ⋈c∈≥,≤ _c∈\≥,≤\ gives the inequality direction, tc∈ℝt_c is the threshold, and κc>0 _c>0 is an optional scale factor. Define δc=+1 _c=+1 when ⋈c=≥ _c=≥ and δc=−1 _c=-1 when ⋈c=≤ _c=≤. The normalized margin of constraint c and the aggregate score are mc(Z) m_c(Z) =δc(CEicjc(Z)−tc)κc, = _c (CE_i_cj_c(Z)-t_c ) _c, (5) h(Z) h_C(Z) =min1≤c≤Cmc(Z). = _1≤ c≤ Cm_c(Z). The margin mc(Z)m_c(Z) is nonnegative exactly when the c-th constraint is satisfied. Therefore h(Z)≥0h_C(Z)≥ 0 means that all constraints are satisfied, so the joint event is ℰ=Z:h(Z)≥0E_C=\Z:h_C(Z)≥ 0\. Interval constraints can be represented by including both a lower-bound and an upper-bound inequality. 4.2 Adaptive multilevel splitting for rare posterior events Extreme-effect events can have very small posterior probability. A direct posterior sampler with M samples will produce only about Mp⋆Mp_ constrained samples on average, which may be close to zero if p⋆p_ is small. Adaptive multilevel splitting addresses this by replacing one difficult rare-event problem with a sequence of easier conditional problems. Intuitively, particles are first asked to reach a moderately high score level, then a higher one, and so on until they reach the target level λ⋆ _ . At each stage, particles that have made sufficient progress are retained and resampled, while MCMC mutation restores diversity within the current truncated posterior. Let λ0<λ1<⋯<λK=λ⋆ _0< _1<·s< _K= _ be increasing score levels, with λ0=−∞ _0=-∞ corresponding to the unconstrained posterior. Define k=Z:h(Z)≥λkA_k=\Z:h(Z)≥ _k\. Because these events are nested, the rare-event probability decomposes as p⋆=π(K)=∏k=0K−1π(k+1∣k).p_ =π(A_K)= _k=0^K-1π(A_k+1 _k). (6) The advantage is that each factor in this product can be much larger than p⋆p_ itself, making it estimable with a moderate number of particles. We use the adaptive version of multilevel splitting [20, 18, 6], which chooses intermediate levels from the particle population rather than requiring them to be fixed in advance. At level k, the particle population k=Zk,nn=1NP_k=\Z_k,n\_n=1^N is first mutated so that it approximately follows the level-truncated posterior πλk(Z)∝π(Z)h(Z)≥λk. _ _k(Z) π(Z)1\h(Z)≥ _k\. (7) We then evaluate the scores h(Zk,n)h(Z_k,n) and choose the next level as an empirical quantile: λ~k+1 λ_k+1 =Q1−ρ(h(Zk,n)n=1N), =Q_1-ρ\! (\h(Z_k,n)\_n=1^N ), (8) λk+1 _k+1 =minλ~k+1,λ⋆. = \ λ_k+1, _ \. Here ρ∈(0,1)ρ∈(0,1) is the survival fraction and Q1−ρQ_1-ρ is the empirical (1−ρ)(1-ρ)-quantile. The estimated conditional factor at this level is the fraction of particles that survive the new threshold, β^k+1=1N∑n=1Nh(Zk,n)≥λk+1. β_k+1= 1N _n=1^N1\h(Z_k,n)≥ _k+1\. (9) The survivors are resampled with replacement to form the starting population for the next level. The final probability estimate is p^⋆=∏k=0K−1β^k+1 p_ = _k=0^K-1 β_k+1, and after the final mutation step the particles approximate samples from π⋆(Z) _ (Z). Algorithm 1 Adaptive multilevel splitting for conditional causal discovery 1: Posterior π(Z)π(Z), score h, target level λ⋆ _ , initial particles P, particle size N, survival fraction ρ, mutation steps m, maximum levels KmaxK_ . 2: Tail-probability estimate p^⋆ p_ and constrained particles P. 3: Initialize =Znn=1NP=\Z_n\_n=1^N using an approximate posterior sampler or warm-start routine targeting π(Z)π(Z). 4: Set λ0←−∞ _0←-∞ and p^⋆←1 p_ ← 1. 5: for k=0,…,Kmax−1k=0,…,K_ -1 do 6: ←Mutateλk(,m)P _ _k(P;m). 7: Compute scores sn←h(Zn)s_n← h(Z_n) for all Zn∈Z_n . 8: Set λk+1←minQ1−ρ(s1,…,sN),λ⋆ _k+1← \Q_1-ρ(s_1,…,s_N), _ \. 9: Let ←Zn∈:sn≥λk+1S←\Z_n :s_n≥ _k+1\. 10: Estimate β^k+1←||/N β_k+1←|S|/N and update p^⋆←p^⋆β^k+1 p_ ← p_ β_k+1. 11: if ||=0|S|=0 then 12: return 0,∅0, . 13: end if 14: Resample N particles from S with replacement to form P. 15: if λk+1=λ⋆ _k+1= _ then 16: ←Mutateλ⋆(,m)P _ _ (P;m). 17: return p^⋆, p_ ,P. 18: end if 19: end for 20: return p^⋆, p_ ,P. 21: 22: Here Mutateλ(,m)Mutate_λ(P;m) denotes m Metropolis–Hastings 23: steps per particle targeting πλ(Z)∝π(Z)h(Z)≥λ _λ(Z) π(Z)1\h(Z)≥λ\. 4.3 Particle initialization The splitting procedure should be viewed as a rare-event wrapper around Bayesian causal discovery rather than as a replacement for it. In principle, the initial population can be obtained from any posterior sampler that approximately targets π(Z)π(Z), including a long MCMC chain, an informed DAG-space sampler, or another Bayesian structure-learning method. The role of multilevel splitting is then to take these ordinary posterior samples and concentrate computation on the score-level region λ⋆A_ _ . In our implementation, we use a simple warm-start scheme. We first draw sparse DAGs G0,nG_0,n from the structural prior p(G)p(G), using an Erdős–Rényi DAG prior in the experiments with edge inclusion probability pedgep_edge chosen to match a target expected number of edges per node. Given each initial graph, we draw weights from the conjugate posterior B0,n∼p(B∣G0,n,)B_0,n p(B G_0,n,D) and set Z0,n=(G0,n,B0,n)Z_0,n=(G_0,n,B_0,n). Because these raw initial particles are not assumed to be exact samples from π(Z)π(Z), we use the inner MCMC kernel at the initial level λ0=−∞ _0=-∞ as a posterior warm-up before adapting the first nontrivial threshold. This makes the first population used in the splitting product approximate the unconstrained posterior, while retaining a simple and scalable initialization procedure. 4.4 Inner Metropolis–Hastings over joint (G,B)(G,B) At each score level λk _k, mutation must preserve the level-truncated target πλk(Z) _ _k(Z) from Eq. equation 7. We use a blocked Metropolis–Hastings kernel over the joint state Z=(G,B)Z=(G,B), alternating between graph updates and coefficient updates. This is important because the rare-event constraint depends on both structure and weights: changing only G can leave the sampler stuck at a fixed set of coefficients, while changing only B cannot explore alternative pathways. With probability pstructp_struct, we propose a structure move. The new graph G′G is drawn from a proposal qG(⋅∣G)q_G(· G), instantiated either as Structure-MCMC, which adds, deletes, or reverses a single edge while rejecting cyclic graphs, or as PARNI-DAG [24], which uses locally informed adaptive neighborhoods guided by posterior edge information. After proposing G′G , we refresh only the coefficient blocks whose parent sets changed. For each affected node v, we sample BpaG′(v),v′∼p(BpaG′(v),v∣G′,),B _pa_G (v),v p\! (B_pa_G (v),v G ,D ), (10) set coefficients for absent edges to zero, and copy all unchanged blocks from B. With probability 1−pstruct1-p_struct, we propose a weight move at fixed structure. We keep G′=G =G, select a node v with at least one parent, and resample its incoming coefficient block, BpaG(v),v′∼p(BpaG(v),v∣G,),B _pa_G(v),v p\! (B_pa_G(v),v G,D ), (11) leaving all other entries unchanged. These blocked refreshes are inexpensive under the conjugate linear-Gaussian model because each conditional coefficient posterior is a multivariate t distribution. Let q(Z′∣Z)q(Z Z) denote the complete proposal density, including the selected move type, the graph proposal when applicable, and the coefficient-refresh density. Since the current state already satisfies h(Z)≥λkh(Z)≥ _k, the level-λk _k acceptance probability is aλk(Z,Z′) a_ _k(Z,Z ) =h(Z′)≥λk =1\h(Z )≥ _k\ (12) ×min1,π(Z′)q(Z∣Z′)π(Z)q(Z′∣Z). × \! \1, π(Z )q(Z Z )π(Z)q(Z Z) \. Thus, proposals that violate the current score-level constraint are rejected immediately, while feasible proposals are accepted according to the usual Metropolis–Hastings ratio for the unconstrained posterior and the proposal probabilities. Repeating these moves after each resampling step helps remove duplicate particles and produces a more representative approximation to the constrained posterior. Zs 5 Experiments 5.1 Experiment Setup We evaluate the proposed conditional causal discovery procedure on synthetic linear-Gaussian datasets generated from random Erdős–Rényi DAGs with d∈4,8,16,32d∈\4,8,16,32\. For each dimension, we first sample a data-generating DAG G⋆G using an Erdős–Rényi DAG generator with target sparsity approximately 2d2d edges. The generator samples an acyclic ordering, draws only order-compatible directed edges, and returns an adjacency matrix G⋆G . The exact edge-sampling rule, including the edge probability or edge budget used for each d and the small-graph handling for d=4d=4, is provided in Appendix A.1. Given G⋆G , we draw raw edge weights independently as B~ij∼(0,1) B_ij (0,1) and mask them by the sampled adjacency, so that B⋆=B~⊙G⋆B = B G . Equivalently, Bij⋆=B~ijB _ij= B_ij if Gij⋆=1G _ij=1 and Bij⋆=0B _ij=0 otherwise. We then generate an observational dataset D of size nobs=1000n_obs=1000 from the linear-Gaussian SEM specified by (G⋆,B⋆)(G ,B ). We evaluate our multilevel splitting framework instantiated with either the Structure-MCMC kernel or the PARNI-DAG kernel [24]. We compare against four baselines: (i) exhaustive enumeration, which is a gold-standard baseline feasible for d=4d=4; (i) DiBS [25]; (i) OrderSPN [37]; and (iv) long single-chain MCMC using the same structural kernels as the multilevel splitting framework. Detailed sampler hyperparameters, thresholds, and evaluation settings are provided in Appendix A.1. The experiments are organized around three questions: small-graph accuracy, multi-effect conditioning behavior, and scalability beyond the enumerable setting. We describe each question in the corresponding subsection below. 5.2 Correctness of Conditional Causal Discovery The correctness experiments are designed to validate two complementary aspects of the method. The single-effect experiment checks the numerical accuracy of posterior tail-probability estimation against an exact enumerative reference. The multi-effect experiment checks the conditioning operator itself: when the conditioning event is constructed around ground-truth causal effects, the conditional posterior should assign more mass to graphs and weights close to the data-generating mechanism, and this concentration should increase as the constraint set becomes more informative. 5.2.1 Single-effect conditioning: validation on d=4d=4 We construct 10 test cases by sampling 10 independent d=4d=4 graphs and selecting one ordered node pair from each graph. For each pair, we evaluate a one-sided causal-effect tail query, represented in the score-level form λ⋆=Z:h(Z)≥λ⋆A_ _ =\Z:h(Z)≥ _ \ from Sec. 4. For every test case, we run each stochastic method 10 times and report the mean estimate together with across-run variability. Figure 2: Single-effect one-sided tail-probability estimation on d=4d=4 across 10 test cases, with exhaustive enumeration as the reference posterior probability. Error bars show across-run variability over 10 independent runs. DiBS is omitted from the plot because its estimates are orders of magnitude larger on these tail events, which would obscure the comparison among the remaining methods; the corresponding numerical values are reported in Table 3. Figure 2 compares all methods against exhaustive enumeration. We make two observations. First, all methods except DiBS are broadly consistent with the enumeration baseline across the d=4d=4 cases, with only minor deviations; the corresponding numerical values are reported in Table 3. Second, DiBS overestimates tail probabilities by orders of magnitude on these tail events. We therefore exclude DiBS from subsequent experiments and focus on methods that give reliable small-graph estimates. 5.2.2 Multi-effect conditioning: validating the conditioning operator on d=4d=4 and d=8d=8 The previous experiment validates tail-probability estimation for a single causal-effect event. We next test whether the method correctly conditions the joint posterior over graphs and weights. To make this test interpretable, we construct constraints from the known data-generating state and ask whether the resulting conditional posterior concentrates toward that state. Multi-constraint conditioning event. We select C ordered node pairs (ic,jc)c=1C\(i_c,j_c)\_c=1^C. Let ec⋆=CEicjc(Z⋆)e_c =CE_i_cj_c(Z ) denote the ground-truth causal effect for pair (ic,jc)(i_c,j_c) under the data-generating state Z⋆=(G⋆,B⋆)Z =(G ,B ). For each pair, we impose a two-sided interval constraint around ec⋆e_c : CEicjc(Z)∈[ec⋆−ε,ec⋆+ε],c=1,…,C.CE_i_cj_c(Z)∈[e_c - ,\,e_c + ], c=1,…,C. (13) Each interval is represented by two one-sided inequalities, and the conjunction of all inequalities defines the conditioning event. Following Sec. 4, we aggregate these inequalities using the minimum normalized margin score h(Z)h_C(Z), so that the target event is ℰ=Z:h(Z)≥0E_C=\Z:h_C(Z)≥ 0\. We compare a weak constraint set, with fewer constrained pairs and hence a looser event, against a strong constraint set, with more constrained pairs and hence a tighter event. All posterior summaries are computed from the final-stage multilevel-splitting particles and therefore condition on the target event by construction. To visualize the conditional posterior, we aggregate samples across runs to compute edge-frequency heatmaps, corresponding to posterior marginal edge probabilities, and average edge-weight heatmaps, corresponding to posterior mean weights. We also quantify structural recovery using the Structural Hamming Distance (SHD) between sampled DAGs and the ground-truth graph G⋆G . Table 1: d=4d=4 (PARNI) structural accuracy under multi-effect conditioning. Constraints Runs Samples/run Mean SHD ↓ Pr(SHD=0)↑ (SHD=0) Weak 10 20 0.325±0.1640.325± 0.164 0.790±0.1170.790± 0.117 Strong 10 20 0.115±0.0530.115± 0.053 0.885±0.0530.885± 0.053 Table 2: d=8d=8 (PARNI) structural accuracy under multi-effect conditioning. We report mean ± standard deviation of the per-run mean SHD, averaged over sampled DAGs within each run, and the fraction of samples with SHD=0SHD=0. Constraints Runs Samples/run Mean SHD ↓ Pr(SHD=0)↑ (SHD=0) Weak 10 50 5.342±0.4425.342± 0.442 0.040±0.0270.040± 0.027 Strong 10 50 2.290±0.1992.290± 0.199 0.108±0.0660.108± 0.066 Multi-effect conditioning results. As demonstrated by Figs. 6 and 7 in Appendix A.7, across both d=4d=4 and d=8d=8, stronger, more informative multi-effect constraints consistently concentrate the conditional posterior toward the ground-truth mechanism. Quantitatively, Tables 1 and 2 show lower SHD under strong constraints than under weak constraints, together with a higher fraction of exact structural matches. Qualitatively, the corresponding edge-frequency and edge-weight heatmaps show the same pattern: posterior mass becomes sharper and closer to the ground truth as the constraint set becomes stronger. Together, these SHD improvements and posterior-summary concentration patterns provide evidence that the framework correctly conditions on multiple causal-effect constraints in the joint graph–parameter space. 5.3 Scalability The scalability experiments examine whether the rare-event estimator remains stable when exact enumeration is no longer available. We focus on one-sided single-effect tail-probability estimation on d∈8,16,32d∈\8,16,32\ using multilevel splitting with the Structure-MCMC and PARNI-DAG kernels. These experiments test three aspects of scalability: whether estimated tail probabilities decrease smoothly as the target level λ⋆ _ becomes more stringent; whether independent runs give reproducible tail curves; and whether informed structure proposals become more important as the graph dimension increases. Figure 3: Single-effect one-sided tail-probability estimation on d∈8,16,32d∈\8,16,32\. Figure 3 summarizes one-sided single-effect tail-probability estimation on d∈8,16,32d∈\8,16,32\. On d=8d=8, we additionally include unconditional baselines, OrderSPN and two single-chain samplers, for comparison. Both multilevel-splitting variants produce smooth, monotone-decaying tail curves as the target level λ⋆ _ becomes more stringent. In contrast, the unconditional methods quickly fail to generate samples satisfying the more extreme targets, leading to degenerate probability estimates in the rare-event region. This highlights the practical advantage of multilevel splitting for conditional causal discovery on larger graphs, where rare-event conditioning makes unconditional sampling increasingly inefficient. The same curves also provide a threshold-sensitivity check. The final threshold λ⋆ _ defines the scientific query, so the conditional posterior should change as the threshold changes. Numerically, however, the estimated posterior mass should vary smoothly and reproducibly as the threshold becomes more extreme. Across independent runs, both structure kernels yield stable and monotone tail curves in the moderate dimensions, supporting this expected behavior beyond the enumerable d=4d=4 setting. Accordingly, for d=16d=16 and d=32d=32, we primarily assess robustness via internal consistency, because exact gold-standard posterior enumeration is unavailable at these scales. For d=16d=16, both multilevel-splitting instantiations produce stable and consistent probability estimates across runs, indicating that the framework remains well behaved as graph size grows. At d=32d=32, the two instantiations separate more clearly: PARNI-DAG retains relatively stable run-to-run behavior, whereas Structure-MCMC exhibits noticeably larger variability in the extreme tail. This suggests that reliable deep-tail estimation in larger graphs benefits from more informed structure proposals, as well as appropriate hyperparameter settings. 6 Case Study We now study an application of our conditional causal discovery framework to a real world protein-signalling dataset [31] commonly used in causal discovery benchmarks. In particular, we will ask the following question: Which causal structures and directed pathways most plausibly realize a specified extreme-effect event? Rather than claiming new biological insights, our goal is to show that conditional sampling can identify and summarize plausible mechanisms in an extreme-effect region that is rare under the unconditional posterior. Data, reference network, and target pairs. We use the Sachs protein-signaling dataset, containing n=7466n=7466 measurements of d=11d=11 proteins [31]. As a qualitative reference, we use the 20-edge consensus network in Fig. 4 [21, 31]. We study two target pairs: (i) PIP3→PIP2PIP3 2 (indices 6→56→ 5), which admits a clean SEM decomposition into a direct route PIP3→PIP2PIP3 2 and a mediated route PIP3→Plcg→PIP2PIP3 2; (i) Erk→AktErk (indices 1→101→ 10), a widely discussed cross-module influence in the Sachs signaling system. 11 1 Node order used throughout: [Raf, Erk, Plcg, PKC, PKA, PIP2, PIP3, Mek, P38, Jnk, Akt]. In particular, PIP3PIP3 is index 6, PIP2PIP2 is index 5, PlcgPlcg is index 2, ErkErk is index 1, and AktAkt is index 10. A published SEM mediation example on the same dataset reports a total effect of approximately 0.56650.5665 for PIP3→PIP2PIP3 2, with an explicit direct/indirect decomposition via PlcgPlcg, which we use as an external numerical anchor for the effect scale [26]. We treat such literature values as sanity checks on sign and magnitude rather than strict targets, since estimators and preprocessing differ. Conditioning events. Using the notation from Sec. 3, let Z=(G,B)Z=(G,B) denote a graph–weight state sampled from the posterior π(Z)=p(G,B∣)π(Z)=p(G,B ), and let CEij(Z)CE_ij(Z) denote the signed linear-SEM total effect. For the two target pairs, define the events ℰPIP(t)=Z:CEPIP3,PIP2(Z)≥tE_PIP(t)= \Z:CE_PIP3,PIP2(Z)≥ t \, andℰERK(t)=Z:CEErk,Akt(Z)≥tE_ERK(t)= \Z:CE_Erk,Akt(Z)≥ t \. We compare four posterior conditions: • Unconditioned: Z∼π(Z)=p(G,B∣)Z π(Z)=p(G,B ). • Cond-PIP: Z∼p(Z∣,ℰPIP(tPIP))Z p(Z ,E_PIP(t_PIP)), with tPIP=0.77t_PIP=0.77. • Cond-ERK: Z∼p(Z∣,ℰERK(tERK))Z p(Z ,E_ERK(t_ERK)), with tERK=0.67t_ERK=0.67. • Cond-Joint: Z∼p(Z∣,ℰPIP(tPIP′)∩ℰERK(tERK′)),Z p\! (Z ,\,E_PIP(t _PIP) _ERK(t _ERK) ), with (tPIP′,tERK′)=(0.74,0.65)(t _PIP,t _ERK)=(0.74,0.65), targeting a similarly rare posterior region. The estimated posterior masses of these events are reported in Appendix A.6 and Table 6. These conditioning events do not assert that the dataset corresponds to a single “ground-truth” extreme state. Instead, they define posterior queries of the form p(Z∣,ℰ)p(Z ,E), asking which graph–weight states and pathways remain plausible under the observed data when a specified causal effect, or a pair of causal effects, is unusually large. Why unconditioned summaries can disagree with the consensus graph. We run experiments on each condition with settings listed in Sec. A.1. Table 7 shows substantial unconditioned mass on zero total effect, corresponding to sampled DAGs with no directed path from the source to the target, even though the consensus network suggests nonzero coupling. This is expected: the consensus graph is a qualitative reference rather than a uniquely identified ground truth, and the pooled perturbation data do not uniquely determine reachability under our model, so many near-equivalent posterior graphs omit these paths [31, 21, 12]. Mechanistic analysis under conditioning. • Cond-PIP: The conditional posterior yields a highly concentrated explanation for PIP3→PIP2PIP3 2: the direct route and the PlcgPlcg-mediated route appear in essentially all conditional samples, and the conditional mean effect is approximately 0.7760.776 (calculated in Sec. A.6), with an increased PlcgPlcg-mediated share (Table 8). • Cond-ERK: Erk→AktErk remains direct-dominated, while PIP3→PIP2PIP3 2 becomes almost always reachable but is supported by many alternative, partly cancelling paths. This contrasts with the two-path concentration under Cond-PIP (Tables 7 and 8). • Non-symmetry: Conditioning on PIP3→PIP2PIP3 2 does not force Erk→AktErk to become extreme; for example, Pr(CEErk,Akt(Z)>0∣,ℰPIP(tPIP))≈0.498 \! (CE_Erk,Akt(Z)>0 ,E_PIP(t_PIP) )≈ 0.498. • Cond-Joint: Joint conditioning increases the frequency of the mediated Erk→Plcg→AktErk route relative to Cond-ERK, highlighting PlcgPlcg as a shared mediator under co-extreme coupling, beyond what either single-condition run reveals alone (Table 8). Graph-level shifts under conditioning. Figure 5 illustrates these structural shifts, summarizing the edge frequencies and mean edge weights under the unconditioned posterior and the three conditional posteriors. Takeaway. Without conditioning, posterior samples often imply zero or near-zero total effect, or distribute the effect across many paths that partially cancel each other, so pathway summaries are hard to determine. Single-effect constraints make the dominant routes clearer: Cond-PIP explains PIP3→PIP2PIP3 2 mainly through two routes, whereas Cond-ERK keeps Erk→AktErk mostly direct but yields a broader set of PIP3→PIP2PIP3 2 paths with partial cancellation (Table 8). Joint conditioning is not just the overlap of the two single-effect results: compared to Cond-ERK, it more often highlights the mediated route Erk→Plcg→AktErk , pointing to PlcgPlcg as a shared mediator that emerges under the joint extreme-effect query (Table 8). 7 Conclusions We introduced conditional causal discovery, a framework for posterior inference over causal graph structures and edge weights under user-specified causal-effect constraints. The framework is designed to answer two questions simultaneously: how likely a specified extreme-effect event is under the posterior, and which graph–weight configurations remain plausible when that event occurs. To make such inference practical in rare-event regions, we developed an adaptive multilevel splitting framework with an MCMC kernel over the joint graph–weight space. By combining structure proposals, such as PARNI-DAG or Structure-MCMC, with blocked weight moves, the method estimates one-sided signed causal-effect tail probabilities and produces representative samples from the corresponding constrained posterior. Empirically, we validated the estimator against exhaustive enumeration in a four-node setting and showed that the approach remains effective as problem size grows, where unconditional baselines often degenerate in the rare-event region. In a case study with the Sachs dataset, conditional samples provided compact pathway-level explanations and highlighted coherent mechanisms under both single-effect and joint-effect queries. Limitations and future work. Our current implementation assumes linear-Gaussian SEMs under causal sufficiency. Extending the framework to nonlinear mechanisms, latent confounding, and interventional data remains an important direction. At the algorithmic level, the multilevel splitting procedure only requires a scalar score function defining the event of interest. For nonlinear structural causal models, the main changes would be to replace the BGe and conjugate posterior components with suitable nonlinear posterior inference modules, and to replace the closed-form linear causal-effect evaluator with an appropriate effect-estimation procedure. More broadly, incorporating richer constraint families, such as path-specific effects or qualitative monotonicity constraints, could enable deeper mechanistic analysis of complex scientific systems. References Anand et al. (2023) T. V. Anand, A. H. Ribeiro, J. Tian, and E. Bareinboim Causal effect identification in cluster dags. In Proceedings of the AAAI Conference on Artificial Intelligence, Cited by: §2. Annadani et al. (2021) Y. Annadani, J. Rothfuss, A. Lacoste, N. Scherrer, A. Goyal, Y. Bengio, and S. Bauer Variational causal networks: approximate bayesian inference over causal structures. arXiv preprint arXiv:2106.07635. Cited by: §2. Bodik et al. (2024) J. Bodik, M. Paluš, and Z. Pawlas Causality in extremes of time series. Extremes 27 (1), p. 67–121. Cited by: §2. Borboudakis and Tsamardinos (2012) G. Borboudakis and I. Tsamardinos Incorporating causal prior knowledge as path-constraints in bayesian networks and maximal ancestral graphs. In Proceedings of the 29th International Coference on International Conference on Machine Learning, p. 427–434. Cited by: §2. Brouillard et al. (2020) P. Brouillard, S. Lachapelle, A. Lacoste, S. Lacoste-Julien, and A. Drouin Differentiable causal discovery from interventional data. Advances in Neural Information Processing Systems 33, p. 21865–21877. Cited by: §2. Cérou and Guyader (2007) F. Cérou and A. Guyader Adaptive multilevel splitting for rare event analysis. Stochastic Analysis and Applications 25 (2), p. 417–443. Cited by: §1, §2, §4.2. Chen et al. (2016) E. Y. Chen, Y. Shen, A. Choi, and A. Darwiche Learning bayesian networks with ancestral constraints. Advances in Neural Information Processing Systems 29. Cited by: §2. Chen and Darwiche (2024) Y. Chen and A. Darwiche Identifying causal effects under functional dependencies. Advances in Neural Information Processing Systems 37, p. 87013–87041. Cited by: §2. Cundy et al. (2021) C. Cundy, A. Grover, and S. Ermon Bcd nets: scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems 34, p. 7095–7110. Cited by: §2. Deleu et al. (2022) T. Deleu, A. Góis, C. Emezue, M. Rankawat, S. Lacoste-Julien, S. Bauer, and Y. Bengio Bayesian structure learning with generative flow networks. In Uncertainty in Artificial Intelligence, p. 518–528. Cited by: §2. Emezue et al. (2023) C. C. Emezue, A. Drouin, T. Deleu, S. Bauer, and Y. Bengio Benchmarking bayesian causal discovery methods for downstream treatment effect estimation. arXiv preprint arXiv:2307.04988. Cited by: §2. Friedman and Koller (2003) N. Friedman and D. Koller Being bayesian about network structure. a bayesian approach to structure discovery in bayesian networks. Machine learning 50, p. 95–125. Cited by: §2, §6. Geiger and Heckerman (1994) D. Geiger and D. Heckerman Learning gaussian networks. In Uncertainty in Artificial Intelligence, p. 235–243. Cited by: §3. Geiger and Heckerman (2002) D. Geiger and D. Heckerman Parameter priors for directed acyclic graphical models and the characterization of several probability distributions. The Annals of Statistics 30 (5), p. 1412–1440. Cited by: §3. Giudice et al. (2023) E. Giudice, J. Kuipers, and G. Moffa A bayesian take on gaussian process networks. Advances in Neural Information Processing Systems 36, p. 56602–56614. Cited by: §2. Glymour et al. (2019) C. Glymour, K. Zhang, and P. Spirtes Review of causal discovery methods based on graphical models. Frontiers in genetics 10, p. 524. Cited by: §3. Gnecco et al. (2021) N. Gnecco, N. Meinshausen, J. Peters, and S. Engelke Causal discovery in heavy-tailed models. The Annals of Statistics 49 (3), p. 1755–1778. Cited by: §2. Guyader et al. (2011) A. Guyader, N. Hengartner, and E. Matzner-Løber Simulation and estimation of extreme quantiles and extreme probabilities. Applied Mathematics & Optimization 64 (2), p. 171–196. Cited by: §4.2. Hauser and Bühlmann (2015) A. Hauser and P. Bühlmann Jointly interventional and observational data: estimation of interventional markov equivalence classes of directed acyclic graphs. Journal of the Royal Statistical Society Series B: Statistical Methodology 77 (1), p. 291–318. Cited by: §2. Kahn and Harris (1951) H. Kahn and T. E. Harris Estimation of particle transmission by random sampling. National Bureau of Standards applied mathematics series 12, p. 27–30. Cited by: §4.2. Koch et al. (2009) M. Koch, B. M. Broom, and D. Subramanian Learning robust cell signalling models from high throughput proteomic data. International journal of bioinformatics research and applications 5 (3), p. 241–253. Cited by: §6, §6. Koller and Friedman (2009) D. Koller and N. Friedman Probabilistic graphical models: principles and techniques. MIT press. Cited by: §3. Kuipers and Moffa (2017) J. Kuipers and G. Moffa Partition mcmc for inference on acyclic digraphs. Journal of the American Statistical Association 112 (517), p. 282–299. Cited by: §2. Liang et al. (2023) X. Liang, A. Caron, S. Livingstone, and J. Griffin Structure learning with adaptive random neighborhood informed mcmc. Advances in Neural Information Processing Systems 36, p. 40760–40772. Cited by: item 4, item 6, §2, §4.4, §5.1. Lorch et al. (2021) L. Lorch, J. Rothfuss, B. Schölkopf, and A. Krause Dibs: differentiable bayesian structure learning. Advances in Neural Information Processing Systems 34, p. 24111–24123. Cited by: item 3, §2, §5.1. Madhanagopal and Amrhein (2019) B. Madhanagopal and J. Amrhein Analyzing structural causal models using the calis procedure. In SAS Global Forum, p. 3240–2019. Cited by: §6. Meek (1995) C. Meek Causal inference and causal explanation with background knowledge. In Proceedings of the Eleventh conference on Uncertainty in artificial intelligence, p. 403–410. Cited by: §2. Pasche et al. (2023) O. C. Pasche, V. Chavez-Demoulin, and A. C. Davison Causal modelling of heavy-tailed variables and confounders with application to river flow. Extremes 26 (3), p. 573–594. Cited by: §2. Pearl (2009) J. Pearl Causality: models, reasoning and inference. 2nd edition, Cambridge University Press, USA. External Links: ISBN 052189560X Cited by: §3. Rittel and Tschiatschek (2023) S. Rittel and S. Tschiatschek Specifying prior beliefs over dags in deep bayesian causal structure learning. In ECAI 2023: 26th European Conference on Artificial Intelligence, including 12th Conference on Prestigious Applications of Intelligent Systems, PAIS 2023-Proceedings, p. 1962–1969. Cited by: §2. Sachs et al. (2005) K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan Causal protein-signaling networks derived from multiparameter single-cell data. Science 308 (5721), p. 523–529. Cited by: §1, §6, §6, §6. Sobel (1990) M. E. Sobel Effect analysis and causation in linear structural equation models. Psychometrika 55 (3), p. 495–515. Cited by: §3. Spirtes et al. (2000) P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman Causation, prediction, and search. MIT press. Cited by: §3. Toth et al. (2024) C. Toth, C. Knoll, F. Pernkopf, and R. Peharz Effective bayesian causal inference via structural marginalisation and autoregressive orders. arXiv preprint arXiv:2402.14781. Cited by: §2. Toth et al. (2022) C. Toth, L. Lorch, C. Knoll, A. Krause, F. Pernkopf, R. Peharz, and J. Von Kügelgen Active bayesian causal inference. Advances in Neural Information Processing Systems 35, p. 16261–16275. Cited by: §2. Viinikka et al. (2020) J. Viinikka, A. Hyttinen, J. Pensar, and M. Koivisto Towards scalable bayesian learning of causal dags. Advances in Neural Information Processing Systems 33, p. 6584–6594. Cited by: §2, §3. Wang et al. (2022) B. Wang, M. R. Wicker, and M. Kwiatkowska Tractable uncertainty for structure learning. In International Conference on Machine Learning, p. 23131–23150. Cited by: item 2, §2, §5.1. Appendix A Appendix A.1 Experiment Setting For all synthetic experiments, we generate the ground-truth DAG G⋆G using an ordered Erdős–Rényi construction. Given dimension d and edge budget parameter edges_per_node=2 edges\_per\_node=2, we set pedge=minedges_per_node⋅d(d−1)/2, 0.5.p_edge= \! \ edges\_per\_node· dd(d-1)/2,\,0.5 \. We then sample a Bernoulli adjacency matrix with edge probability pedgep_edge, keep only the strict lower-triangular part to ensure acyclicity under the current ordering, and finally apply a random permutation of the node labels. The effective edge probabilities are therefore pedge=0.5p_edge=0.5 for d=4d=4, pedge=0.5p_edge=0.5 for d=8d=8, pedge=4/15p_edge=4/15 for d=16d=16, and pedge=4/31p_edge=4/31 for d=32d=32. The cap at 0.50.5 is the small-graph handling used by the generator: for d=4d=4, the nominal value 4/34/3 is clipped to 0.50.5, giving an expected 33 edges instead of an overly dense four-node DAG. The corresponding expected edge counts are 33, 1414, 3232, and 6464 for d=4,8,16,32d=4,8,16,32, respectively. Conditional on G⋆G , absent edges are assigned weight zero and nonzero edge weights are sampled independently from a Gaussian distribution, as in the data generation script. Observations are then generated from the linear-Gaussian SEM using nobs=1000n_obs=1000 training samples. Below are the experiment settings for each method: 1. Exhaustive enumeration (d=4d=4 only): enumerate all DAGs (|d=4|=543|G_d=4|=543) and draw 2000020000 weight samples per graph. 2. OrderSPN [37]: sample 50005000 graphs from the learned circuit and draw 100100 weight samples per graph. 3. DiBS [25]: draw 500000500000 joint graph–weight samples in total. 4. PARNI-DAG single chain [24]: run a single long MCMC chain over (G,B)(G,B) using PARNI-DAG as the structure proposal kernel for 500000500000 iterations with 10%10\% burn-in. 5. Structure-MCMC single chain: run a single long MCMC chain over (G,B)(G,B) using Structure-MCMC as the structure proposal kernel for 500000500000 iterations with 10%10\% burn-in. 6. MLS-PARNI-DAG [24]: use our MLS framework with PARNI-DAG as the structure proposal kernel. We use N=200N=200 particles, m=2000m=2000 MCMC iterations per level, and Kmax=10K_ =10 levels for all dimensions. 7. MLS-Structure-MCMC: use our MLS framework with Structure-MCMC as the structure proposal kernel, with the same MLS hyperparameters as above. A.2 Target-threshold construction We construct target score levels to evaluate each method’s ability to estimate deep-tail probabilities under increasingly extreme one-sided causal-effect constraints. Let e0=CEij(Z⋆)=CEij(G⋆,B⋆)e_0=CE_ij(Z )=CE_ij(G ,B ) denote the ground-truth signed total causal effect for a queried ordered pair (i,j)(i,j), where Z⋆=(G⋆,B⋆)Z =(G ,B ) is the data-generating graph–weight state. For synthetic experiments, we choose the tail direction using the sign of e0e_0. Equivalently, define s=sign(e0)s=sign(e_0) and use the scalar score h(Z)=sCEij(Z).h(Z)=s\,CE_ij(Z). Then larger values of h(Z)h(Z) always correspond to more extreme effects in the selected direction. A target score level λ defines the event λ=Z:h(Z)≥λ.A_λ=\Z:h(Z)≥λ\. When e0>0e_0>0, this is the right-tail event ℰij+(λ)=Z:CEij(Z)≥λE_ij^+(λ)=\Z:CE_ij(Z)≥λ\. When e0<0e_0<0, it is the left-tail event ℰij−(λ)=Z:CEij(Z)≤−λE_ij^-(λ)=\Z:CE_ij(Z)≤-λ\. This is the same score-level representation used by the MLS sampler in Sec. 4. 1. Initial score level. We initialize the threshold grid at the magnitude of the ground-truth effect, λ1=|e0| _1=|e_0|. 2. Pilot run. We run a short, low-budget adaptive multilevel-splitting pilot to identify a more extreme target score level λT _T such that the corresponding tail probability is already in the rare-event regime, approximately 10−610^-6 or smaller. 3. Threshold grid. We form a monotone increasing sequence of T target score levels λ1<λ2<⋯<λT _1< _2<·s< _T. Larger λ always corresponds to a more extreme one-sided effect in the selected direction. A.3 Detailed result for d=4 baseline experiment Table 3 lists statistics for the d=4d=4 baseline experiments across seven methods. Table 3: d=4d=4 per-case probability estimates across runs (mean ± SD). MLS refers to multilevel splitting, and SC refers to single-chain sampling. Case Enumeration p^±SD p PARNI-MLS p^±SD p Structure-MCMC-MLS p^±SD p SPN p^±SD p PARNI-SC p^±SD p Structure-MCMC-SC p^±SD p DIBS p^±SD p 1 1.87×10−51.87× 10^-5 ± 7.25×10−6±\,7.25× 10^-6 1.83×10−51.83× 10^-5 ± 1.07×10−5±\,1.07× 10^-5 1.9×10−51.9× 10^-5 ± 6.09×10−6±\,6.09× 10^-6 1.42×10−51.42× 10^-5 ± 6×10−6±\,6× 10^-6 1.84×10−51.84× 10^-5 ± 5.32×10−6±\,5.32× 10^-6 1.4×10−51.4× 10^-5 ± 1.01×10−5±\,1.01× 10^-5 2.32×10−22.32× 10^-2 ± 1.35×10−3±\,1.35× 10^-3 2 2.3×10−42.3× 10^-4 ± 1.13×10−4±\,1.13× 10^-4 2.5×10−42.5× 10^-4 ± 1.07×10−4±\,1.07× 10^-4 2.41×10−42.41× 10^-4 ± 8.04×10−5±\,8.04× 10^-5 1.37×10−21.37× 10^-2 ± 8.62×10−4±\,8.62× 10^-4 2.28×10−42.28× 10^-4 ± 7.49×10−5±\,7.49× 10^-5 2.53×10−42.53× 10^-4 ± 5.78×10−5±\,5.78× 10^-5 1.35×10−21.35× 10^-2 ± 7.72×10−4±\,7.72× 10^-4 3 1.65×10−41.65× 10^-4 ± 5.96×10−5±\,5.96× 10^-5 1.7×10−41.7× 10^-4 ± 5.47×10−5±\,5.47× 10^-5 1.61×10−41.61× 10^-4 ± 1.41×10−5±\,1.41× 10^-5 1.83×10−11.83× 10^-1 ± 7.92×10−3±\,7.92× 10^-3 1.67×10−41.67× 10^-4 ± 2.42×10−5±\,2.42× 10^-5 1.79×10−41.79× 10^-4 ± 2.31×10−5±\,2.31× 10^-5 1.81×10−11.81× 10^-1 ± 8.38×10−3±\,8.38× 10^-3 4 6.17×10−56.17× 10^-5 ± 4.07×10−5±\,4.07× 10^-5 5.89×10−55.89× 10^-5 ± 3.67×10−5±\,3.67× 10^-5 6.79×10−56.79× 10^-5 ± 2.13×10−5±\,2.13× 10^-5 2.67×10−22.67× 10^-2 ± 1.76×10−3±\,1.76× 10^-3 7.06×10−57.06× 10^-5 ± 2.85×10−5±\,2.85× 10^-5 7.21×10−57.21× 10^-5 ± 3.06×10−5±\,3.06× 10^-5 2.65×10−22.65× 10^-2 ± 1.66×10−3±\,1.66× 10^-3 5 1.6×10−41.6× 10^-4 ± 5.56×10−5±\,5.56× 10^-5 1.71×10−41.71× 10^-4 ± 5.9×10−5±\,5.9× 10^-5 1.62×10−41.62× 10^-4 ± 1.74×10−5±\,1.74× 10^-5 3.34×10−23.34× 10^-2 ± 1.88×10−3±\,1.88× 10^-3 1.55×10−41.55× 10^-4 ± 4.54×10−5±\,4.54× 10^-5 1.69×10−41.69× 10^-4 ± 4.88×10−5±\,4.88× 10^-5 3.34×10−23.34× 10^-2 ± 1.63×10−3±\,1.63× 10^-3 6 2.15×10−32.15× 10^-3 ± 7.08×10−4±\,7.08× 10^-4 2.21×10−32.21× 10^-3 ± 9.14×10−4±\,9.14× 10^-4 2.16×10−32.16× 10^-3 ± 7.4×10−4±\,7.4× 10^-4 2.27×10−22.27× 10^-2 ± 1.61×10−3±\,1.61× 10^-3 2.09×10−32.09× 10^-3 ± 1.38×10−4±\,1.38× 10^-4 2.15×10−32.15× 10^-3 ± 1.78×10−4±\,1.78× 10^-4 2.3×10−22.3× 10^-2 ± 1.55×10−3±\,1.55× 10^-3 7 2.44×10−42.44× 10^-4 ± 1.13×10−4±\,1.13× 10^-4 2.4×10−42.4× 10^-4 ± 1.06×10−4±\,1.06× 10^-4 2.19×10−42.19× 10^-4 ± 2.48×10−5±\,2.48× 10^-5 9.71×10−29.71× 10^-2 ± 6.54×10−3±\,6.54× 10^-3 2.22×10−42.22× 10^-4 ± 2.25×10−5±\,2.25× 10^-5 2.42×10−42.42× 10^-4 ± 2.5×10−5±\,2.5× 10^-5 9.65×10−29.65× 10^-2 ± 6.06×10−3±\,6.06× 10^-3 8 3.09×10−53.09× 10^-5 ± 2.84×10−5±\,2.84× 10^-5 4.48×10−54.48× 10^-5 ± 2.71×10−5±\,2.71× 10^-5 3.1×10−53.1× 10^-5 ± 6.72×10−6±\,6.72× 10^-6 5.76×10−35.76× 10^-3 ± 6.65×10−4±\,6.65× 10^-4 2.82×10−52.82× 10^-5 ± 8.1×10−6±\,8.1× 10^-6 3.4×10−53.4× 10^-5 ± 1.35×10−5±\,1.35× 10^-5 5.86×10−35.86× 10^-3 ± 5.4×10−4±\,5.4× 10^-4 9 3.22×10−43.22× 10^-4 ± 1.34×10−4±\,1.34× 10^-4 3.25×10−43.25× 10^-4 ± 1.33×10−4±\,1.33× 10^-4 3.22×10−43.22× 10^-4 ± 5.19×10−5±\,5.19× 10^-5 1.32×10−21.32× 10^-2 ± 1.46×10−3±\,1.46× 10^-3 3.09×10−43.09× 10^-4 ± 6.17×10−5±\,6.17× 10^-5 3.3×10−43.3× 10^-4 ± 4.04×10−5±\,4.04× 10^-5 1.33×10−21.33× 10^-2 ± 1.45×10−3±\,1.45× 10^-3 10 1.33×10−51.33× 10^-5 ± 1.02×10−5±\,1.02× 10^-5 1.98×10−51.98× 10^-5 ± 1.78×10−5±\,1.78× 10^-5 1.37×10−51.37× 10^-5 ± 1.11×10−5±\,1.11× 10^-5 1.82×10−21.82× 10^-2 ± 1.55×10−3±\,1.55× 10^-3 1.68×10−51.68× 10^-5 ± 1.52×10−5±\,1.52× 10^-5 1.21×10−51.21× 10^-5 ± 1.58×10−5±\,1.58× 10^-5 1.8×10−21.8× 10^-2 ± 1.72×10−3±\,1.72× 10^-3 A.4 Time Complexity We separate the cost of the outer adaptive multilevel-splitting loop from the cost of the inner MCMC mutation kernel. Let N be the number of particles, m the number of MCMC mutation steps per particle per level, K the number of splitting levels, and CMHC_MH the average cost of one Metropolis– Hastings proposal, including proposal generation, posterior-ratio evaluation, and causal-effect score evaluation. At each level, the algorithm mutates N particles for m steps, computes N scores, sorts or partially sorts the scores to choose the next empirical quantile, and resamples the survivors. The overall cost is therefore O(KNmCMH+KNlogN),O\! (KNmC_MH+KN N ), where the KNlogNKN N term comes from quantile selection. If selection is implemented by a linear-time order-statistic routine, this sorting term can be reduced to O(KN)O(KN) and the mutation cost dominates. The average proposal cost depends on the mixture of structure and weight moves. If pstructp_struct is the probability of proposing a structure move, then CMH≈pstructCstruct+(1−pstruct)Cweight+Ch,C_MH≈ p_structC_struct+(1-p_struct)C_weight+C_h, where CstructC_struct is the cost of proposing and scoring a graph update, CweightC_weight is the cost of a blocked coefficient refresh, and ChC_h is the cost of evaluating the scalar score h(Z)h(Z). In the linear-Gaussian implementation, the coefficient refreshes use node-wise conjugate posterior updates, while h(Z)h(Z) is computed from the linear total-effect matrix. A direct matrix inverse gives Ch=O(d3)C_h=O(d^3), although for a DAG and a single target pair one can exploit the topological structure or solve a triangular system to reduce this cost in practice. The key rare-event advantage is the dependence on the event probability p⋆p_ . A direct posterior sampler needs about 1/p⋆1/p_ samples to see one sample from an event of posterior probability p⋆p_ , and O(1/(r2p⋆))O(1/(r^2p_ )) samples to estimate that probability with fixed relative error r. By contrast, adaptive multilevel splitting replaces the single rare event with a product of moderate conditional survival probabilities. If the empirical survival fraction is approximately ρ at each level, then K grows roughly like log(p⋆)/log(ρ) (p_ )/ (ρ), so the leading cost grows approximately logarithmically in 1/p⋆1/p_ rather than linearly in 1/p⋆1/p_ . This is why multilevel splitting remains useful in the deep-tail regimes where unconditional sampling degenerates. Runtime measurements. Table 4 reports wall-clock running times for the d=4d=4 case 1 experiment. Table 5 reports the average time per outer MLS loop for the two structure kernels across dimensions. PARNI-DAG is more expensive per loop because its locally informed proposal requires additional neighborhood construction and scoring, while Structure-MCMC has cheaper local edge proposals. The benefit of PARNI-DAG is not lower per-iteration cost, but more stable deep-tail behavior in larger graphs, as shown in the scalability experiments. Table 4: Average running time for the d=4d=4 case 1 experiment over 10 runs. We report mean ± standard deviation in seconds. Method Time (s) MLS-PARNI-DAG 192.3±4.8192.3± 4.8 MLS-Structure-MCMC 56.1±2.156.1± 2.1 PARNI-DAG 268.2±5.3268.2± 5.3 Structure-MCMC 105.5±3.0105.5± 3.0 Enumeration-Baseline 301.6±4.6301.6± 4.6 OrderSPN 166.8±3.9166.8± 3.9 DiBS 434.2±6.7434.2± 6.7 Table 5: Average running time per outer MLS loop in seconds. Standard deviations over 10 runs are shown in parentheses. Kernel d=4d=4 d=8d=8 d=16d=16 d=32d=32 PARNI-DAG 38.4 (1.1) 56.8 (1.6) 92.4 (2.8) 172.5 (4.9) Structure-MCMC 11.2 (0.5) 19.5 (0.8) 34.8 (1.2) 58.5 (2.1) A.5 Extension to nonlinear Gaussian mechanisms The outer MLS framework is not specific to linear-Gaussian SEMs. It only requires three ingredients: an unconstrained posterior target π(Z)π(Z), a scalar score h(Z)h(Z) whose super-level set defines the event of interest, and an MCMC kernel that approximately preserves the level-truncated posterior π(Z)h(Z)≥λπ(Z)1\h(Z)≥λ\. The linear-Gaussian assumptions used in the main experiments provide convenient closed forms for the BGe marginal likelihood, the node-wise coefficient posterior, and the total causal effect, but they are not required by the splitting principle itself. For a nonlinear Gaussian structural causal model, one could write Xj=fj(XpaG(j),ϕj)+ϵj,ϵj∼(0,σj2),X_j=f_j\! (X_pa_G(j); _j )+ _j, _j (0, _j^2), and replace the graph–weight state Z=(G,B)Z=(G,B) by a graph–mechanism state Z=(G,ϕ,σ)Z=(G,φ,σ). The posterior target would become p(G,ϕ,σ∣)p(G,φ,σ ), obtained using an appropriate nonlinear mechanism class such as splines, Gaussian processes, or neural networks. The inner mutation kernel would then update graph structure and mechanism parameters instead of graph structure and linear coefficients. The causal-effect score would also be replaced by a nonlinear effect evaluator. For example, for a scalar intervention one may define an average interventional contrast, h(Z)=Z[Xj∣do(Xi=x+Δ)]−Z[Xj∣do(Xi=x)],h(Z)=E_Z\! [X_j (X_i=x+ ) ]-E_Z\! [X_j (X_i=x) ], or an averaged derivative when such derivatives are well-defined. These expectations can be estimated by ancestral simulation under the proposed nonlinear SCM. The outer splitting loop would remain unchanged; the additional computational burden would come from evaluating nonlinear posterior ratios and interventional scores. Therefore, the main technical requirement for such an extension is not a new rare-event algorithm, but reliable posterior inference and causal-effect evaluation for the chosen nonlinear mechanism class. A.6 Sachs experiment calculations and tables This section provides supplementary numerical details for the Sachs case study in Sec. 6, including (i) the qualitative reference network, (i) estimated posterior masses of the tail events, and (i) pathway-level decompositions of conditional effects. Reference network. Figure 4 shows the 11-node, 20-edge consensus signaling network used as a qualitative reference in Sec. 6. Figure 4: Consensus signaling network (Sachs) used as a qualitative reference. Posterior mass of tail events. We run MLS with 10 independent repetitions for each conditioning regime. Table 6 reports (i) the geometric mean estimate of Pr(ℰ∣) (E ) (i.e., exp([logp^]) (E[ p]) across runs) and (i) the mean ± sd of −logp^- p across runs. Table 6: Estimated posterior mass of Sachs tail events over 10 MLS runs. Condition Pr(ℰ∣) (E ) (geo-mean) −logp^- p (mean ± sd) Cond-PIP 2.2×10−72.2× 10^-7 15.32±0.4415.32± 0.44 Cond-ERK 1.3×10−71.3× 10^-7 15.84±0.5315.84± 0.53 Cond-Joint 5.2×10−75.2× 10^-7 14.47±0.6714.47± 0.67 Path decomposition under Cond-PIP. Across samples drawn from the conditional posterior p(Z∣,ℰPIP(tPIP)),p\! (Z ,E_PIP(t_PIP) ), the two canonical pathways PIP3→PIP2PIP3\!→\!PIP2 and PIP3→Plcg→PIP2PIP3\!→\!Plcg\!→\!PIP2 are present in 100% of samples. These two routes account for essentially the entire conditional mean effect: [CEPIP3,PIP2(Z)∣,ℰPIP(tPIP)] \! [CE_PIP3,PIP2(Z) ,E_PIP(t_PIP) ] (14) ≈0.6142⏟direct+0.1614⏟via Plcg+2⋅10−4⏟other paths≈0.7758. ≈ 0.6142_direct+ 0.1614_via Plcg+ 2· 10^-4_other paths≈ 0.7758. The contribution from all remaining directed paths is numerically negligible, confirming that conditioning isolates a highly concentrated two-path mechanism. Summary tables. Table 7 reports posterior effect statistics (e.g., Pr(CE>0) (CE>0), conditional means, and path counts) across the four regimes. Table 8 further decomposes the conditional mean effects into direct and mediated contributions for the two target pairs (PIP3→ 2 and Erk→ ), together with the frequency of the canonical mediated paths. (a) Unconditioned posterior (b) Cond-PIP: CEPIP3→PIP2>tPIPCE_PIP3 2>t_PIP (c) Cond-ERK: CEErk→Akt>tERKCE_Erk >t_ERK (d) Cond-Joint: both constraints Figure 5: Heatmaps summarizing posterior graph features under four conditioning settings. Each panel contains two heatmaps: edge frequency (probability an edge appears in sampled graphs; left) and mean edge weight (conditional mean of BuvB_uv given the edge is present; right). Table 7: Comparison across four conditioning settings for two target pairs. Pr(CE>0) (CE>0) is the posterior mass with a directed influence (i.e., at least one directed path), [CE∣CE>0]E[CE >0] and quantiles are computed on the nonzero subset, and #paths is the median number of directed paths among samples with CE>0CE>0. We run PARNI-DAG single chain with 200,000200,000 iterations and collect 180000180000 samples for the unconditioned case. For each of the other conditions, we collect 500500 survivor samples from 10 runs with 5050 samples for each run. PIP3→ 2 Erk→ Condition Pr(CE>0) (CE>0) [CE∣CE>0]E[CE >0] #pathsmed_med Pr(CE>0) (CE>0) [CE∣CE>0]E[CE >0] #pathsmed_med Unconditioned 0.109 0.600 [0.543, 0.681] 4 0.085 0.569 [0.532, 0.617] 1 Cond-PIP (CEPIP>tPIPCE_PIP>t_PIP) 1.000 0.776 [0.771, 0.782] 2 0.498 0.596 [0.538, 0.634] 1 Cond-ERK (CEERK>tERKCE_ERK>t_ERK) 0.984 0.655 [0.562, 0.710] 10 1.000 0.673 [0.670, 0.676] 1 Cond-Joint (CEPIP>tPIP′∧CEERK>tERK′CE_PIP>t _PIP _ERK>t _ERK) 1.000 0.747 [0.742, 0.758] 3 1.000 0.656 [0.652, 0.662] 1 Table 8: Mechanism decomposition across conditions, using the two pathways most relevant for interpretation: for PIP3→ 2 we use the direct path (PIP3→ 2) and the Plcg-mediated path (PIP3→ → 2); for Erk→ we use the direct path (Erk→ ) and the Plcg-mediated path (Erk→ → ). For unconditioned samples, we select those with CE>0CE>0, and calculate [CE∣CE>0]E[CE >0]. The last column reports the frequency that the canonical mediated path edges are present (over all samples under the condition). Pair Condition [CE∣CE>0]E[CE >0] Direct contrib Via-Plcg contrib Pr(Via-Plcg path) (Via-Plcg path) PIP3→ 2 Unconditioned 0.600 0.560 (93.4%) 0.076 (12.7%) 0.064 Cond-PIP 0.776 0.614 (79.2%) 0.161 (20.8%) 1.000 Cond-ERK 0.655 0.559 (85.2%) 0.140 (21.4%) 0.975 Cond-Joint 0.747 0.598 (79.9%) 0.154 (20.5%) 1.000 Erk→ Unconditioned 0.569 0.563 (98.9%) 0.000 (0.0%) 0.0004 Cond-PIP 0.596 0.606 (101.6%) 0.001 (0.2%) 0.026 Cond-ERK 0.673 0.668 (99.3%) 0.004 (0.6%) 0.177 Cond-Joint 0.656 0.650 (99.1%) 0.005 (0.8%) 0.244 Note: Percentages are (path contribution)/(mean effect). They may exceed 100% or be negative due to path cancellations. A.7 Synthetic multi-effect heat maps Figures 6 and 7 provide the detailed heat-map visualizations for the multi-effect conditioning experiments discussed in Sec. 5.2.2. Each row displays the aggregated edge frequency and average edge weight summaries for one posterior condition. These figures are included in the appendix because they are visually dense; the main text reports the more compact SHD summaries. (a) Ground-truth graph. (b) Conditional posterior under strong multi-effect constraints. (c) Conditional posterior under weak multi-effect constraints. Figure 6: d=4d=4 multi-effect conditioning results. Each panel summarizes edge frequencies and average edge weights. Stronger constraints concentrate posterior mass more tightly around the ground-truth mechanism. (a) Ground-truth graph. (b) Conditional posterior under strong multi-effect constraints. (c) Conditional posterior under weak multi-effect constraints. Figure 7: d=8d=8 multi-effect conditioning results, with the same layout as Fig. 6. The strong constraints produce sharper edge and weight summaries than the weak constraints. A.8 Pseudo-code Algorithm 2 summarizes the adaptive multilevel-splitting loop, and Algorithm 3 summarizes the level-truncated MCMC mutation kernel. The notation matches Sec. 4: h is the scalar score, λ⋆ _ is the target level, and λ=Z:h(Z)≥λA_λ=\Z:h(Z)≥λ\ is the corresponding score-level event. Algorithm 2 Adaptive Multilevel Splitting for Conditional Causal Discovery 1: Data D; posterior density π(Z)=p(G,B∣)π(Z)=p(G,B ); score h; target level λ⋆ _ ; particle size N; survival fraction ρ; mutation steps m per level; maximum number of levels KmaxK_ ; structure kernel Kernel. 2: Tail-probability estimate p^⋆ p_ and approximately constrained particles P. 3: Initialize particles =Zn=(Gn,Bn)n=1NP=\Z_n=(G_n,B_n)\_n=1^N by drawing Gn∼p(G)G_n p(G) and Bn∼p(B∣Gn,)B_n p(B G_n,D). 4: Set λ0←−∞ _0←-∞, p^⋆←1 p_ ← 1, and k←0k← 0. 5: ←MCMCMutation(,,h,λ0,m,)P← MCMCMutation(P,D,h, _0,m, Kernel). 6: while k<Kmaxk<K_ do 7: Compute scores sn←h(Zn)s_n← h(Z_n) for all Zn∈Z_n . 8: λ~k+1←Q1−ρ(snn=1N) λ_k+1← Q_1-ρ(\s_n\_n=1^N). 9: λk+1←minλ~k+1,λ⋆ _k+1← \ λ_k+1, _ \. 10: ←Zn∈:h(Zn)≥λk+1S←\Z_n :h(Z_n)≥ _k+1\. 11: β^k+1←||/N β_k+1←|S|/N. 12: p^⋆←p^⋆β^k+1 p_ ← p_ \, β_k+1. 13: if ||=0|S|=0 then 14: return (0,∅)(0, ). 15: end if 16: Resample N particles from S with replacement to form P. 17: ←MCMCMutation(,,h,λk+1,m,)P← MCMCMutation(P,D,h, _k+1,m, Kernel). 18: if λk+1=λ⋆ _k+1= _ then 19: return (p^⋆,)( p_ ,P). 20: end if 21: k←k+1k← k+1. 22: end while 23: return (p^⋆,)( p_ ,P). Algorithm 3 Level-Truncated MCMC Mutation over Graph–Weight States 1: Population P; data D; score h; level λ; mutation steps m; structure kernel Kernel; structure-move probability pstructp_struct. 2: Mutated population ′P . 3: ′←∅P ← . 4: for all Z=(G,B)∈Z=(G,B) do 5: for r=1r=1 to m do 6: Draw u∼Uniform(0,1)u (0,1). 7: if u<pstructu<p_struct then 8: Propose a graph G′∼qG(⋅∣G)G q_G(· G) using Kernel. 9: Refresh affected coefficient blocks from p(B∣G′,)p(B G ,D) and copy unchanged blocks to obtain B′B . 10: else 11: Set G′←G ← G and refresh one coefficient block from p(B∣G,)p(B G,D) to obtain B′B . 12: end if 13: Set Z′←(G′,B′)Z ←(G ,B ) and let q(Z′∣Z)q(Z Z) denote the complete proposal density. 14: if h(Z′)<λh(Z )<λ then 15: Reject Z′Z and continue. 16: else 17: α←min1,π(Z′)q(Z∣Z′)π(Z)q(Z′∣Z)α← \! \1, π(Z )q(Z Z )π(Z)q(Z Z) \. 18: Draw a∼Uniform(0,1)a (0,1). 19: if a≤αa≤α then 20: Z←Z′Z← Z . 21: end if 22: end if 23: end for 24: ′←′∪ZP ∪\Z\. 25: end for 26: return ′P .