Interpretable Causal Discovery via Causal-Effect Constraints

# Interpretable Causal Discovery via Causal-Effect Constraints
Source: [https://arxiv.org/html/2608.12640](https://arxiv.org/html/2608.12640)
Guy Van den BroeckAffiliation:Computer Science Dept\. University of California, Los Angeles Los Angeles, California, USABenjie WangAffiliation: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](https://github.com/ZCX031116/MLS-Framework)\.## 1Introduction

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](https://arxiv.org/html/2608.12640#bib.bib12)\], 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\\mid\\mathcal\{D\}\): many plausible DAGsiijjiijjCausal\-effect constraintℰi​j\+​\(t\)=\{\(G,B\):CEi​j​\(G,B\)≥t\}\\mathcal\{E\}\_\{ij\}^\{\+\}\(t\)=\\\{\(G,B\):\\mathrm\{CE\}\_\{ij\}\(G,B\)\\geq t\\\}ttrare posterior eventConditional posteriorp\(G,B∣𝒟,ℰi​j\(t\)\)p\(G,B\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{ij\}\(t\)\)iijjhigh\-frequency pathwaysexplain the eventscore byCEi​j\\mathrm\{CE\}\_\{ij\}MLS\+ 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 nodeiito nodejj\. 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 \(GG\) and parameters \(θ\\theta\)\. Our target is the posterior conditioned on an event determined by the graph and parameters, such as the event that the signed causal effect fromiitojjexceeds a thresholdtt\. 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 thresholdttis 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](https://arxiv.org/html/2608.12640#bib.bib13)\]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\\in\{4,8,16,32\}\) and present a Sachs dataset case study showing how extreme\-effect conditioning enables pathway\-level interpretation and hypothesis generation\.

## 2Related 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](https://arxiv.org/html/2608.12640#bib.bib27),[4](https://arxiv.org/html/2608.12640#bib.bib29),[1](https://arxiv.org/html/2608.12640#bib.bib28)\]\. Prior work has considered both identifiability and algorithmic procedures for learning graphs subject to these constraints\[[7](https://arxiv.org/html/2608.12640#bib.bib31),[8](https://arxiv.org/html/2608.12640#bib.bib30)\]\. A related line of work uses interventional data as an additional source of information beyond pure observations\[[19](https://arxiv.org/html/2608.12640#bib.bib32),[5](https://arxiv.org/html/2608.12640#bib.bib33)\]\. 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](https://arxiv.org/html/2608.12640#bib.bib26),[28](https://arxiv.org/html/2608.12640#bib.bib35),[3](https://arxiv.org/html/2608.12640#bib.bib36)\]\. 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](https://arxiv.org/html/2608.12640#bib.bib14),[23](https://arxiv.org/html/2608.12640#bib.bib15),[36](https://arxiv.org/html/2608.12640#bib.bib6),[15](https://arxiv.org/html/2608.12640#bib.bib16)\]\. More recent work has developed efficient DAG\-space MCMC samplers with locally informed and adaptive proposals, including PARNI\-DAG\[[24](https://arxiv.org/html/2608.12640#bib.bib23)\], 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](https://arxiv.org/html/2608.12640#bib.bib34)\]\.

Alternative approaches approximate the posterior via variational inference over graphs\[[2](https://arxiv.org/html/2608.12640#bib.bib17),[25](https://arxiv.org/html/2608.12640#bib.bib1),[9](https://arxiv.org/html/2608.12640#bib.bib18),[37](https://arxiv.org/html/2608.12640#bib.bib11),[10](https://arxiv.org/html/2608.12640#bib.bib19),[30](https://arxiv.org/html/2608.12640#bib.bib21),[34](https://arxiv.org/html/2608.12640#bib.bib20)\]\. These methods aim either to generate graph samples reflecting structural uncertainty or to support Bayesian model averaging for causal inference\[[35](https://arxiv.org/html/2608.12640#bib.bib22)\]\.

Our computational approach depends on the rare\-event simulation literature, particularly adaptive multilevel splitting\[[6](https://arxiv.org/html/2608.12640#bib.bib13)\]\. 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\.

## 3Preliminaries

##### Causal Bayesian Networks

A Bayesian network \(BN\)\(G,θ\)\(G,\\theta\)is a probabilistic modelp⁡\(𝑿\)p\(\\bm\{X\}\)overddvariables𝑿=X1,…,Xd\\bm\{X\}=\{X\_\{1\},\\ldots,X\_\{d\}\}, specified by a directed acyclic graph \(DAG\)GGand mechanism parametersθ\\theta\. The graph encodes conditional independencies, whileθi\\theta\_\{i\}parameterizes the conditional distribution ofXiX\_\{i\}given its parents\. The joint distribution factorizes as

p⁡\(𝑿∣G,θ\)=∏i=1dp⁡\(Xi∣paG​\(Xi\),θi\),p\(\\bm\{X\}\\mid G,\\theta\)=\\prod\_\{i=1\}^\{d\}p\\left\(X\_\{i\}\\mid\\mathrm\{pa\}\_\{G\}\(X\_\{i\}\),\\theta\_\{i\}\\right\),wherepaG​\(Xi\)\\mathrm\{pa\}\_\{G\}\(X\_\{i\}\)denotes the parents ofXiX\_\{i\}inGG\.

In this paper, the generic mechanism parametersθ\\thetaare instantiated by a linear\-Gaussian structural equation model\. Specifically, the observed variables satisfy𝑿=𝑿​B\+ϵ\\bm\{X\}=\\bm\{X\}B\+\\bm\{\\epsilon\}, whereB∈ℝd×dB\\in\\mathbb\{R\}^\{d\\times d\}is the weighted adjacency matrix andϵ∼𝒩⁡\(𝒃,Σ\)\\bm\{\\epsilon\}\\sim\\mathcal\{N\}\(\\bm\{b\},\\Sigma\), with𝒃∈ℝd\\bm\{b\}\\in\\mathbb\{R\}^\{d\}and diagonalΣ∈ℝ≥0d×d\\Sigma\\in\\mathbb\{R\}\_\{\\geq 0\}^\{d\\times d\}\. For a given DAGGG, we imposeBi​j=0B\_\{ij\}=0wheneveriiis not a parent ofjjinGG\.

Causal Bayesian networks\[[33](https://arxiv.org/html/2608.12640#bib.bib2),[29](https://arxiv.org/html/2608.12640#bib.bib3)\]add a causal interpretation to the directed edges inGG: they describe how the joint distribution changes under interventions\. In the linear\-Gaussian SEM, the signed total causal effect fromXiX\_\{i\}toXjX\_\{j\}is the derivative of the post\-intervention mean ofXjX\_\{j\}with respect to an intervention onXiX\_\{i\}\. It has the closed form

CEi​j​\(G,B\)≜CE⁡\(i→j∣G,B\)=\[\(I−B\)−1\]i​j\\mathrm\{CE\}\_\{ij\}\(G,B\)\\triangleq\\mathrm\{CE\}\(i\\to j\\mid G,B\)=\\left\[\(I\-B\)^\{\-1\}\\right\]\_\{ij\}\(1\)following the standard total\-effect formula for linear structural equation models\[[32](https://arxiv.org/html/2608.12640#bib.bib37)\]\. SinceGGis acyclic,BBis nilpotent after a topological ordering, so\(I−B\)−1=I\+B\+B2\+⋯\+Bd−1\(I\-B\)^\{\-1\}=I\+B\+B^\{2\}\+\\cdots\+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 fromXiX\_\{i\}toXjX\_\{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 asZ=\(G,B\)Z=\(G,B\), and abbreviate the signed causal effect asCEi​j​\(Z\)=CEi​j​\(G,B\)\\mathrm\{CE\}\_\{ij\}\(Z\)=\\mathrm\{CE\}\_\{ij\}\(G,B\)\. For a target ordered pair\(i,j\)\(i,j\)and thresholdt\>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

ℰi​j\+​\(t\)\\displaystyle\\mathcal\{E\}\_\{ij\}^\{\+\}\(t\)=\{Z:CEi​j​\(Z\)≥t\},\\displaystyle=\\left\\\{Z:\\mathrm\{CE\}\_\{ij\}\(Z\)\\geq t\\right\\\},\(2\)ℰi​j−​\(t\)\\displaystyle\\mathcal\{E\}\_\{ij\}^\{\-\}\(t\)=\{Z:CEi​j​\(Z\)≤−t\}\.\\displaystyle=\\left\\\{Z:\\mathrm\{CE\}\_\{ij\}\(Z\)\\leq\-t\\right\\\}\.When the sign is clear from context, we writeℰi​j​\(t\)\\mathcal\{E\}\_\{ij\}\(t\)for either tail event\. These events play two roles in our inference problem\. First, we estimate their posterior probability, such asℙ⁡\(ℰi​j±​\(t\)∣𝒟\)\\mathbb\{P\}\(\\mathcal\{E\}\_\{ij\}^\{\\pm\}\(t\)\\mid\\mathcal\{D\}\)\. Second, we use them to define the constrained posteriorp\(G,B∣𝒟,ℰi​j±\(t\)\)p\\left\(G,B\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{ij\}^\{\\pm\}\(t\)\\right\), from which we draw graph–weight samples for pathway\-level summaries\.

##### Bayesian Causal Discovery

Causal discovery\[[22](https://arxiv.org/html/2608.12640#bib.bib4),[16](https://arxiv.org/html/2608.12640#bib.bib5)\]is the problem of inferring the DAGGGresponsible for generating an observed dataset𝒟\\mathcal\{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 priorp⁡\(G\)p\(G\)over DAGs and use the BGe marginal likelihoodp⁡\(𝒟∣G\)p\(\\mathcal\{D\}\\mid G\)for linear\-Gaussian models\[[13](https://arxiv.org/html/2608.12640#bib.bib7),[14](https://arxiv.org/html/2608.12640#bib.bib8)\]\. 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](https://arxiv.org/html/2608.12640#A1.SS1)\. GivenGG, the posterior over edge weights factorizes by node:

p⁡\(B∣G,𝒟\)=∏j=1dp⁡\(BpaG​\(Xj\),j∣G,𝒟\),p\(B\\mid G,\\mathcal\{D\}\)=\\prod\_\{j=1\}^\{d\}p\\left\(B\_\{\\mathrm\{pa\}\_\{G\}\(X\_\{j\}\),j\}\\mid G,\\mathcal\{D\}\\right\),where coefficients outside the parent set ofXjX\_\{j\}are fixed to zero and each nonzero incoming\-coefficient block follows a multivariatett\-distribution\[[36](https://arxiv.org/html/2608.12640#bib.bib6)\]\. 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\\mid\\mathcal\{D\}\)\\propto p\(G\),p\(\\mathcal\{D\}\\mid G\),p\(B\\mid G,\\mathcal\{D\}\)\.\(3\)We writeπ⁡\(Z\)=π⁡\(G,B\)=p⁡\(G,B∣𝒟\)\\pi\(Z\)=\\pi\(G,B\)=p\(G,B\\mid\\mathcal\{D\}\)for this unconstrained posterior in the method section for compactness\.

## 4Conditional 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 scoreh⁡\(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\.1Conditional causal discovery as a score\-level problem

LetZ=\(G,B\)Z=\(G,B\)denote a graph–weight state and letπ⁡\(Z\)=p⁡\(G,B∣𝒟\)\\pi\(Z\)=p\(G,B\\mid\\mathcal\{D\}\)be the unconstrained posterior from Eq\. equation[3](https://arxiv.org/html/2608.12640#S3.E3)\. Throughout this section, we useπ⁡\(𝒜\)\\pi\(\\mathcal\{A\}\)to denote the posterior mass of an event𝒜\\mathcal\{A\}\. A conditional query is specified by a score functionh:𝒵→ℝh:\\mathcal\{Z\}\\to\\mathbb\{R\}and a final levelλ⋆\\lambda\_\{\\star\}\. The score\-level event, its posterior probability, and the corresponding constrained posterior are

𝒜λ⋆\\displaystyle\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}=\{Z:h⁡\(Z\)≥λ⋆\},\\displaystyle=\\\{Z:h\(Z\)\\geq\\lambda\_\{\\star\}\\\},\(4\)p⋆\\displaystyle p\_\{\\star\}=π⁡\(𝒜λ⋆\),\\displaystyle=\\pi\(\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}\),π⋆​\(Z\)\\displaystyle\\pi\_\{\\star\}\(Z\)=π⁡\(Z∣𝒜λ⋆\)\.\\displaystyle=\\pi\\\!\\left\(Z\\mid\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}\\right\)\.This formulation separates the scientific query from the computational procedure: once the query has been written throughhhandλ⋆\\lambda\_\{\\star\}, 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 thresholdt\>0t\>0, the right\-tail eventℰi​j\+​\(t\)\\mathcal\{E\}\_\{ij\}^\{\+\}\(t\)is represented by choosingh​\(Z\)=CEi​j​\(Z\)h\(Z\)=\\mathrm\{CE\}\_\{ij\}\(Z\)andλ⋆=t\\lambda\_\{\\star\}=t\. The left\-tail eventℰi​j−​\(t\)\\mathcal\{E\}\_\{ij\}^\{\-\}\(t\)is represented by choosingh​\(Z\)=−CEi​j​\(Z\)h\(Z\)=\-\\mathrm\{CE\}\_\{ij\}\(Z\)andλ⋆=t\\lambda\_\{\\star\}=t\. Thus both tails are written as the same score\-level event𝒜λ⋆=\{Z:h⁡\(Z\)≥λ⋆\}\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}=\\\{Z:h\(Z\)\\geq\\lambda\_\{\\star\}\\\}, 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=1C\\mathcal\{C\}=\\\{\(i\_\{c\},j\_\{c\},\\bowtie\_\{c\},t\_\{c\},\\kappa\_\{c\}\)\\\}\_\{c=1\}^\{C\}be a collection ofCCconstraints, where⋈c∈\{≥,≤\}\\bowtie\_\{c\}\\in\\\{\\geq,\\leq\\\}gives the inequality direction,tc∈ℝt\_\{c\}\\in\\mathbb\{R\}is the threshold, andκc\>0\\kappa\_\{c\}\>0is an optional scale factor\. Defineδc=\+1\\delta\_\{c\}=\+1when⋈c=≥\\bowtie\_\{c\}=\\geqandδc=−1\\delta\_\{c\}=\-1when⋈c=≤\\bowtie\_\{c\}=\\leq\. The normalized margin of constraintccand the aggregate score are

mc​\(Z\)\\displaystyle m\_\{c\}\(Z\)=δc​\(CEic​jc​\(Z\)−tc\)κc,\\displaystyle=\\frac\{\\delta\_\{c\}\\bigl\(\\mathrm\{CE\}\_\{i\_\{c\}j\_\{c\}\}\(Z\)\-t\_\{c\}\\bigr\)\}\{\\kappa\_\{c\}\},\(5\)h𝒞​\(Z\)\\displaystyle h\_\{\\mathcal\{C\}\}\(Z\)=min1≤c≤C⁡mc​\(Z\)\.\\displaystyle=\\min\_\{1\\leq c\\leq C\}m\_\{c\}\(Z\)\.The marginmc​\(Z\)m\_\{c\}\(Z\)is nonnegative exactly when thecc\-th constraint is satisfied\. Thereforeh𝒞​\(Z\)≥0h\_\{\\mathcal\{C\}\}\(Z\)\\geq 0means that all constraints are satisfied, so the joint event isℰ𝒞=\{Z:h𝒞​\(Z\)≥0\}\\mathcal\{E\}\_\{\\mathcal\{C\}\}=\\\{Z:h\_\{\\mathcal\{C\}\}\(Z\)\\geq 0\\\}\. Interval constraints can be represented by including both a lower\-bound and an upper\-bound inequality\.

### 4\.2Adaptive multilevel splitting for rare posterior events

Extreme\-effect events can have very small posterior probability\. A direct posterior sampler withMMsamples will produce only aboutM​p⋆Mp\_\{\\star\}constrained samples on average, which may be close to zero ifp⋆p\_\{\\star\}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λ⋆\\lambda\_\{\\star\}\. 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=λ⋆\\lambda\_\{0\}<\\lambda\_\{1\}<\\cdots<\\lambda\_\{K\}=\\lambda\_\{\\star\}be increasing score levels, withλ0=−∞\\lambda\_\{0\}=\-\\inftycorresponding to the unconstrained posterior\. Define𝒜k=\{Z:h⁡\(Z\)≥λk\}\\mathcal\{A\}\_\{k\}=\\\{Z:h\(Z\)\\geq\\lambda\_\{k\}\\\}\. Because these events are nested, the rare\-event probability decomposes as

p⋆=π⁡\(𝒜K\)=∏k=0K−1π⁡\(𝒜k\+1∣𝒜k\)\.p\_\{\\star\}=\\pi\(\\mathcal\{A\}\_\{K\}\)=\\prod\_\{k=0\}^\{K\-1\}\\pi\(\\mathcal\{A\}\_\{k\+1\}\\mid\\mathcal\{A\}\_\{k\}\)\.\(6\)The advantage is that each factor in this product can be much larger thanp⋆p\_\{\\star\}itself, making it estimable with a moderate number of particles\.

We use the adaptive version of multilevel splitting\[[20](https://arxiv.org/html/2608.12640#bib.bib9),[18](https://arxiv.org/html/2608.12640#bib.bib10),[6](https://arxiv.org/html/2608.12640#bib.bib13)\], which chooses intermediate levels from the particle population rather than requiring them to be fixed in advance\. At levelkk, the particle population𝒫k=\{Zk,n\}n=1N\\mathcal\{P\}\_\{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\}\.\\pi\_\{\\lambda\_\{k\}\}\(Z\)\\propto\\pi\(Z\)\\mathbf\{1\}\\\{h\(Z\)\\geq\\lambda\_\{k\}\\\}\.\(7\)We then evaluate the scoresh⁡\(Zk,n\)h\(Z\_\{k,n\}\)and choose the next level as an empirical quantile:

λ~k\+1\\displaystyle\\widetilde\{\\lambda\}\_\{k\+1\}=Q1−ρ​\(\{h⁡\(Zk,n\)\}n=1N\),\\displaystyle=Q\_\{1\-\\rho\}\\\!\\left\(\\\{h\(Z\_\{k,n\}\)\\\}\_\{n=1\}^\{N\}\\right\),\(8\)λk\+1\\displaystyle\\lambda\_\{k\+1\}=min⁡\{λ~k\+1,λ⋆\}\.\\displaystyle=\\min\\\{\\widetilde\{\\lambda\}\_\{k\+1\},\\lambda\_\{\\star\}\\\}\.Hereρ∈\(0,1\)\\rho\\in\(0,1\)is the survival fraction andQ1−ρQ\_\{1\-\\rho\}is the empirical\(1−ρ\)\(1\-\\rho\)\-quantile\. The estimated conditional factor at this level is the fraction of particles that survive the new threshold,

β^k\+1=1N∑n=1N𝟏\{h\(Zk,n\)≥λk\+1\}\.\\widehat\{\\beta\}\_\{k\+1\}=\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\\mathbf\{1\}\\\{h\(Z\_\{k,n\}\)\\geq\\lambda\_\{k\+1\}\\\}\.\(9\)The survivors are resampled with replacement to form the starting population for the next level\. The final probability estimate isp^⋆=∏k=0K−1β^k\+1\\widehat\{p\}\_\{\\star\}=\\prod\_\{k=0\}^\{K\-1\}\\widehat\{\\beta\}\_\{k\+1\}, and after the final mutation step the particles approximate samples fromπ⋆​\(Z\)\\pi\_\{\\star\}\(Z\)\.

Algorithm 1Adaptive multilevel splitting for conditional causal discovery1:Posterior

π⁡\(Z\)\\pi\(Z\), score

hh, target level

λ⋆\\lambda\_\{\\star\}, initial particles

𝒫\\mathcal\{P\}, particle size

NN, survival fraction

ρ\\rho, mutation steps

mm, maximum levels

KmaxK\_\{\\max\}\.

2:Tail\-probability estimate

p^⋆\\widehat\{p\}\_\{\\star\}and constrained particles

𝒫\\mathcal\{P\}\.

3:Initialize

𝒫=\{Zn\}n=1N\\mathcal\{P\}=\\\{Z\_\{n\}\\\}\_\{n=1\}^\{N\}using an approximate posterior sampler or warm\-start routine targeting

π⁡\(Z\)\\pi\(Z\)\.

4:Set

λ0←−∞\\lambda\_\{0\}\\leftarrow\-\\inftyand

p^⋆←1\\widehat\{p\}\_\{\\star\}\\leftarrow 1\.

5:for

k=0,…,Kmax−1k=0,\\ldots,K\_\{\\max\}\-1do

6:

𝒫←Mutateλk​\(𝒫,m\)\\mathcal\{P\}\\leftarrow\\mathrm\{Mutate\}\_\{\\lambda\_\{k\}\}\(\\mathcal\{P\};m\)\.

7:Compute scores

sn←h⁡\(Zn\)s\_\{n\}\\leftarrow h\(Z\_\{n\}\)for all

Zn∈𝒫Z\_\{n\}\\in\\mathcal\{P\}\.

8:Set

λk\+1←min⁡\{Q1−ρ​\(s1,…,sN\),λ⋆\}\\lambda\_\{k\+1\}\\leftarrow\\min\\\{Q\_\{1\-\\rho\}\(s\_\{1\},\\ldots,s\_\{N\}\),\\lambda\_\{\\star\}\\\}\.

9:Let

𝒮←\{Zn∈𝒫:sn≥λk\+1\}\\mathcal\{S\}\\leftarrow\\\{Z\_\{n\}\\in\\mathcal\{P\}:s\_\{n\}\\geq\\lambda\_\{k\+1\}\\\}\.

10:Estimate

β^k\+1←\|𝒮\|/N\\widehat\{\\beta\}\_\{k\+1\}\\leftarrow\|\\mathcal\{S\}\|/Nand update

p^⋆←p^⋆​β^k\+1\\widehat\{p\}\_\{\\star\}\\leftarrow\\widehat\{p\}\_\{\\star\}\\widehat\{\\beta\}\_\{k\+1\}\.

11:if

\|𝒮\|=0\|\\mathcal\{S\}\|=0then

12:return

0,∅0,\\emptyset\.

13:endif

14:Resample

NNparticles from

𝒮\\mathcal\{S\}with replacement to form

𝒫\\mathcal\{P\}\.

15:if

λk\+1=λ⋆\\lambda\_\{k\+1\}=\\lambda\_\{\\star\}then

16:

𝒫←Mutateλ⋆​\(𝒫,m\)\\mathcal\{P\}\\leftarrow\\mathrm\{Mutate\}\_\{\\lambda\_\{\\star\}\}\(\\mathcal\{P\};m\)\.

17:return

p^⋆,𝒫\\widehat\{p\}\_\{\\star\},\\mathcal\{P\}\.

18:endif

19:endfor

20:return

p^⋆,𝒫\\widehat\{p\}\_\{\\star\},\\mathcal\{P\}\.

21:

22:HereMutateλ​\(𝒫,m\)\\mathrm\{Mutate\}\_\{\\lambda\}\(\\mathcal\{P\};m\)denotesmmMetropolis–Hastings

23:steps per particle targetingπλ\(Z\)∝π\(Z\)𝟏\{h\(Z\)≥λ\}\\pi\_\{\\lambda\}\(Z\)\\propto\\pi\(Z\)\\mathbf\{1\}\\\{h\(Z\)\\geq\\lambda\\\}\.

### 4\.3Particle 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\)\\pi\(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𝒜λ⋆\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}\.

In our implementation, we use a simple warm\-start scheme\. We first draw sparse DAGsG0,nG\_\{0,n\}from the structural priorp⁡\(G\)p\(G\), using an Erdős–Rényi DAG prior in the experiments with edge inclusion probabilitypedgep\_\{\\mathrm\{edge\}\}chosen to match a target expected number of edges per node\. Given each initial graph, we draw weights from the conjugate posteriorB0,n∼p⁡\(B∣G0,n,𝒟\)B\_\{0,n\}\\sim p\(B\\mid G\_\{0,n\},\\mathcal\{D\}\)and setZ0,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\)\\pi\(Z\), we use the inner MCMC kernel at the initial levelλ0=−∞\\lambda\_\{0\}=\-\\inftyas 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\.4Inner Metropolis–Hastings over joint\(G,B\)\(G,B\)

At each score levelλk\\lambda\_\{k\}, mutation must preserve the level\-truncated targetπλk​\(Z\)\\pi\_\{\\lambda\_\{k\}\}\(Z\)from Eq\. equation[7](https://arxiv.org/html/2608.12640#S4.E7)\. We use a blocked Metropolis–Hastings kernel over the joint stateZ=\(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 onlyGGcan leave the sampler stuck at a fixed set of coefficients, while changing onlyBBcannot explore alternative pathways\.

With probabilitypstructp\_\{\\mathrm\{struct\}\}, we propose a structure move\. The new graphG′G^\{\\prime\}is drawn from a proposalqG\(⋅∣G\)q\_\{G\}\(\\cdot\\mid G\), instantiated either as Structure\-MCMC, which adds, deletes, or reverses a single edge while rejecting cyclic graphs, or as PARNI\-DAG\[[24](https://arxiv.org/html/2608.12640#bib.bib23)\], which uses locally informed adaptive neighborhoods guided by posterior edge information\. After proposingG′G^\{\\prime\}, we refresh only the coefficient blocks whose parent sets changed\. For each affected nodevv, we sample

BpaG′​\(v\),v′∼p⁡\(BpaG′​\(v\),v∣G′,𝒟\),B^\{\\prime\}\_\{\\mathrm\{pa\}\_\{G^\{\\prime\}\}\(v\),v\}\\sim p\\\!\\left\(B\_\{\\mathrm\{pa\}\_\{G^\{\\prime\}\}\(v\),v\}\\mid G^\{\\prime\},\\mathcal\{D\}\\right\),\(10\)set coefficients for absent edges to zero, and copy all unchanged blocks fromBB\.

With probability1−pstruct1\-p\_\{\\mathrm\{struct\}\}, we propose a weight move at fixed structure\. We keepG′=GG^\{\\prime\}=G, select a nodevvwith at least one parent, and resample its incoming coefficient block,

BpaG​\(v\),v′∼p⁡\(BpaG​\(v\),v∣G,𝒟\),B^\{\\prime\}\_\{\\mathrm\{pa\}\_\{G\}\(v\),v\}\\sim p\\\!\\left\(B\_\{\\mathrm\{pa\}\_\{G\}\(v\),v\}\\mid G,\\mathcal\{D\}\\right\),\(11\)leaving all other entries unchanged\. These blocked refreshes are inexpensive under the conjugate linear\-Gaussian model because each conditional coefficient posterior is a multivariatettdistribution\.

Letq⁡\(Z′∣Z\)q\(Z^\{\\prime\}\\mid 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 satisfiesh⁡\(Z\)≥λkh\(Z\)\\geq\\lambda\_\{k\}, the level\-λk\\lambda\_\{k\}acceptance probability is

aλk​\(Z,Z′\)\\displaystyle a\_\{\\lambda\_\{k\}\}\(Z,Z^\{\\prime\}\)=𝟏\{h\(Z′\)≥λk\}\\displaystyle=\\mathbf\{1\}\\\{h\(Z^\{\\prime\}\)\\geq\\lambda\_\{k\}\\\}\(12\)×min⁡\{1,π⁡\(Z′\)​q​\(Z∣Z′\)π⁡\(Z\)​q​\(Z′∣Z\)\}\.\\displaystyle\\times\\min\\\!\\left\\\{1,\\frac\{\\pi\(Z^\{\\prime\}\)q\(Z\\mid Z^\{\\prime\}\)\}\{\\pi\(Z\)q\(Z^\{\\prime\}\\mid Z\)\}\\right\\\}\.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

## 5Experiments

### 5\.1Experiment Setup

We evaluate the proposed conditional causal discovery procedure on synthetic linear\-Gaussian datasets generated from random Erdős–Rényi DAGs withd∈\{4,8,16,32\}d\\in\\\{4,8,16,32\\\}\. For each dimension, we first sample a data\-generating DAGG⋆G^\{\\star\}using an Erdős–Rényi DAG generator with target sparsity approximately2​d2dedges\. The generator samples an acyclic ordering, draws only order\-compatible directed edges, and returns an adjacency matrixG⋆G^\{\\star\}\. The exact edge\-sampling rule, including the edge probability or edge budget used for eachddand the small\-graph handling ford=4d=4, is provided in Appendix[A\.1](https://arxiv.org/html/2608.12640#A1.SS1)\.

GivenG⋆G^\{\\star\}, we draw raw edge weights independently asB~i​j∼𝒩⁡\(0,1\)\\widetilde\{B\}\_\{ij\}\\sim\\mathcal\{N\}\(0,1\)and mask them by the sampled adjacency, so thatB⋆=B~⊙G⋆B^\{\\star\}=\\widetilde\{B\}\\odot G^\{\\star\}\. Equivalently,Bi​j⋆=B~i​jB^\{\\star\}\_\{ij\}=\\widetilde\{B\}\_\{ij\}ifGi​j⋆=1G^\{\\star\}\_\{ij\}=1andBi​j⋆=0B^\{\\star\}\_\{ij\}=0otherwise\. We then generate an observational dataset𝒟\\mathcal\{D\}of sizenobs=1000n\_\{\\mathrm\{obs\}\}=1000from the linear\-Gaussian SEM specified by\(G⋆,B⋆\)\(G^\{\\star\},B^\{\\star\}\)\.

We evaluate our multilevel splitting framework instantiated with either the Structure\-MCMC kernel or the PARNI\-DAG kernel\[[24](https://arxiv.org/html/2608.12640#bib.bib23)\]\. We compare against four baselines: \(i\) exhaustive enumeration, which is a gold\-standard baseline feasible ford=4d=4; \(ii\) DiBS\[[25](https://arxiv.org/html/2608.12640#bib.bib1)\]; \(iii\) OrderSPN\[[37](https://arxiv.org/html/2608.12640#bib.bib11)\]; 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](https://arxiv.org/html/2608.12640#A1.SS1)\.

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\.2Correctness 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\.1Single\-effect conditioning: validation ond=4d=4

We construct 10 test cases by sampling 10 independentd=4d=4graphs 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\)≥λ⋆\}\\mathcal\{A\}\_\{\\lambda\_\{\\star\}\}=\\\{Z:h\(Z\)\\geq\\lambda\_\{\\star\}\\\}from Sec\.[4](https://arxiv.org/html/2608.12640#S4)\. For every test case, we run each stochastic method 10 times and report the mean estimate together with across\-run variability\.

![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/7model_comparison_d=4_with_errorbar_improved.png)Figure 2:Single\-effect one\-sided tail\-probability estimation ond=4d=4across 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](https://arxiv.org/html/2608.12640#A1.T3)\.Figure[2](https://arxiv.org/html/2608.12640#S5.F2)compares all methods against exhaustive enumeration\. We make two observations\. First, all methods except DiBS are broadly consistent with the enumeration baseline across thed=4d=4cases, with only minor deviations; the corresponding numerical values are reported in Table[3](https://arxiv.org/html/2608.12640#A1.T3)\. 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\.2Multi\-effect conditioning: validating the conditioning operator ond=4d=4andd=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 selectCCordered node pairs\{\(ic,jc\)\}c=1C\\\{\(i\_\{c\},j\_\{c\}\)\\\}\_\{c=1\}^\{C\}\. Letec⋆=CEic​jc​\(Z⋆\)e\_\{c\}^\{\\star\}=\\mathrm\{CE\}\_\{i\_\{c\}j\_\{c\}\}\(Z^\{\\star\}\)denote the ground\-truth causal effect for pair\(ic,jc\)\(i\_\{c\},j\_\{c\}\)under the data\-generating stateZ⋆=\(G⋆,B⋆\)Z^\{\\star\}=\(G^\{\\star\},B^\{\\star\}\)\. For each pair, we impose a two\-sided interval constraint aroundec⋆e\_\{c\}^\{\\star\}:

CEic​jc\(Z\)∈\[ec⋆−ε,ec⋆\+ε\],c=1,…,C\.\\mathrm\{CE\}\_\{i\_\{c\}j\_\{c\}\}\(Z\)\\in\[e\_\{c\}^\{\\star\}\-\\varepsilon,\\,e\_\{c\}^\{\\star\}\+\\varepsilon\],\\quad c=1,\\ldots,C\.\(13\)Each interval is represented by two one\-sided inequalities, and the conjunction of all inequalities defines the conditioning event\. Following Sec\.[4](https://arxiv.org/html/2608.12640#S4), we aggregate these inequalities using the minimum normalized margin scoreh𝒞​\(Z\)h\_\{\\mathcal\{C\}\}\(Z\), so that the target event isℰ𝒞=\{Z:h𝒞​\(Z\)≥0\}\\mathcal\{E\}\_\{\\mathcal\{C\}\}=\\\{Z:h\_\{\\mathcal\{C\}\}\(Z\)\\geq 0\\\}\.

We compare aweakconstraint set, with fewer constrained pairs and hence a looser event, against astrongconstraint 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 graphG⋆G^\{\\star\}\.

Table 1:d=4d=4\(PARNI\)structural accuracy under multi\-effect conditioning\.ConstraintsRunsSamples/runMean SHD↓\\downarrowPr⁡\(SHD=0\)↑\\Pr\(\\mathrm\{SHD\}=0\)\\uparrowWeak10200\.325±0\.1640\.325\\pm 0\.1640\.790±0\.1170\.790\\pm 0\.117Strong10200\.115±0\.0530\.115\\pm 0\.0530\.885±0\.0530\.885\\pm 0\.053

Table 2:d=8d=8\(PARNI\)structural accuracy under multi\-effect conditioning\. We report mean±\\pmstandard deviation of the per\-run mean SHD, averaged over sampled DAGs within each run, and the fraction of samples withSHD=0\\mathrm\{SHD\}=0\.ConstraintsRunsSamples/runMean SHD↓\\downarrowPr⁡\(SHD=0\)↑\\Pr\(\\mathrm\{SHD\}=0\)\\uparrowWeak10505\.342±0\.4425\.342\\pm 0\.4420\.040±0\.0270\.040\\pm 0\.027Strong10502\.290±0\.1992\.290\\pm 0\.1990\.108±0\.0660\.108\\pm 0\.066

##### Multi\-effect conditioning results\.

As demonstrated by Figs\.[6](https://arxiv.org/html/2608.12640#A1.F6)and[7](https://arxiv.org/html/2608.12640#A1.F7)in Appendix[A\.7](https://arxiv.org/html/2608.12640#A1.SS7), across bothd=4d=4andd=8d=8, stronger, more informative multi\-effect constraints consistently concentrate the conditional posterior toward the ground\-truth mechanism\. Quantitatively, Tables[1](https://arxiv.org/html/2608.12640#S5.T1)and[2](https://arxiv.org/html/2608.12640#S5.T2)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\.3Scalability

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 ond∈\{8,16,32\}d\\in\\\{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λ⋆\\lambda\_\{\\star\}becomes more stringent; whether independent runs give reproducible tail curves; and whether informed structure proposals become more important as the graph dimension increases\.

![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/single_effect_scalability_d8_d16_d32.png)Figure 3:Single\-effect one\-sided tail\-probability estimation ond∈\{8,16,32\}d\\in\\\{8,16,32\\\}\.Figure[3](https://arxiv.org/html/2608.12640#S5.F3)summarizes one\-sided single\-effect tail\-probability estimation ond∈\{8,16,32\}d\\in\\\{8,16,32\\\}\. Ond=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λ⋆\\lambda\_\{\\star\}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λ⋆\\lambda\_\{\\star\}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 enumerabled=4d=4setting\.

Accordingly, ford=16d=16andd=32d=32, we primarily assess robustness via internal consistency, because exact gold\-standard posterior enumeration is unavailable at these scales\. Ford=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\. Atd=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\.

## 6Case Study

We now study an application of our conditional causal discovery framework to a real world protein\-signalling dataset\[[31](https://arxiv.org/html/2608.12640#bib.bib12)\]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, containingn=7466n=7466measurements ofd=11d=11proteins\[[31](https://arxiv.org/html/2608.12640#bib.bib12)\]\. As a qualitative reference, we use the 20\-edge consensus network in Fig\.[4](https://arxiv.org/html/2608.12640#A1.F4)\[[21](https://arxiv.org/html/2608.12640#bib.bib25),[31](https://arxiv.org/html/2608.12640#bib.bib12)\]\. We study two target pairs:

\(i\)PIP3→PIP2\\mathrm\{PIP3\}\\rightarrow\\mathrm\{PIP2\}\(indices6→56\\to 5\), which admits a clean SEM decomposition into a direct routePIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}and a mediated routePIP3→Plcg→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{Plcg\}\\to\\mathrm\{PIP2\};

\(ii\)Erk→Akt\\mathrm\{Erk\}\\rightarrow\\mathrm\{Akt\}\(indices1→101\\to 10\), a widely discussed cross\-module influence in the Sachs signaling system\.111Node order used throughout: \[Raf,Erk,Plcg,PKC,PKA,PIP2,PIP3,Mek,P38,Jnk,Akt\]\. In particular,PIP3\\mathrm\{PIP3\}is index 6,PIP2\\mathrm\{PIP2\}is index 5,Plcg\\mathrm\{Plcg\}is index 2,Erk\\mathrm\{Erk\}is index 1, andAkt\\mathrm\{Akt\}is index 10\.

A published SEM mediation example on the same dataset reports a total effect of approximately0\.56650\.5665forPIP3→PIP2\\mathrm\{PIP3\}\\rightarrow\\mathrm\{PIP2\}, with an explicit direct/indirect decomposition viaPlcg\\mathrm\{Plcg\}, which we use as an external numerical anchor for the effect scale\[[26](https://arxiv.org/html/2608.12640#bib.bib24)\]\. 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](https://arxiv.org/html/2608.12640#S3), letZ=\(G,B\)Z=\(G,B\)denote a graph–weight state sampled from the posteriorπ⁡\(Z\)=p⁡\(G,B∣𝒟\)\\pi\(Z\)=p\(G,B\\mid\\mathcal\{D\}\), and letCEi​j​\(Z\)\\mathrm\{CE\}\_\{ij\}\(Z\)denote the signed linear\-SEM total effect\. For the two target pairs, define the eventsℰPIP​\(t\)=\{Z:CEPIP3,PIP2​\(Z\)≥t\}\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t\)=\\left\\\{Z:\\mathrm\{CE\}\_\{\\mathrm\{PIP3\},\\mathrm\{PIP2\}\}\(Z\)\\geq t\\right\\\}, andℰERK​\(t\)=\{Z:CEErk,Akt​\(Z\)≥t\}\\mathcal\{E\}\_\{\\mathrm\{ERK\}\}\(t\)=\\left\\\{Z:\\mathrm\{CE\}\_\{\\mathrm\{Erk\},\\mathrm\{Akt\}\}\(Z\)\\geq t\\right\\\}\. We compare four posterior conditions:

- •Unconditioned:Z∼π⁡\(Z\)=p⁡\(G,B∣𝒟\)Z\\sim\\pi\(Z\)=p\(G,B\\mid\\mathcal\{D\}\)\.
- •Cond\-PIP:Z∼p⁡\(Z∣𝒟,ℰPIP​\(tPIP\)\)Z\\sim p\(Z\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t\_\{\\mathrm\{PIP\}\}\)\), withtPIP=0\.77t\_\{\\mathrm\{PIP\}\}=0\.77\.
- •Cond\-ERK:Z∼p⁡\(Z∣𝒟,ℰERK​\(tERK\)\)Z\\sim p\(Z\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{\\mathrm\{ERK\}\}\(t\_\{\\mathrm\{ERK\}\}\)\), withtERK=0\.67t\_\{\\mathrm\{ERK\}\}=0\.67\.
- •Cond\-Joint: Z∼p⁡\(Z∣𝒟,ℰPIP​\(tPIP′\)∩ℰERK​\(tERK′\)\),Z\\sim p\\\!\\left\(Z\\mid\\mathcal\{D\},\\,\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t^\{\\prime\}\_\{\\mathrm\{PIP\}\}\)\\cap\\mathcal\{E\}\_\{\\mathrm\{ERK\}\}\(t^\{\\prime\}\_\{\\mathrm\{ERK\}\}\)\\right\),with\(tPIP′,tERK′\)=\(0\.74,0\.65\)\(t^\{\\prime\}\_\{\\mathrm\{PIP\}\},t^\{\\prime\}\_\{\\mathrm\{ERK\}\}\)=\(0\.74,0\.65\), targeting a similarly rare posterior region\.

The estimated posterior masses of these events are reported in Appendix[A\.6](https://arxiv.org/html/2608.12640#A1.SS6)and Table[6](https://arxiv.org/html/2608.12640#A1.T6)\. These conditioning events do*not*assert that the dataset corresponds to a single “ground\-truth” extreme state\. Instead, they define posterior queries of the formp⁡\(Z∣𝒟,ℰ\)p\(Z\\mid\\mathcal\{D\},\\mathcal\{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](https://arxiv.org/html/2608.12640#A1.SS1)\. Table[7](https://arxiv.org/html/2608.12640#A1.T7)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](https://arxiv.org/html/2608.12640#bib.bib12),[21](https://arxiv.org/html/2608.12640#bib.bib25),[12](https://arxiv.org/html/2608.12640#bib.bib14)\]\.

##### Mechanistic analysis under conditioning\.

- •Cond\-PIP:The conditional posterior yields a highly concentrated explanation forPIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}: the direct route and thePlcg\\mathrm\{Plcg\}\-mediated route appear in essentially all conditional samples, and the conditional mean effect is approximately0\.7760\.776\(calculated in Sec\.[A\.6](https://arxiv.org/html/2608.12640#A1.SS6)\), with an increasedPlcg\\mathrm\{Plcg\}\-mediated share \(Table[8](https://arxiv.org/html/2608.12640#A1.T8)\)\.
- •Cond\-ERK:Erk→Akt\\mathrm\{Erk\}\\to\\mathrm\{Akt\}remains direct\-dominated, whilePIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}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](https://arxiv.org/html/2608.12640#A1.T7)and[8](https://arxiv.org/html/2608.12640#A1.T8)\)\.
- •Non\-symmetry:Conditioning onPIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}does not forceErk→Akt\\mathrm\{Erk\}\\to\\mathrm\{Akt\}to become extreme; for example,Pr⁡\(CEErk,Akt​\(Z\)\>0∣𝒟,ℰPIP​\(tPIP\)\)≈0\.498\\Pr\\\!\\left\(\\mathrm\{CE\}\_\{\\mathrm\{Erk\},\\mathrm\{Akt\}\}\(Z\)\>0\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t\_\{\\mathrm\{PIP\}\}\)\\right\)\\approx 0\.498\.
- •Cond\-Joint:Joint conditioning increases the frequency of the mediatedErk→Plcg→Akt\\mathrm\{Erk\}\\to\\mathrm\{Plcg\}\\to\\mathrm\{Akt\}route relative to Cond\-ERK, highlightingPlcg\\mathrm\{Plcg\}as a shared mediator under co\-extreme coupling, beyond what either single\-condition run reveals alone \(Table[8](https://arxiv.org/html/2608.12640#A1.T8)\)\.

##### Graph\-level shifts under conditioning\.

Figure[5](https://arxiv.org/html/2608.12640#A1.F5)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\-PIPexplainsPIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}mainly through two routes, whereasCond\-ERKkeepsErk→Akt\\mathrm\{Erk\}\\to\\mathrm\{Akt\}mostly direct but yields a broader set ofPIP3→PIP2\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}paths with partial cancellation \(Table[8](https://arxiv.org/html/2608.12640#A1.T8)\)\. Joint conditioning is not just the overlap of the two single\-effect results: compared toCond\-ERK, it more often highlights the mediated routeErk→Plcg→Akt\\mathrm\{Erk\}\\to\\mathrm\{Plcg\}\\to\\mathrm\{Akt\}, pointing toPlcg\\mathrm\{Plcg\}as a shared mediator that emerges under the joint extreme\-effect query \(Table[8](https://arxiv.org/html/2608.12640#A1.T8)\)\.

## 7Conclusions

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

- Anandet al\.\(2023\)T\. V\. Anand, A\. H\. Ribeiro, J\. Tian, and E\. BareinboimCausal effect identification in cluster dags\.InProceedings of the AAAI Conference on Artificial Intelligence,Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Annadaniet al\.\(2021\)Y\. Annadani, J\. Rothfuss, A\. Lacoste, N\. Scherrer, A\. Goyal, Y\. Bengio, and S\. BauerVariational causal networks: approximate bayesian inference over causal structures\.arXiv preprint arXiv:2106\.07635\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Bodiket al\.\(2024\)J\. Bodik, M\. Paluš, and Z\. PawlasCausality in extremes of time series\.Extremes27\(1\),pp\. 67–121\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p2.1)\.
- Borboudakis and Tsamardinos \(2012\)G\. Borboudakis and I\. TsamardinosIncorporating causal prior knowledge as path\-constraints in bayesian networks and maximal ancestral graphs\.InProceedings of the 29th International Coference on International Conference on Machine Learning,pp\. 427–434\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Brouillardet al\.\(2020\)P\. Brouillard, S\. Lachapelle, A\. Lacoste, S\. Lacoste\-Julien, and A\. DrouinDifferentiable causal discovery from interventional data\.Advances in Neural Information Processing Systems33,pp\. 21865–21877\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Cérou and Guyader \(2007\)F\. Cérou and A\. GuyaderAdaptive multilevel splitting for rare event analysis\.Stochastic Analysis and Applications25\(2\),pp\. 417–443\.Cited by:[§1](https://arxiv.org/html/2608.12640#S1.p4.1),[§2](https://arxiv.org/html/2608.12640#S2.p5.1),[§4\.2](https://arxiv.org/html/2608.12640#S4.SS2.p3.2)\.
- Chenet al\.\(2016\)E\. Y\. Chen, Y\. Shen, A\. Choi, and A\. DarwicheLearning bayesian networks with ancestral constraints\.Advances in Neural Information Processing Systems29\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Chen and Darwiche \(2024\)Y\. Chen and A\. DarwicheIdentifying causal effects under functional dependencies\.Advances in Neural Information Processing Systems37,pp\. 87013–87041\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Cundyet al\.\(2021\)C\. Cundy, A\. Grover, and S\. ErmonBcd nets: scalable variational approaches for bayesian causal discovery\.Advances in Neural Information Processing Systems34,pp\. 7095–7110\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Deleuet al\.\(2022\)T\. Deleu, A\. Góis, C\. Emezue, M\. Rankawat, S\. Lacoste\-Julien, S\. Bauer, and Y\. BengioBayesian structure learning with generative flow networks\.InUncertainty in Artificial Intelligence,pp\. 518–528\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Emezueet al\.\(2023\)C\. C\. Emezue, A\. Drouin, T\. Deleu, S\. Bauer, and Y\. BengioBenchmarking bayesian causal discovery methods for downstream treatment effect estimation\.arXiv preprint arXiv:2307\.04988\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p3.1)\.
- Friedman and Koller \(2003\)N\. Friedman and D\. KollerBeing bayesian about network structure\. a bayesian approach to structure discovery in bayesian networks\.Machine learning50,pp\. 95–125\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p3.1),[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px3.p1.1)\.
- Geiger and Heckerman \(1994\)D\. Geiger and D\. HeckermanLearning gaussian networks\.InUncertainty in Artificial Intelligence,pp\. 235–243\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px2.p2.1)\.
- Geiger and Heckerman \(2002\)D\. Geiger and D\. HeckermanParameter priors for directed acyclic graphical models and the characterization of several probability distributions\.The Annals of Statistics30\(5\),pp\. 1412–1440\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px2.p2.1)\.
- Giudiceet al\.\(2023\)E\. Giudice, J\. Kuipers, and G\. MoffaA bayesian take on gaussian process networks\.Advances in Neural Information Processing Systems36,pp\. 56602–56614\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p3.1)\.
- Glymouret al\.\(2019\)C\. Glymour, K\. Zhang, and P\. SpirtesReview of causal discovery methods based on graphical models\.Frontiers in genetics10,pp\. 524\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px2.p1.1)\.
- Gneccoet al\.\(2021\)N\. Gnecco, N\. Meinshausen, J\. Peters, and S\. EngelkeCausal discovery in heavy\-tailed models\.The Annals of Statistics49\(3\),pp\. 1755–1778\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p2.1)\.
- Guyaderet al\.\(2011\)A\. Guyader, N\. Hengartner, and E\. Matzner\-LøberSimulation and estimation of extreme quantiles and extreme probabilities\.Applied Mathematics & Optimization64\(2\),pp\. 171–196\.Cited by:[§4\.2](https://arxiv.org/html/2608.12640#S4.SS2.p3.2)\.
- Hauser and Bühlmann \(2015\)A\. Hauser and P\. BühlmannJointly interventional and observational data: estimation of interventional markov equivalence classes of directed acyclic graphs\.Journal of the Royal Statistical Society Series B: Statistical Methodology77\(1\),pp\. 291–318\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Kahn and Harris \(1951\)H\. Kahn and T\. E\. HarrisEstimation of particle transmission by random sampling\.National Bureau of Standards applied mathematics series12,pp\. 27–30\.Cited by:[§4\.2](https://arxiv.org/html/2608.12640#S4.SS2.p3.2)\.
- Kochet al\.\(2009\)M\. Koch, B\. M\. Broom, and D\. SubramanianLearning robust cell signalling models from high throughput proteomic data\.International journal of bioinformatics research and applications5\(3\),pp\. 241–253\.Cited by:[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px1.p1.1),[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px3.p1.1)\.
- Koller and Friedman \(2009\)D\. Koller and N\. FriedmanProbabilistic graphical models: principles and techniques\.MIT press\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px2.p1.1)\.
- Kuipers and Moffa \(2017\)J\. Kuipers and G\. MoffaPartition mcmc for inference on acyclic digraphs\.Journal of the American Statistical Association112\(517\),pp\. 282–299\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p3.1)\.
- Lianget al\.\(2023\)X\. Liang, A\. Caron, S\. Livingstone, and J\. GriffinStructure learning with adaptive random neighborhood informed mcmc\.Advances in Neural Information Processing Systems36,pp\. 40760–40772\.Cited by:[item 4](https://arxiv.org/html/2608.12640#A1.I1.i4.p1.1.1),[item 6](https://arxiv.org/html/2608.12640#A1.I1.i6.p1.1.1),[§2](https://arxiv.org/html/2608.12640#S2.p3.1),[§4\.4](https://arxiv.org/html/2608.12640#S4.SS4.p2.1),[§5\.1](https://arxiv.org/html/2608.12640#S5.SS1.p3.1)\.
- Lorchet al\.\(2021\)L\. Lorch, J\. Rothfuss, B\. Schölkopf, and A\. KrauseDibs: differentiable bayesian structure learning\.Advances in Neural Information Processing Systems34,pp\. 24111–24123\.Cited by:[item 3](https://arxiv.org/html/2608.12640#A1.I1.i3.p1.1.1),[§2](https://arxiv.org/html/2608.12640#S2.p4.1),[§5\.1](https://arxiv.org/html/2608.12640#S5.SS1.p3.1)\.
- Madhanagopal and Amrhein \(2019\)B\. Madhanagopal and J\. AmrheinAnalyzing structural causal models using the calis procedure\.InSAS Global Forum,pp\. 3240–2019\.Cited by:[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px1.p4.1)\.
- Meek \(1995\)C\. MeekCausal inference and causal explanation with background knowledge\.InProceedings of the Eleventh conference on Uncertainty in artificial intelligence,pp\. 403–410\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p1.1)\.
- Pascheet al\.\(2023\)O\. C\. Pasche, V\. Chavez\-Demoulin, and A\. C\. DavisonCausal modelling of heavy\-tailed variables and confounders with application to river flow\.Extremes26\(3\),pp\. 573–594\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p2.1)\.
- Pearl \(2009\)J\. PearlCausality: models, reasoning and inference\.2nd edition,Cambridge University Press,USA\.External Links:ISBN 052189560XCited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px1.p3.1)\.
- Rittel and Tschiatschek \(2023\)S\. Rittel and S\. TschiatschekSpecifying prior beliefs over dags in deep bayesian causal structure learning\.InECAI 2023: 26th European Conference on Artificial Intelligence, including 12th Conference on Prestigious Applications of Intelligent Systems, PAIS 2023\-Proceedings,pp\. 1962–1969\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Sachset al\.\(2005\)K\. Sachs, O\. Perez, D\. Pe’er, D\. A\. Lauffenburger, and G\. P\. NolanCausal protein\-signaling networks derived from multiparameter single\-cell data\.Science308\(5721\),pp\. 523–529\.Cited by:[§1](https://arxiv.org/html/2608.12640#S1.p1.1),[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px1.p1.1),[§6](https://arxiv.org/html/2608.12640#S6.SS0.SSS0.Px3.p1.1),[§6](https://arxiv.org/html/2608.12640#S6.p1.1)\.
- Sobel \(1990\)M\. E\. SobelEffect analysis and causation in linear structural equation models\.Psychometrika55\(3\),pp\. 495–515\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px1.p3.2)\.
- Spirteset al\.\(2000\)P\. Spirtes, C\. N\. Glymour, R\. Scheines, and D\. HeckermanCausation, prediction, and search\.MIT press\.Cited by:[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px1.p3.1)\.
- Tothet al\.\(2024\)C\. Toth, C\. Knoll, F\. Pernkopf, and R\. PeharzEffective bayesian causal inference via structural marginalisation and autoregressive orders\.arXiv preprint arXiv:2402\.14781\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Tothet al\.\(2022\)C\. Toth, L\. Lorch, C\. Knoll, A\. Krause, F\. Pernkopf, R\. Peharz, and J\. Von KügelgenActive bayesian causal inference\.Advances in Neural Information Processing Systems35,pp\. 16261–16275\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p4.1)\.
- Viinikkaet al\.\(2020\)J\. Viinikka, A\. Hyttinen, J\. Pensar, and M\. KoivistoTowards scalable bayesian learning of causal dags\.Advances in Neural Information Processing Systems33,pp\. 6584–6594\.Cited by:[§2](https://arxiv.org/html/2608.12640#S2.p3.1),[§3](https://arxiv.org/html/2608.12640#S3.SS0.SSS0.Px2.p2.2)\.
- Wanget al\.\(2022\)B\. Wang, M\. R\. Wicker, and M\. KwiatkowskaTractable uncertainty for structure learning\.InInternational Conference on Machine Learning,pp\. 23131–23150\.Cited by:[item 2](https://arxiv.org/html/2608.12640#A1.I1.i2.p1.1.1),[§2](https://arxiv.org/html/2608.12640#S2.p4.1),[§5\.1](https://arxiv.org/html/2608.12640#S5.SS1.p3.1)\.

## Appendix AAppendix

### A\.1Experiment Setting

For all synthetic experiments, we generate the ground\-truth DAGG⋆G^\{\\star\}using an ordered Erdős–Rényi construction\. Given dimensionddand edge budget parameteredges\_per\_node=2\\texttt\{edges\\\_per\\\_node\}=2, we set

pedge=min⁡\{edges\_per\_node⋅dd⁡\(d−1\)/2,0\.5\}\.p\_\{\\mathrm\{edge\}\}=\\min\\\!\\left\\\{\\frac\{\\texttt\{edges\\\_per\\\_node\}\\cdot d\}\{d\(d\-1\)/2\},\\,0\.5\\right\\\}\.We then sample a Bernoulli adjacency matrix with edge probabilitypedgep\_\{\\mathrm\{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 thereforepedge=0\.5p\_\{\\mathrm\{edge\}\}=0\.5ford=4d=4,pedge=0\.5p\_\{\\mathrm\{edge\}\}=0\.5ford=8d=8,pedge=4/15p\_\{\\mathrm\{edge\}\}=4/15ford=16d=16, andpedge=4/31p\_\{\\mathrm\{edge\}\}=4/31ford=32d=32\. The cap at0\.50\.5is the small\-graph handling used by the generator: ford=4d=4, the nominal value4/34/3is clipped to0\.50\.5, giving an expected33edges instead of an overly dense four\-node DAG\. The corresponding expected edge counts are33,1414,3232, and6464ford=4,8,16,32d=4,8,16,32, respectively\.

Conditional onG⋆G^\{\\star\}, 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 usingnobs=1000n\_\{\\mathrm\{obs\}\}=1000training samples\.

Below are the experiment settings for each method:

1. 1\.Exhaustive enumeration \(d=4d=4only\):enumerate all DAGs \(\|𝒢d=4\|=543\|\\mathcal\{G\}\_\{d=4\}\|=543\) and draw2000020000weight samples per graph\.
2. 2\.OrderSPN\[[37](https://arxiv.org/html/2608.12640#bib.bib11)\]:sample50005000graphs from the learned circuit and draw100100weight samples per graph\.
3. 3\.DiBS\[[25](https://arxiv.org/html/2608.12640#bib.bib1)\]:draw500000500000joint graph–weight samples in total\.
4. 4\.PARNI\-DAG single chain\[[24](https://arxiv.org/html/2608.12640#bib.bib23)\]:run a single long MCMC chain over\(G,B\)\(G,B\)using PARNI\-DAG as the structure proposal kernel for500000500000iterations with10%10\\%burn\-in\.
5. 5\.Structure\-MCMC single chain:run a single long MCMC chain over\(G,B\)\(G,B\)using Structure\-MCMC as the structure proposal kernel for500000500000iterations with10%10\\%burn\-in\.
6. 6\.MLS\-PARNI\-DAG\[[24](https://arxiv.org/html/2608.12640#bib.bib23)\]:use our MLS framework with PARNI\-DAG as the structure proposal kernel\. We useN=200N=200particles,m=2000m=2000MCMC iterations per level, andKmax=10K\_\{\\max\}=10levels for all dimensions\.
7. 7\.MLS\-Structure\-MCMC:use our MLS framework with Structure\-MCMC as the structure proposal kernel, with the same MLS hyperparameters as above\.

### A\.2Target\-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=CEi​j​\(Z⋆\)=CEi​j​\(G⋆,B⋆\)e\_\{0\}=\\mathrm\{CE\}\_\{ij\}\(Z^\{\\star\}\)=\\mathrm\{CE\}\_\{ij\}\(G^\{\\star\},B^\{\\star\}\)denote the ground\-truth signed total causal effect for a queried ordered pair\(i,j\)\(i,j\), whereZ⋆=\(G⋆,B⋆\)Z^\{\\star\}=\(G^\{\\star\},B^\{\\star\}\)is the data\-generating graph–weight state\.

For synthetic experiments, we choose the tail direction using the sign ofe0e\_\{0\}\. Equivalently, defines=sign⁡\(e0\)s=\\operatorname\{sign\}\(e\_\{0\}\)and use the scalar score

h⁡\(Z\)=s​CEi​j​\(Z\)\.h\(Z\)=s\\,\\mathrm\{CE\}\_\{ij\}\(Z\)\.Then larger values ofh⁡\(Z\)h\(Z\)always correspond to more extreme effects in the selected direction\. A target score levelλ\\lambdadefines the event

𝒜λ=\{Z:h⁡\(Z\)≥λ\}\.\\mathcal\{A\}\_\{\\lambda\}=\\\{Z:h\(Z\)\\geq\\lambda\\\}\.Whene0\>0e\_\{0\}\>0, this is the right\-tail eventℰi​j\+​\(λ\)=\{Z:CEi​j​\(Z\)≥λ\}\\mathcal\{E\}\_\{ij\}^\{\+\}\(\\lambda\)=\\\{Z:\\mathrm\{CE\}\_\{ij\}\(Z\)\\geq\\lambda\\\}\. Whene0<0e\_\{0\}<0, it is the left\-tail eventℰi​j−​\(λ\)=\{Z:CEi​j​\(Z\)≤−λ\}\\mathcal\{E\}\_\{ij\}^\{\-\}\(\\lambda\)=\\\{Z:\\mathrm\{CE\}\_\{ij\}\(Z\)\\leq\-\\lambda\\\}\. This is the same score\-level representation used by the MLS sampler in Sec\.[4](https://arxiv.org/html/2608.12640#S4)\.

1. 1\.Initial score level\.We initialize the threshold grid at the magnitude of the ground\-truth effect,λ1=\|e0\|\\lambda\_\{1\}=\|e\_\{0\}\|\.
2. 2\.Pilot run\.We run a short, low\-budget adaptive multilevel\-splitting pilot to identify a more extreme target score levelλT\\lambda\_\{T\}such that the corresponding tail probability is already in the rare\-event regime, approximately10−610^\{\-6\}or smaller\.
3. 3\.Threshold grid\.We form a monotone increasing sequence ofTTtarget score levelsλ1<λ2<⋯<λT\\lambda\_\{1\}<\\lambda\_\{2\}<\\cdots<\\lambda\_\{T\}\. Largerλ\\lambdaalways corresponds to a more extreme one\-sided effect in the selected direction\.

### A\.3Detailed result for d=4 baseline experiment

Table[3](https://arxiv.org/html/2608.12640#A1.T3)lists statistics for thed=4d=4baseline experiments across seven methods\.

Table 3:d=4d=4per\-case probability estimates across runs \(mean±\\pmSD\)\. MLS refers to multilevel splitting, and SC refers to single\-chain sampling\.CaseEnumerationp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}PARNI\-MLSp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}Structure\-MCMC\-MLSp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}SPNp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}PARNI\-SCp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}Structure\-MCMC\-SCp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}DIBSp^±SD\\hat\{p\}\\pm\\mathrm\{SD\}11\.87×10−51\.87\\times 10^\{\-5\}±7\.25×10−6\\pm\\,7\.25\\times 10^\{\-6\}1\.83×10−51\.83\\times 10^\{\-5\}±1\.07×10−5\\pm\\,1\.07\\times 10^\{\-5\}1\.9×10−51\.9\\times 10^\{\-5\}±6\.09×10−6\\pm\\,6\.09\\times 10^\{\-6\}1\.42×10−51\.42\\times 10^\{\-5\}±6×10−6\\pm\\,6\\times 10^\{\-6\}1\.84×10−51\.84\\times 10^\{\-5\}±5\.32×10−6\\pm\\,5\.32\\times 10^\{\-6\}1\.4×10−51\.4\\times 10^\{\-5\}±1\.01×10−5\\pm\\,1\.01\\times 10^\{\-5\}2\.32×10−22\.32\\times 10^\{\-2\}±1\.35×10−3\\pm\\,1\.35\\times 10^\{\-3\}22\.3×10−42\.3\\times 10^\{\-4\}±1\.13×10−4\\pm\\,1\.13\\times 10^\{\-4\}2\.5×10−42\.5\\times 10^\{\-4\}±1\.07×10−4\\pm\\,1\.07\\times 10^\{\-4\}2\.41×10−42\.41\\times 10^\{\-4\}±8\.04×10−5\\pm\\,8\.04\\times 10^\{\-5\}1\.37×10−21\.37\\times 10^\{\-2\}±8\.62×10−4\\pm\\,8\.62\\times 10^\{\-4\}2\.28×10−42\.28\\times 10^\{\-4\}±7\.49×10−5\\pm\\,7\.49\\times 10^\{\-5\}2\.53×10−42\.53\\times 10^\{\-4\}±5\.78×10−5\\pm\\,5\.78\\times 10^\{\-5\}1\.35×10−21\.35\\times 10^\{\-2\}±7\.72×10−4\\pm\\,7\.72\\times 10^\{\-4\}31\.65×10−41\.65\\times 10^\{\-4\}±5\.96×10−5\\pm\\,5\.96\\times 10^\{\-5\}1\.7×10−41\.7\\times 10^\{\-4\}±5\.47×10−5\\pm\\,5\.47\\times 10^\{\-5\}1\.61×10−41\.61\\times 10^\{\-4\}±1\.41×10−5\\pm\\,1\.41\\times 10^\{\-5\}1\.83×10−11\.83\\times 10^\{\-1\}±7\.92×10−3\\pm\\,7\.92\\times 10^\{\-3\}1\.67×10−41\.67\\times 10^\{\-4\}±2\.42×10−5\\pm\\,2\.42\\times 10^\{\-5\}1\.79×10−41\.79\\times 10^\{\-4\}±2\.31×10−5\\pm\\,2\.31\\times 10^\{\-5\}1\.81×10−11\.81\\times 10^\{\-1\}±8\.38×10−3\\pm\\,8\.38\\times 10^\{\-3\}46\.17×10−56\.17\\times 10^\{\-5\}±4\.07×10−5\\pm\\,4\.07\\times 10^\{\-5\}5\.89×10−55\.89\\times 10^\{\-5\}±3\.67×10−5\\pm\\,3\.67\\times 10^\{\-5\}6\.79×10−56\.79\\times 10^\{\-5\}±2\.13×10−5\\pm\\,2\.13\\times 10^\{\-5\}2\.67×10−22\.67\\times 10^\{\-2\}±1\.76×10−3\\pm\\,1\.76\\times 10^\{\-3\}7\.06×10−57\.06\\times 10^\{\-5\}±2\.85×10−5\\pm\\,2\.85\\times 10^\{\-5\}7\.21×10−57\.21\\times 10^\{\-5\}±3\.06×10−5\\pm\\,3\.06\\times 10^\{\-5\}2\.65×10−22\.65\\times 10^\{\-2\}±1\.66×10−3\\pm\\,1\.66\\times 10^\{\-3\}51\.6×10−41\.6\\times 10^\{\-4\}±5\.56×10−5\\pm\\,5\.56\\times 10^\{\-5\}1\.71×10−41\.71\\times 10^\{\-4\}±5\.9×10−5\\pm\\,5\.9\\times 10^\{\-5\}1\.62×10−41\.62\\times 10^\{\-4\}±1\.74×10−5\\pm\\,1\.74\\times 10^\{\-5\}3\.34×10−23\.34\\times 10^\{\-2\}±1\.88×10−3\\pm\\,1\.88\\times 10^\{\-3\}1\.55×10−41\.55\\times 10^\{\-4\}±4\.54×10−5\\pm\\,4\.54\\times 10^\{\-5\}1\.69×10−41\.69\\times 10^\{\-4\}±4\.88×10−5\\pm\\,4\.88\\times 10^\{\-5\}3\.34×10−23\.34\\times 10^\{\-2\}±1\.63×10−3\\pm\\,1\.63\\times 10^\{\-3\}62\.15×10−32\.15\\times 10^\{\-3\}±7\.08×10−4\\pm\\,7\.08\\times 10^\{\-4\}2\.21×10−32\.21\\times 10^\{\-3\}±9\.14×10−4\\pm\\,9\.14\\times 10^\{\-4\}2\.16×10−32\.16\\times 10^\{\-3\}±7\.4×10−4\\pm\\,7\.4\\times 10^\{\-4\}2\.27×10−22\.27\\times 10^\{\-2\}±1\.61×10−3\\pm\\,1\.61\\times 10^\{\-3\}2\.09×10−32\.09\\times 10^\{\-3\}±1\.38×10−4\\pm\\,1\.38\\times 10^\{\-4\}2\.15×10−32\.15\\times 10^\{\-3\}±1\.78×10−4\\pm\\,1\.78\\times 10^\{\-4\}2\.3×10−22\.3\\times 10^\{\-2\}±1\.55×10−3\\pm\\,1\.55\\times 10^\{\-3\}72\.44×10−42\.44\\times 10^\{\-4\}±1\.13×10−4\\pm\\,1\.13\\times 10^\{\-4\}2\.4×10−42\.4\\times 10^\{\-4\}±1\.06×10−4\\pm\\,1\.06\\times 10^\{\-4\}2\.19×10−42\.19\\times 10^\{\-4\}±2\.48×10−5\\pm\\,2\.48\\times 10^\{\-5\}9\.71×10−29\.71\\times 10^\{\-2\}±6\.54×10−3\\pm\\,6\.54\\times 10^\{\-3\}2\.22×10−42\.22\\times 10^\{\-4\}±2\.25×10−5\\pm\\,2\.25\\times 10^\{\-5\}2\.42×10−42\.42\\times 10^\{\-4\}±2\.5×10−5\\pm\\,2\.5\\times 10^\{\-5\}9\.65×10−29\.65\\times 10^\{\-2\}±6\.06×10−3\\pm\\,6\.06\\times 10^\{\-3\}83\.09×10−53\.09\\times 10^\{\-5\}±2\.84×10−5\\pm\\,2\.84\\times 10^\{\-5\}4\.48×10−54\.48\\times 10^\{\-5\}±2\.71×10−5\\pm\\,2\.71\\times 10^\{\-5\}3\.1×10−53\.1\\times 10^\{\-5\}±6\.72×10−6\\pm\\,6\.72\\times 10^\{\-6\}5\.76×10−35\.76\\times 10^\{\-3\}±6\.65×10−4\\pm\\,6\.65\\times 10^\{\-4\}2\.82×10−52\.82\\times 10^\{\-5\}±8\.1×10−6\\pm\\,8\.1\\times 10^\{\-6\}3\.4×10−53\.4\\times 10^\{\-5\}±1\.35×10−5\\pm\\,1\.35\\times 10^\{\-5\}5\.86×10−35\.86\\times 10^\{\-3\}±5\.4×10−4\\pm\\,5\.4\\times 10^\{\-4\}93\.22×10−43\.22\\times 10^\{\-4\}±1\.34×10−4\\pm\\,1\.34\\times 10^\{\-4\}3\.25×10−43\.25\\times 10^\{\-4\}±1\.33×10−4\\pm\\,1\.33\\times 10^\{\-4\}3\.22×10−43\.22\\times 10^\{\-4\}±5\.19×10−5\\pm\\,5\.19\\times 10^\{\-5\}1\.32×10−21\.32\\times 10^\{\-2\}±1\.46×10−3\\pm\\,1\.46\\times 10^\{\-3\}3\.09×10−43\.09\\times 10^\{\-4\}±6\.17×10−5\\pm\\,6\.17\\times 10^\{\-5\}3\.3×10−43\.3\\times 10^\{\-4\}±4\.04×10−5\\pm\\,4\.04\\times 10^\{\-5\}1\.33×10−21\.33\\times 10^\{\-2\}±1\.45×10−3\\pm\\,1\.45\\times 10^\{\-3\}101\.33×10−51\.33\\times 10^\{\-5\}±1\.02×10−5\\pm\\,1\.02\\times 10^\{\-5\}1\.98×10−51\.98\\times 10^\{\-5\}±1\.78×10−5\\pm\\,1\.78\\times 10^\{\-5\}1\.37×10−51\.37\\times 10^\{\-5\}±1\.11×10−5\\pm\\,1\.11\\times 10^\{\-5\}1\.82×10−21\.82\\times 10^\{\-2\}±1\.55×10−3\\pm\\,1\.55\\times 10^\{\-3\}1\.68×10−51\.68\\times 10^\{\-5\}±1\.52×10−5\\pm\\,1\.52\\times 10^\{\-5\}1\.21×10−51\.21\\times 10^\{\-5\}±1\.58×10−5\\pm\\,1\.58\\times 10^\{\-5\}1\.8×10−21\.8\\times 10^\{\-2\}±1\.72×10−3\\pm\\,1\.72\\times 10^\{\-3\}
### A\.4Time Complexity

We separate the cost of the outer adaptive multilevel\-splitting loop from the cost of the inner MCMC mutation kernel\. LetNNbe the number of particles,mmthe number of MCMC mutation steps per particle per level,KKthe number of splitting levels, andCMHC\_\{\\mathrm\{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 mutatesNNparticles formmsteps, computesNNscores, sorts or partially sorts the scores to choose the next empirical quantile, and resamples the survivors\. The overall cost is therefore

O⁡\(K​N​m​CMH\+K​N​log⁡N\),O\\\!\\left\(KNmC\_\{\\mathrm\{MH\}\}\+KN\\log N\\right\),where theK​N​log⁡NKN\\log Nterm comes from quantile selection\. If selection is implemented by a linear\-time order\-statistic routine, this sorting term can be reduced toO⁡\(K​N\)O\(KN\)and the mutation cost dominates\.

The average proposal cost depends on the mixture of structure and weight moves\. Ifpstructp\_\{\\mathrm\{struct\}\}is the probability of proposing a structure move, then

CMH≈pstruct​Cstruct\+\(1−pstruct\)​Cweight\+Ch,C\_\{\\mathrm\{MH\}\}\\approx p\_\{\\mathrm\{struct\}\}C\_\{\\mathrm\{struct\}\}\+\(1\-p\_\{\\mathrm\{struct\}\}\)C\_\{\\mathrm\{weight\}\}\+C\_\{h\},whereCstructC\_\{\\mathrm\{struct\}\}is the cost of proposing and scoring a graph update,CweightC\_\{\\mathrm\{weight\}\}is the cost of a blocked coefficient refresh, andChC\_\{h\}is the cost of evaluating the scalar scoreh⁡\(Z\)h\(Z\)\. In the linear\-Gaussian implementation, the coefficient refreshes use node\-wise conjugate posterior updates, whileh⁡\(Z\)h\(Z\)is computed from the linear total\-effect matrix\. A direct matrix inverse givesCh=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 probabilityp⋆p\_\{\\star\}\. A direct posterior sampler needs about1/p⋆1/p\_\{\\star\}samples to see one sample from an event of posterior probabilityp⋆p\_\{\\star\}, andO⁡\(1/\(r2​p⋆\)\)O\(1/\(r^\{2\}p\_\{\\star\}\)\)samples to estimate that probability with fixed relative errorrr\. 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ρ\\rhoat each level, thenKKgrows roughly likelog⁡\(p⋆\)/log⁡\(ρ\)\\log\(p\_\{\\star\}\)/\\log\(\\rho\), so the leading cost grows approximately logarithmically in1/p⋆1/p\_\{\\star\}rather than linearly in1/p⋆1/p\_\{\\star\}\. This is why multilevel splitting remains useful in the deep\-tail regimes where unconditional sampling degenerates\.

##### Runtime measurements\.

Table[4](https://arxiv.org/html/2608.12640#A1.T4)reports wall\-clock running times for thed=4d=4case 1 experiment\. Table[5](https://arxiv.org/html/2608.12640#A1.T5)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 thed=4d=4case 1 experiment over 10 runs\. We report mean±\\pmstandard deviation in seconds\.MethodTime \(s\)MLS\-PARNI\-DAG192\.3±4\.8192\.3\\pm 4\.8MLS\-Structure\-MCMC56\.1±2\.156\.1\\pm 2\.1PARNI\-DAG268\.2±5\.3268\.2\\pm 5\.3Structure\-MCMC105\.5±3\.0105\.5\\pm 3\.0Enumeration\-Baseline301\.6±4\.6301\.6\\pm 4\.6OrderSPN166\.8±3\.9166\.8\\pm 3\.9DiBS434\.2±6\.7434\.2\\pm 6\.7Table 5:Average running time per outer MLS loop in seconds\. Standard deviations over 10 runs are shown in parentheses\.Kerneld=4d=4d=8d=8d=16d=16d=32d=32PARNI\-DAG38\.4 \(1\.1\)56\.8 \(1\.6\)92\.4 \(2\.8\)172\.5 \(4\.9\)Structure\-MCMC11\.2 \(0\.5\)19\.5 \(0\.8\)34\.8 \(1\.2\)58\.5 \(2\.1\)

### A\.5Extension 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\)\\pi\(Z\), a scalar scoreh⁡\(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\)≥λ\}\\pi\(Z\)\\mathbf\{1\}\\\{h\(Z\)\\geq\\lambda\\\}\. 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\}\\\!\\left\(X\_\{\\mathrm\{pa\}\_\{G\}\(j\)\};\\phi\_\{j\}\\right\)\+\\epsilon\_\{j\},\\qquad\\epsilon\_\{j\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{j\}^\{2\}\),and replace the graph–weight stateZ=\(G,B\)Z=\(G,B\)by a graph–mechanism stateZ=\(G,ϕ,σ\)Z=\(G,\\phi,\\sigma\)\. The posterior target would becomep⁡\(G,ϕ,σ∣𝒟\)p\(G,\\phi,\\sigma\\mid\\mathcal\{D\}\), 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\)=\\mathbb\{E\}\_\{Z\}\\\!\\left\[X\_\{j\}\\mid\\mathrm\{do\}\(X\_\{i\}=x\+\\Delta\)\\right\]\-\\mathbb\{E\}\_\{Z\}\\\!\\left\[X\_\{j\}\\mid\\mathrm\{do\}\(X\_\{i\}=x\)\\right\],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\.6Sachs experiment calculations and tables

This section provides supplementary numerical details for the Sachs case study in Sec\.[6](https://arxiv.org/html/2608.12640#S6), including \(i\) the qualitative reference network, \(ii\) estimated posterior masses of the tail events, and \(iii\) pathway\-level decompositions of conditional effects\.

##### Reference network\.

Figure[4](https://arxiv.org/html/2608.12640#A1.F4)shows the 11\-node, 20\-edge consensus signaling network used as a qualitative reference in Sec\.[6](https://arxiv.org/html/2608.12640#S6)\.

![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/sachs_graph.png)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](https://arxiv.org/html/2608.12640#A1.T6)reports \(i\) the*geometric mean*estimate ofPr⁡\(ℰ∣𝒟\)\\Pr\(\\mathcal\{E\}\\mid\\mathcal\{D\}\)\(i\.e\.,exp⁡\(𝔼⁡\[log⁡p^\]\)\\exp\(\\mathbb\{E\}\[\\log\\hat\{p\}\]\)across runs\) and \(ii\) the mean±\\pmsd of−log⁡p^\-\\log\\hat\{p\}across runs\.

Table 6:Estimated posterior mass of Sachs tail events over 10 MLS runs\.ConditionPr⁡\(ℰ∣𝒟\)\\Pr\(\\mathcal\{E\}\\mid\\mathcal\{D\}\)\(geo\-mean\)−log⁡p^\-\\log\\hat\{p\}\(mean±\\pmsd\)Cond\-PIP2\.2×10−72\.2\\times 10^\{\-7\}15\.32±0\.4415\.32\\pm 0\.44Cond\-ERK1\.3×10−71\.3\\times 10^\{\-7\}15\.84±0\.5315\.84\\pm 0\.53Cond\-Joint5\.2×10−75\.2\\times 10^\{\-7\}14\.47±0\.6714\.47\\pm 0\.67
##### Path decomposition under Cond\-PIP\.

Across samples drawn from the conditional posterior

p⁡\(Z∣𝒟,ℰPIP​\(tPIP\)\),p\\\!\\left\(Z\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t\_\{\\mathrm\{PIP\}\}\)\\right\),the two canonical pathwaysPIP3→PIP2\\mathrm\{PIP3\}\\\!\\to\\\!\\mathrm\{PIP2\}andPIP3→Plcg→PIP2\\mathrm\{PIP3\}\\\!\\to\\\!\\mathrm\{Plcg\}\\\!\\to\\\!\\mathrm\{PIP2\}are present in*100%*of samples\. These two routes account for essentially the entire conditional mean effect:

𝔼\[CEPIP3,PIP2\(Z\)∣𝒟,ℰPIP\(tPIP\)\]\\displaystyle\\mathbb\{E\}\\\!\\left\[\\mathrm\{CE\}\_\{\\mathrm\{PIP3\},\\mathrm\{PIP2\}\}\(Z\)\\mid\\mathcal\{D\},\\mathcal\{E\}\_\{\\mathrm\{PIP\}\}\(t\_\{\\mathrm\{PIP\}\}\)\\right\]\(14\)≈0\.6142⏟direct\+0\.1614⏟via Plcg\+2⋅10−4⏟other paths≈0\.7758\.\\displaystyle\\approx\\underbrace\{0\.6142\}\_\{\\text\{direct\}\}\+\\underbrace\{0\.1614\}\_\{\\text\{via Plcg\}\}\+\\underbrace\{2\\cdot 10^\{\-4\}\}\_\{\\text\{other paths\}\}\\approx 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](https://arxiv.org/html/2608.12640#A1.T7)reports posterior effect statistics \(e\.g\.,Pr⁡\(CE\>0\)\\Pr\(\\mathrm\{CE\}\>0\), conditional means, and path counts\) across the four regimes\. Table[8](https://arxiv.org/html/2608.12640#A1.T8)further decomposes the conditional mean effects into direct and mediated contributions for the two target pairs \(PIP3→\\toPIP2 and Erk→\\toAkt\), together with the frequency of the canonical mediated paths\.

![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/Sachs_uncondition_graph.png)\(a\)Unconditioned posterior
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/Sachs_parni_case1_graph.png)\(b\)Cond\-PIP:CEPIP3→PIP2\>tPIP\\mathrm\{CE\}\_\{\\mathrm\{PIP3\}\\to\\mathrm\{PIP2\}\}\>t\_\{\\mathrm\{PIP\}\}
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/Sachs_parni_case2_graph.png)\(c\)Cond\-ERK:CEErk→Akt\>tERK\\mathrm\{CE\}\_\{\\mathrm\{Erk\}\\to\\mathrm\{Akt\}\}\>t\_\{\\mathrm\{ERK\}\}
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/Sachs_parni_multi_graph.png)\(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\) andmean edge weight\(conditional mean ofBu​vB\_\{uv\}given the edge is present; right\)\.Table 7:Comparison across four conditioning settings for two target pairs\.Pr⁡\(CE\>0\)\\Pr\(\\mathrm\{CE\}\>0\)is the posterior mass with a directed influence \(i\.e\., at least one directed path\),𝔼⁡\[CE∣CE\>0\]\\mathbb\{E\}\[\\mathrm\{CE\}\\mid\\mathrm\{CE\}\>0\]and quantiles are computed on the nonzero subset, and \#paths is the median number of directed paths among samples withCE\>0\\mathrm\{CE\}\>0\. We run PARNI\-DAG single chain with200,000200,000iterations and collect180000180000samples for the unconditioned case\. For each of the other conditions, we collect500500survivor samples from 10 runs with5050samples for each run\.PIP3→\\toPIP2Erk→\\toAktConditionPr⁡\(CE\>0\)\\Pr\(\\mathrm\{CE\}\>0\)𝔼⁡\[CE∣CE\>0\]\\mathbb\{E\}\[\\mathrm\{CE\}\\mid\\mathrm\{CE\}\>0\]\#pathsmed\{\}\_\{\\text\{med\}\}Pr⁡\(CE\>0\)\\Pr\(\\mathrm\{CE\}\>0\)𝔼⁡\[CE∣CE\>0\]\\mathbb\{E\}\[\\mathrm\{CE\}\\mid\\mathrm\{CE\}\>0\]\#pathsmed\{\}\_\{\\text\{med\}\}Unconditioned0\.1090\.600\[0\.543, 0\.681\]40\.0850\.569\[0\.532, 0\.617\]1Cond\-PIP\(CEPIP\>tPIP\\mathrm\{CE\}\_\{\\mathrm\{PIP\}\}\>t\_\{\\mathrm\{PIP\}\}\)1\.0000\.776\[0\.771, 0\.782\]20\.4980\.596\[0\.538, 0\.634\]1Cond\-ERK\(CEERK\>tERK\\mathrm\{CE\}\_\{\\mathrm\{ERK\}\}\>t\_\{\\mathrm\{ERK\}\}\)0\.9840\.655\[0\.562, 0\.710\]101\.0000\.673\[0\.670, 0\.676\]1Cond\-Joint\(CEPIP\>tPIP′∧CEERK\>tERK′\\mathrm\{CE\}\_\{\\mathrm\{PIP\}\}\>t^\{\\prime\}\_\{\\mathrm\{PIP\}\}\\wedge\\mathrm\{CE\}\_\{\\mathrm\{ERK\}\}\>t^\{\\prime\}\_\{\\mathrm\{ERK\}\}\)1\.0000\.747\[0\.742, 0\.758\]31\.0000\.656\[0\.652, 0\.662\]1Table 8:Mechanism decomposition across conditions, using the two pathways most relevant for interpretation: for PIP3→\\toPIP2 we use the direct path \(PIP3→\\toPIP2\) and the Plcg\-mediated path \(PIP3→\\toPlcg→\\toPIP2\); for Erk→\\toAkt we use the direct path \(Erk→\\toAkt\) and the Plcg\-mediated path \(Erk→\\toPlcg→\\toAkt\)\. For unconditioned samples, we select those withCE\>0\\mathrm\{CE\}\>0, and calculate𝔼⁡\[CE∣CE\>0\]\\mathbb\{E\}\[\\mathrm\{CE\}\\mid\\mathrm\{CE\}\>0\]\. The last column reports the frequency that the canonical mediated path edges are present \(over all samples under the condition\)\.PairCondition𝔼⁡\[CE∣CE\>0\]\\mathbb\{E\}\[\\mathrm\{CE\}\\mid\\mathrm\{CE\}\>0\]Direct contribVia\-Plcg contribPr⁡\(Via\-Plcg path\)\\Pr\(\\text\{Via\-Plcg path\}\)PIP3→\\toPIP2Unconditioned0\.6000\.560 \(93\.4%\)0\.076 \(12\.7%\)0\.064Cond\-PIP0\.7760\.614 \(79\.2%\)0\.161 \(20\.8%\)1\.000Cond\-ERK0\.6550\.559 \(85\.2%\)0\.140 \(21\.4%\)0\.975Cond\-Joint0\.7470\.598 \(79\.9%\)0\.154 \(20\.5%\)1\.000Erk→\\toAktUnconditioned0\.5690\.563 \(98\.9%\)0\.000 \(0\.0%\)0\.0004Cond\-PIP0\.5960\.606 \(101\.6%\)0\.001 \(0\.2%\)0\.026Cond\-ERK0\.6730\.668 \(99\.3%\)0\.004 \(0\.6%\)0\.177Cond\-Joint0\.6560\.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\.7Synthetic multi\-effect heat maps

Figures[6](https://arxiv.org/html/2608.12640#A1.F6)and[7](https://arxiv.org/html/2608.12640#A1.F7)provide the detailed heat\-map visualizations for the multi\-effect conditioning experiments discussed in Sec\.[5\.2\.2](https://arxiv.org/html/2608.12640#S5.SS2.SSS2)\. 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\.

![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=4_Ground_Truth_Graph.png)\(a\)Ground\-truth graph\.
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=4_PARNI-DAG_Samples_strong_constraints.png)\(b\)Conditional posterior under strong multi\-effect constraints\.
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=4_PARNI-DAG_Samples_under_weak_constraints.png)\(c\)Conditional posterior under weak multi\-effect constraints\.

Figure 6:d=4d=4multi\-effect conditioning results\. Each panel summarizes edge frequencies and average edge weights\. Stronger constraints concentrate posterior mass more tightly around the ground\-truth mechanism\.![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=8_Ground_Truth_Graph.png)\(a\)Ground\-truth graph\.
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=8_PARNI-DAG_Samples_under_strong_constraints.png)\(b\)Conditional posterior under strong multi\-effect constraints\.
![Refer to caption](https://arxiv.org/html/2608.12640v1/figures/d=8_PARNI-DAG_Samples_under_weak_constraints.png)\(c\)Conditional posterior under weak multi\-effect constraints\.

Figure 7:d=8d=8multi\-effect conditioning results, with the same layout as Fig\.[6](https://arxiv.org/html/2608.12640#A1.F6)\. The strong constraints produce sharper edge and weight summaries than the weak constraints\.
### A\.8Pseudo\-code

Algorithm[2](https://arxiv.org/html/2608.12640#alg2)summarizes the adaptive multilevel\-splitting loop, and Algorithm[3](https://arxiv.org/html/2608.12640#alg3)summarizes the level\-truncated MCMC mutation kernel\. The notation matches Sec\.[4](https://arxiv.org/html/2608.12640#S4):hhis the scalar score,λ⋆\\lambda\_\{\\star\}is the target level, and𝒜λ=\{Z:h⁡\(Z\)≥λ\}\\mathcal\{A\}\_\{\\lambda\}=\\\{Z:h\(Z\)\\geq\\lambda\\\}is the corresponding score\-level event\.

Algorithm 2Adaptive Multilevel Splitting for Conditional Causal Discovery1:Data

𝒟\\mathcal\{D\}; posterior density

π⁡\(Z\)=p⁡\(G,B∣𝒟\)\\pi\(Z\)=p\(G,B\\mid\\mathcal\{D\}\); score

hh; target level

λ⋆\\lambda\_\{\\star\}; particle size

NN; survival fraction

ρ\\rho; mutation steps

mmper level; maximum number of levels

KmaxK\_\{\\max\}; structure kernel

𝖪𝖾𝗋𝗇𝖾𝗅\\mathsf\{Kernel\}\.

2:Tail\-probability estimate

p^⋆\\widehat\{p\}\_\{\\star\}and approximately constrained particles

𝒫\\mathcal\{P\}\.

3:Initialize particles

𝒫=\{Zn=\(Gn,Bn\)\}n=1N\\mathcal\{P\}=\\\{Z\_\{n\}=\(G\_\{n\},B\_\{n\}\)\\\}\_\{n=1\}^\{N\}by drawing

Gn∼p⁡\(G\)G\_\{n\}\\sim p\(G\)and

Bn∼p⁡\(B∣Gn,𝒟\)B\_\{n\}\\sim p\(B\\mid G\_\{n\},\\mathcal\{D\}\)\.

4:Set

λ0←−∞\\lambda\_\{0\}\\leftarrow\-\\infty,

p^⋆←1\\widehat\{p\}\_\{\\star\}\\leftarrow 1, and

k←0k\\leftarrow 0\.

5:

𝒫←MCMCMutation​\(𝒫,𝒟,h,λ0,m,𝖪𝖾𝗋𝗇𝖾𝗅\)\\mathcal\{P\}\\leftarrow\\textsc\{MCMCMutation\}\(\\mathcal\{P\},\\mathcal\{D\},h,\\lambda\_\{0\},m,\\mathsf\{Kernel\}\)\.

6:while

k<Kmaxk<K\_\{\\max\}do

7:Compute scores

sn←h⁡\(Zn\)s\_\{n\}\\leftarrow h\(Z\_\{n\}\)for all

Zn∈𝒫Z\_\{n\}\\in\\mathcal\{P\}\.

8:

λ~k\+1←Q1−ρ​\(\{sn\}n=1N\)\\widetilde\{\\lambda\}\_\{k\+1\}\\leftarrow Q\_\{1\-\\rho\}\(\\\{s\_\{n\}\\\}\_\{n=1\}^\{N\}\)\.

9:

λk\+1←min⁡\{λ~k\+1,λ⋆\}\\lambda\_\{k\+1\}\\leftarrow\\min\\\{\\widetilde\{\\lambda\}\_\{k\+1\},\\lambda\_\{\\star\}\\\}\.

10:

𝒮←\{Zn∈𝒫:h⁡\(Zn\)≥λk\+1\}\\mathcal\{S\}\\leftarrow\\\{Z\_\{n\}\\in\\mathcal\{P\}:h\(Z\_\{n\}\)\\geq\\lambda\_\{k\+1\}\\\}\.

11:

β^k\+1←\|𝒮\|/N\\widehat\{\\beta\}\_\{k\+1\}\\leftarrow\|\\mathcal\{S\}\|/N\.

12:

p^⋆←p^⋆​β^k\+1\\widehat\{p\}\_\{\\star\}\\leftarrow\\widehat\{p\}\_\{\\star\}\\,\\widehat\{\\beta\}\_\{k\+1\}\.

13:if

\|𝒮\|=0\|\\mathcal\{S\}\|=0then

14:return

\(0,∅\)\(0,\\emptyset\)\.

15:endif

16:Resample

NNparticles from

𝒮\\mathcal\{S\}with replacement to form

𝒫\\mathcal\{P\}\.

17:

𝒫←MCMCMutation​\(𝒫,𝒟,h,λk\+1,m,𝖪𝖾𝗋𝗇𝖾𝗅\)\\mathcal\{P\}\\leftarrow\\textsc\{MCMCMutation\}\(\\mathcal\{P\},\\mathcal\{D\},h,\\lambda\_\{k\+1\},m,\\mathsf\{Kernel\}\)\.

18:if

λk\+1=λ⋆\\lambda\_\{k\+1\}=\\lambda\_\{\\star\}then

19:return

\(p^⋆,𝒫\)\(\\widehat\{p\}\_\{\\star\},\\mathcal\{P\}\)\.

20:endif

21:

k←k\+1k\\leftarrow k\+1\.

22:endwhile

23:return

\(p^⋆,𝒫\)\(\\widehat\{p\}\_\{\\star\},\\mathcal\{P\}\)\.

Algorithm 3Level\-Truncated MCMC Mutation over Graph–Weight States1:Population

𝒫\\mathcal\{P\}; data

𝒟\\mathcal\{D\}; score

hh; level

λ\\lambda; mutation steps

mm; structure kernel

𝖪𝖾𝗋𝗇𝖾𝗅\\mathsf\{Kernel\}; structure\-move probability

pstructp\_\{\\mathrm\{struct\}\}\.

2:Mutated population

𝒫′\\mathcal\{P\}^\{\\prime\}\.

3:

𝒫′←∅\\mathcal\{P\}^\{\\prime\}\\leftarrow\\emptyset\.

4:for all

Z=\(G,B\)∈𝒫Z=\(G,B\)\\in\\mathcal\{P\}do

5:for

r=1r=1to

mmdo

6:Draw

u∼Uniform⁡\(0,1\)u\\sim\\mathrm\{Uniform\}\(0,1\)\.

7:if

u<pstructu<p\_\{\\mathrm\{struct\}\}then

8:Propose a graph

G′∼qG\(⋅∣G\)G^\{\\prime\}\\sim q\_\{G\}\(\\cdot\\mid G\)using

𝖪𝖾𝗋𝗇𝖾𝗅\\mathsf\{Kernel\}\.

9:Refresh affected coefficient blocks from

p⁡\(B∣G′,𝒟\)p\(B\\mid G^\{\\prime\},\\mathcal\{D\}\)and copy unchanged blocks to obtain

B′B^\{\\prime\}\.

10:else

11:Set

G′←GG^\{\\prime\}\\leftarrow Gand refresh one coefficient block from

p⁡\(B∣G,𝒟\)p\(B\\mid G,\\mathcal\{D\}\)to obtain

B′B^\{\\prime\}\.

12:endif

13:Set

Z′←\(G′,B′\)Z^\{\\prime\}\\leftarrow\(G^\{\\prime\},B^\{\\prime\}\)and let

q⁡\(Z′∣Z\)q\(Z^\{\\prime\}\\mid Z\)denote the complete proposal density\.

14:if

h⁡\(Z′\)<λh\(Z^\{\\prime\}\)<\\lambdathen

15:Reject

Z′Z^\{\\prime\}and continue\.

16:else

17:

α←min⁡\{1,π⁡\(Z′\)​q​\(Z∣Z′\)π⁡\(Z\)​q​\(Z′∣Z\)\}\\alpha\\leftarrow\\min\\\!\\left\\\{1,\\dfrac\{\\pi\(Z^\{\\prime\}\)q\(Z\\mid Z^\{\\prime\}\)\}\{\\pi\(Z\)q\(Z^\{\\prime\}\\mid Z\)\}\\right\\\}\.

18:Draw

a∼Uniform⁡\(0,1\)a\\sim\\mathrm\{Uniform\}\(0,1\)\.

19:if

a≤αa\\leq\\alphathen

20:

Z←Z′Z\\leftarrow Z^\{\\prime\}\.

21:endif

22:endif

23:endfor

24:

𝒫′←𝒫′∪\{Z\}\\mathcal\{P\}^\{\\prime\}\\leftarrow\\mathcal\{P\}^\{\\prime\}\\cup\\\{Z\\\}\.

25:endfor

26:return

𝒫′\\mathcal\{P\}^\{\\prime\}\.