Exogenous Matching: Learning Good Proposals for Tractable Counterfactual Estimation
Authors: Yikang Chen
Dehui du
Lili Tian
Venue: NeurIPS 2024 poster (2024)
Keywords: causality
causal inference
counterfactual estimation
importance sampling
normalizing flows
Primary area: causal_inference

Abstract: We propose an importance sampling method for tractable and efficient estimation of counterfactual expressions in general settings, named Exogenous Matching. By minimizing a common upper bound of counterfactual estimators, we transform the variance minimization problem into a conditional distribution learning problem, enabling its integration with existing conditional distribution modeling approaches. We validate the theoretical results through experiments under various types and settings of Structural Causal Models (SCMs) and demonstrate the outperformance on counterfactual estimation tasks compared to other existing importance sampling methods. We also explore the impact of injecting structural prior knowledge (counterfactual Markov boundaries) on the results. Finally, we apply this method to identifiable proxy SCMs and demonstrate the unbiasedness of the estimates, empirically illustrating the applicability of the method to practical scenarios.

Exogenous Matching: Learning Good Proposals for Tractable Counterfactual Estimation
Yikang Chen', Dehui Dul Lili Tian'
"Shanghai Key Laboratory of Trustworthy Computing, East China Normal University
Abstract
We propose importance sampling method for tractable and efficient estimation of counterfactual expressions in general settings, named Exogenous Matching. By minimizing a common upper bound of counterfactual estimators, we transform the variance minimization problem into a conditional distribution learning problem, enabling its integration with existing conditional distribution modeling approaches. We validate the theoretical results through experiments under various types and settings of Structural Causal Models (SCMs) and demonstrate the outperformance counterfactual estimation tasks compared to other existing importance sampling methods. We also explore the impact of injecting structural prior knowledge (counterfactual Markov boundaries) on the results. Finally, we apply this method identifiable proxy SCMs and demonstrate the unbiasedness of the estimates, empirically illustrating the applicability of the method to practical scenarios
1 Introduction
Counterfactual reasoning, considered one of the advanced cognitive abilities of humans, aims to address questions about hypothetical worlds that have not been observed. Answering these questions, typically in the form of "what-if" or "why." is crucial for attribution and explanation, thus counterfactual reasoning is widely applied in generating explanations [48]
Variables: o. Observatons: H, Intervestions:
submodel SCM under intervention Ma
SCM M
Observational Interventinal Counterfactual
Behavior
o were under A would he under
Natural Language
Cis
would be under K?
47
94 63), mak-
58 and evaluating class of gen- erative models called Structural Causal Models (SCMs) provides a semantic characterization of counterfactuals, consisting of of endoge- nous mechanisms and an exogenous distribu- tion, from which a collection of distributions described by specific formal languages is in- duced. These languages are organized into a progressively refined hierarchical structure, and each is associated with human activities: language Li is about seeing (observational), L2 is about doing (interventional), and L3 is about imagining (counterfactual). This hierarchy referred as the causal ladder (77) or the Pearl Causal Hierarchy (PCH) [5). Fig. provides a brief illustration of the relevant concepts. As the highest level of the PCH, counterfactuals entail the most subtle information
ing decisions [71]
37
95
Formal Language
P
P(o
P(a o.=0)
fairness [56] 113,
105
78
117].
PCH
La
La
SCM
M
M
M
Figure 1: A brief illustration of SCM, PCH and counterfactual concepts.
"Corresponding author (dhdu@sei.ecmueeduuu
'Code available http://///gsstuu.ooo/////////
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

in the hierarchy, making the evaluation of counterfactual languages an appealing objective that has received widespread attention in recent years.
Counterfactual reasoning consists of two subtasks: counterfactual identification and counterfactual estimation. For an expression Q in C3, identification is to determine the uniqueness of the answer to Q [86] 17 or its bounds [T4 25 66 under some assumptions, while estimation is to compute the specific value of 2. A series of literature 86 17 transforms Q into combination of expressions in Li and L2 through identification. However, estimating the transformed expression may be intractable. Another series of literature focuses on identifying and estimating Q specific settings, such as back
door [82 38 115 102], front door [43 109 98]. instrumental variables [14] 100 35], representation
decoupling [1 97 61], categorical exogenous [22]. exponential family 85 domain counterfactuals
[16]. There are also some literature applicable to relaxed cases, such as 22 45 46 57]. but still
less developed.
In recent years, a class of methods has emerged that simulate the structure of the original SCM using neural network architectures and learn causal mechanisms based on observed or intervened distributions, which we refer to as neural proxy SCMs. These methods include the use of VAEs
1T0
55
92 2 GANs [52 T08, normalizing flows [74 50 69 42
DDPMs
Here, 84
we focus only on models that provide identifiability results. Using bijections, 669
demonstrates the identifiability of learned causal mechanisms. However, constrained by bijections, they only provide tractable counterfactual estimation methods for fully observed endogenous variables. The closest method to performing tractable counterfactual estimation in general cases is NCM [107] [08]. which employs an optimization algorithm to concurrently achieve tractable counterfactual estimation and identification (or partial identifiability (93).. However, the counterfactual estimation depends on rejection sampling, thereby lacking scalability.
42
is commonly recognized that performing tractable counterfactual estimation in general cases is challenging, even when the SCM is fully specified. From the perspective of definition [5]. one needs to integrate over the entire exogenous space, which involves determining whether the constraints of different hypothetical worlds are simultaneously satisfied. From the perspective of complexity, although (34] has demonstrated that counterfactual reasoning is tractable if observational and interventional reasoning are tractable, it relies on algorithms with exponential complexity. Furthermore, according to (12), performing marginal inference on SCMs, including parameterized SCMs such as NCM, has been shown to be NP-hard.
This paper presents tractable and efficient counterfactual estimation method in general settings based on importance sampling. Here, by "general settings", we mean that the exogenous variables of the SCM can be discrete continuous, the expressions be estimated can involve an arbitrary finite number of hypothetical worlds, and the observations and interventions can also be arbitrary.
Specifically, the contributions of this paper
follows:
We propose a tractable and efficient importance sampling method for counterfactual estimation in general settings, which is based on optimizing upper bound on the variance of the estimators, and is formulated as a conditional distribution learning problem.
2. We explore the injection of structural prior knowledge (counterfactual Markov boundaries) into neural networks used for parameter optimization.
3. The experiments conducted on different SCM settings demonstrate that the proposed method outperforms existing importance sampling methods. Ablation studies empirically highlight the effectiveness of injecting structural prior knowledge. The experiments counterfactual estimation tasks conducted on two identifiable proxy SCMs illustrate the fcasibility of applying the method in conjunction with proxy SCMs to real-world problems.
2 Preliminaries
In this section, we provide the necessary background and definitions for our work, which are consistent with 5 g 114]. To maintain consistency in notation, we the uppercase letter X denote random variables and the lowercase letter x to denote their values. Bold uppercase letter X represents a set
of variables and x the corresponding set of values. The domain of the variable X is denoted as Ωx, and Qx for X {Xi...., X,,, } represents the product of domains X" Ωx.- We use Px to denote distribution over set of random variables X, and P(X X) (abbreviated as P(Y)) to denote

---

probability when X take values X' S Ωx. The lowercase p(x) represents the probability density function of Px if Q2x is continuous.
Structural causal models An SCM is a 4-uple M = (U.V.F where U is a set of exogenous variables; V is set of endogenous variables; JF, which describes the causal mechanisms of M, is a collection of functions (fv,, fv...... .v.......v.. } such that each endogenous variable V; is determined by function fv, Each fv. E F a mapping from (the domains of) Uv, Pav, to V,, where Uv, S U and Pay. S V VG The entire set F forms a mapping from U to V, while the uncertainty comes from probability distribution Pu over exogenous variables U.
An SCM yields a causal graph G. where each Vi E V is vertex, and the edges come in two types: there is directed edge (V; V2) between each pair V and V, E Pay;; and there is bidirected edge (V; ++ Vj) between each pair V; and VG, if Uv, U A U. This resulting graph is also known as directed mixed graph (DMG). Here, we assume that the SCM recursive, which means that there exists an order in F such that for any pair fi F, if f; < f.. then Vj & Pay. In particular, this implies that given values for exogenous variables (also known as a unit or an individual) u E Du all values of endogenous variables in V are also fixed (i.e., there exists a unique solution for endogenous variables w.r.t. u), and the resulting causal graph an acyclic directed mixed graph (ADMG). We denote Y (u) as the solution for Y C V given u.
In the paradigm of SCM, (perfect) intervention is described as replacing the mechanisms of some variables X with constants x. The model after this intervention is termed a submodel, defined as
Mr V. Pu), where Fx (fi Vi & X] U + x).
(1)
In submodel Mx. the solution of Y C U given u, denoted as Y,,, (u), is termed as the potential response, and is typically abbreviated as (u()).
Counterfactual events and probabilities The random variable Y E V corresponding to the submodel MM is denoted as Yx. termed counterfactual variable, with its domain denoted fy. A set of counterfactual variables from the same submodel MMx is denoted as Yx. and the set of counterfactual variables from & different submodels MM is denoted as Y. = = (((x))+
with the corresponding domain Ωy = X DY Let o-algebra Σy. C 20v,, where 20v, is the power set of fy,, The counterfactual probability space is then a triplet (Ωy NY., Py. where Py. is probability measure PY, EYx - [0, 1], and for any measurable set V. E VY.,, its value equals the Lebesgue integral,
/ (u)
(x.()))
(2)
no
where [u) is an indicator function, equal to if u E Ω and 0 otherwise. The set Cu(V)) =
{u (u) vV,, where Y. (u) denotes the potential responses U Yixx,) (u) of counterfactual variables in Y. w.r.t. u. We refer to the random event Y. E V. as a counterfactual event, and its probability P().)= Px. (D%.) is termed the counterfactual probability.
Counterfactual probability is the cornerstone of the language L3 in PCH. Specifically, all expressions in La consist of inequalities between polynomials (and their Boolean combinations) over terms of the form P()).). ). Here, we assume that the La expressions to be estimated can be represented by a finite number of P()..) terms (including those approximated by Monte Carlo methods). Therefore, the estimation of P()+) will become a focal point of the subsequent discussion.
Normalizing flows The normalizing flow To RRd Rd is transformation parameterized by 0 that maps observable samples x with d features to latent samples z distributed according to a simple distribution Pz with probability density function p(z). The mapping To is a diffeomorphism, which allows us to compute the density p(x) through the change-of variables
vaangell--vriablls formula:
log p(x) = log (x)) + log |det( V TT())))
When the Jacobian determinant of (x() is tractable to compute, we can use samples from the observable distribution and maximum likelihood estimation to train the flow Tg. Additionally, we can easily generate i.idd. samples conforming to the learned observable distribution by first sampling ( from the base distribution Pz. and then evaluating the inverse transformation T (2()))
3

---

3 Exogenous Matching
In this section, we present only the final theoretical conclusions. Detailed derivations of all theories and equations can be found in App. A
Assumptions In addition assuming that the SCM (U, (U, V, Pv) is recursive, the proposed method also assumes: for any X C V.x E 9xx u Ir. Fx(()) computable; ii) Ptr is sampleable and computable. These assumptions imply the feasibility of counterfactual generation, meaning we can draw arbitrary counterfactual samples, with Fx and Po both treated as black boxes. For example, in the context of neural proxy SCMs, Pu. the latent distribution of the generative model, is typically modeled as a simple distribution, while Fx is a neural network.
Importance sampling for Monte Carlo integration A widely used method for estimating the Lebesgue integral is Monte Carlo integration. However, the indicator function (oo (u) in Eq. limits the efficiency. We consider (2.()))
rare event. Since importance sampling is commonly used to handle Monte Carlo estimations involving rare events, we adopt the method of importance sampling on parameterized proposal distribution Qu for tractable Monte Carlo integration, and the estimator can be derived as:
12 Σ oy. ou let oy, (u) p(u) q(u)
P(Y.) = Eu [oy, (u))
()))))
(4)
where for each sample u ~ Qu. it is required that if p(u()) > 0, then (() > 0. The density ratio (qu)/q(u) is also termed as the importance weight and is usually assumed be bounded. In this section, we assume that our derivations are based on an exogenous distribution defined over continuous space. Naturally, similar conclusions hold for discrete exogenous distributions.
3.1
An Optimizable Variance Upper Bound
In importance sampling. a proposal distribution with low variance regarding σy, is crucial, as implies that fewer samples are required to achieve the same effect, leading to higher sampling efficiency. Previous works have attempted to constrain variance, such as providing upper bounds
18 49], introducing multiple proposals [16] [2], and learning from sampled data [72 83 11 65).
Given that the samples are i.i.d., with the same sample size n, a lower variance of ay. implies a lower estimator variance:
n ov.((ul Eu loy, (u)] - ()).,,
Vu-Qu (ay. (u)) = n
(5)
It is evident that this formulation can be directly optimized, and there are several methods employing
similar approaches (e.g., transforming it into cross-entropy or x2 divergence), as shown in [32] 68]
This process is called proposal learning, which ultimately yields a proposal distribution Qu that
suitable for the estimator.
The learned Qu is typically "one-off", rendering it unsuitable for other estimators. This necessitates multiple rounds of proposal learning to estimate an La expression that involves various counterfactual events, which remains inefficient. To address this issue, we propose replacing the proposal distribution Qu with a conditional proposal distribution QUIY.- Intuitively, the proposal distribution Quly. corresponding to y. E Qy, should be concentrated on the support Pu (5(y))). thus allowing for reuse across different estimators.
as illustrated in App. B. Additionally, we find that when certain conditions are met in Vi,, for any y. E V., the proposal distribution Quly. shares a common variance upper bound:
Building on this idea, we extend the methods from 32 68
Theorem (Variance Upper Bound). Let σy. (u) (p(u)/ //////////)))/) ( (u), where q(u yy) denote the density of the proposal distribution Quly.- and let Y. Yu(u) be the potential response w.r.t. u. for any y. E V., there exists k 2 1 such that 1/x & ((u)/q(u ly.) S k holds almost surely on the support lu ()..,, then for any Quly. where y. E V.:
V Vu-Quy. [ay. (u)] < -uu-Pu qu) Y. (u))] + where the constant is solely dependent x and Pu.
(6)

---

It is noted that the expected term in Eq. 6 bears resemblance to cross cross-entropy, albeit not strictly. Intuitively, the original indicator function is now implied as potential responses within the condition, and optimizing this expected term can encourage the density of each proposal distribution at specific locations to match that of the exogenous distribution as closely as possible. Therefore, this approach is named Exogenous Matching (EXOM).
The boundedness condition on the importance weights in Thm. 1 likely be violated for V with an infinite support set. We further discuss this issue in App. A.2 where we: i) provide two relaxed versions of Thm. based on concentration inequalities; ii) demonstrate how to construct guard proposal distribution that necessarily satisfies the boundedness condition.
3.2 Learning and Inference
Since any L3 expression Q is composed of terms P().). and assuming that is expressible by a finite set of P().). we can consider all y/ involved in Q as outcomes of a stochastic process. To allow different forms of V.. we extend Thm. to general cases, making it potentially applicable to the estimation of any counterfactual probability term P().) in an L3 expression 2, with only one proposal learning process required.
Generalize Thm.
Definition 1 (Stochastic Counterfactual Process). The collection of counterfactual variable sets 2) (Y E S] is referred to a stochastic counterfactual process w.r.f. the state space S. A state & E S is set composed of triplets (Y, X, x), where X, Y C V represent the intervened and observed variables, respectively, and denotes the intervened values. Each state corresponds to set of counterfactual variables Y U(Y,X.))== Yx-
Corollary 1 (Expected Variance Upper Bound). Let Qs denote arbitrary distribution defined over the state space S of a stochastic counterfactual process, and let P(- denote an arbitrary distribution defined over the o-algebra 2n/ corresponding to a set of counterfactual variables y given a state
E S. Py. is the joint distribution induced by p(s) and Qs. If the conditions in Thm. are met for
any a and any s E S, then for any Quy- where y+ E y and any S:
E Eyg wpy (Va-Quyy [oy. (w) S
[EwwPo [log ( Y! (u)) + c,
where the constant e is the same as in Thm.
Sampling and optimization Disregarding the constant terms of the expected variance upper bound in Cor. we obtain the optimization objective as follows:
[Eq.
[log Y (w)]]
arg min EanQQ
Based on our assumptions, Pu and Y (u) are known, while the prior distribution of states Qs needs be specifically designed for the particular expression in C3. The problem is now reformulated as conditional distribution learning problem, and we can model this conditional distribution using neural networks. The complete learning algorithm is presented in Fig. in App. A.3]
Inference Given the conditional proposal QUJY, we can take the expectation of the importance sampling results over V, to intuitively further enhance robustness. This leads to another unbiased estimator in the form of multiple importance sampling [16] 271:
P(Y.) By.~Qr. [Eu-QuN [oy. (u)l
Σ o 2. (1110)((1 (u))).
where (y.(y.](u) (p(u)/ (((/y.)) Inoo))(u). and ( ~ Qr u(() ~ Quys. The distribu- tion Qy could be any distribution with support covering y
Conditioning We attempt to find a function that maps y. to the vectorized parameters fy. required for the proposal distribution yu The information provided by yx (i.e., the counterfactual event Y. E {y+}) can be interpreted as set of 4-tuples, where each 4-tuple (Y, X, x) expresses
5

---

m
mi
m2
m
Mask Compute
Forward
Thm. 3
Encuding
Eq. 10
y+
h
9
0y. QUIy-
Overview of the conditioning and masking process. y. the input to the entire process, m represents the inferred mask, and the vectorized parameters fy. of the proposal distribution Quly. are the output. Different colors represent information from different submodels. Both and represent neural networks.
Figure 2:
counterfactual information within the same submodel Mx including the observed variables Y and their values and the intervened variables X and their values x. We represent it as vector:
C(Y y. (y Ux Y UX B (Y) w(X),
(10)
where @ denotes vector concatenation. The functions π and w map inputs to vectors, with each component of the vector corresponding to an exogenous variable X; V, defined as follows:
Zi if X X 0, otherwise and (X)
1,
if X, X 0, otherwise
T. (x; X)
(11)
To encode the entire set {(Y,, V1. x1), 22, X2, x.)...... we introduce two functions, and Here, h serves as an encoder, mapping the vector of a single 4-tuple o xxxx).,x) to latent encoding, while I acts as an aggregator, capturing the permutation invariance of the set. The vectorized parameters 8y, are inferred through the interplay of these two functions:
By. = (f M(x)))x.Xxx)) Y [xx.] YY,)).
(12)
We choose model h using multilayer perceptron (MLP), and model g function that ensures permutation invariance for encoding aggregation, such as summation or weighted summation, where opt for attention [96] to enhance expressiveness.
3.3 Injecting Markov Boundaries
This section will explore how to inject available structural prior knowledge (specifically, counterfac- tual Markov boundaries) into neural networks used for conditional encoding, intuitively improving the quality and efficiency of learning.
Counterfactual Markov boundary
Markov boundaries 75 have been used in feature selection
3 87 [03], since they reveal the local causal structure of variables, where all elements strongly correlated with the variable (2). Just as in feature selection. intuitively we can enhance the performance of distribution learning by masking redundant information. For modeling distributions defined over exogenous variables and conditioned counterfactual variables, the structural prior information used here named counterfactual Markov boundaries:
[I]]] 64 [04] and causal discovery
Definition 2 (Counterfactual Markov Boundary). For exogenous variable U, U and set of counterfactual variables Y, along with their joint distribution P. If U, is independent of Y X given X under i.. U, _ (Y \ X) X, then X is termed a Markov blanket U, on Y. The collection of all Markov blankets of U, on Y is denoted as (Y). If X E B (Y) is a Markov blanket of U, on Y. and there exists no X' Ç such that X' E (Y(Y), then X termed a (counterfactual) Markov boundary of U, on Y, denoted as B,(Y)
Methods for learning Markov boundaries based on independence data distribution (e-g.
3
87 103) cannot be directly applied to learn counterfactual Markov boundaries unless exogenous variables are observable. Moreover, these methods may have low efficiency due to their reliance on independence tests. To obtain the counterfactual Markov boundary, in this work, we assume that

---

another representation of SCM graphs called augmented graphs is known. Compared to causal graphs, augmented graphs incorporate exogenous variables U into the nodes, then remove bidirectional edges and add unidirectional edges (U, V,) for all V, E U and Uj- E Uv.
When the SCM is recursive, the augmented graph of any submodel is necessarily acyclic due to the closure of recursiveness under (perfect) intervention [8], which implies that d-separation [75] always holds. We additionally assume faithfulness (2]. which means that only the variables that are d-separated in the graph are independent in the distribution, allowing us to directly compute counterfactual Markov boundary for any exogenous variable by graph algorithm.
Under the above assumptions, the following theorem indicates that the counterfactual Markov boundary of U is independent across different submodels. This not only reduces the search space of the Markov boundary, but also aligns with Eq. 12, where we encode counterfactual information according to each submodel.
Theorem 2 (Counterfactual Markov Boundary Independence), If = = corresponds to different submodel Mx.. then for each U3 E U. there exists Markov boundary BY.Y.) - B, s(x:)) on Y., where B, (px))) is a Markov boundary on Mx.-
Di Y Yxx and each (/.))
Another theorem demonstrates how to obtain counterfactual Markov boundaries via d-separation, and indirectly proves their uniqueness:
Theorem 3 (Counterfactual Markov Boundary on Graph). For an exogenous variable U, E U and a counterfactual variable set Y from the submodel Mx. the counterfactual variable Yx E (Y if and only if Yx V s2 Un Yx (Yx), i.e., when given Yx \ (Yx), Yx and Uj are not d-separated on 9x. where 9% is the augmented graph induced from the submodel MMx
The graph algorithm to determine counterfactual Markov boundaries B, (Y.) for U, E U is derived simply by combining Thms. 2 and with the d-separation criterion as a subroutine, which can be found in Fig. 7 in App. A.4
Masking We then inject the counterfactual Markov boundary into the neural network used for conditioning by vectorizing it as weight mask, as briefly illustrated in Fig. 2 depicting the entire conditioning and masking process. Specifically, let Oy... i be any parameter corresponding to exoge- nous variable U, E in Quly. According to the description in Thm. 2 the Markov boundary
on Y. U 1=1 Yixx,) can be precisely decomposed into the union of the Markov boundaries on each Y sxaal This corresponds exactly to encoding the information for counterfactual events on each submodel Mx, as depicted Eq. 10 Therefore, we can modify [q. 12 follows:
y 9 (1 hy (CLY, (C.Y,... m,;) Y [|x]] Y.)).
(13)
In particular, m.j = (B, (Y))) and (Y are obtained through algorithm derived by Thm. The correspondence between the parameters 8y.., and U, exists in the specific design of the model, such as the mean and covariance of each component in GMMs or the element- desent-wise transformations in normalizing flows. As used in 31 81 13,, h is an MLP that allows weight masking, such that for all j, the i-th component of Vxhy (X, miy) * if and only if m,j = 1. This ensures that Oy... is only related to the Markov boundary of U, on y..
4 Related Works
Importance sampling with normalizing flows recent years, the combination of normalizing flows and importance sampling has been widely employed in tasks related to Monte Carlo integration.
For instance, a plethora of literature (70, 99 10T 53 60 TO efficiently compute the partition function of energy. Some studies [68 30 36) also employ this
67 utilizes this combination to
method to solve integration problems of complex functions in nigh-dimensional spaces. When the
integrand involves rejection sampling [6 89) particle events [29] 9 88) or failure events [33] 20),
this method can improve the efficiency of sampling rare events. This combination is also applied for 79 as alternative integrating over the intractable distribution.
posterior inference [19
Identifiable neural proxy SCMs Our approach applicable to pre-trained neural proxy SCMs, where causal mechanisms modeled neural networks. Increasing attention has been paid to the
7

---

identifiability of these models. BGM 69 combines bijections to demonstrate the identifiability of
causal mechanisms in several special cases. CausalNF 42 builds upon the conclusions of [066]
to prove the identifiability of its causal mechanisms up to invertible functions. NCM [08] extends the findings of their previous work [107] and proposes a sound and complete algorithm to identify counterfactual queries.
Arbitrary conditioning In dealing with exponentially many conditional distributions, our work establishes connection with another unsupervised learning task, known as arbitrary conditioning. The task of this paper can be viewed as extending arbitrary conditioning from a single observational distribution to multiple counterfactual distributions. Relevant works include VAE-AC [40) based on VAE, NC 77 based on GAN, ACE [90] based on energy models, and ACFlow [59] based on normalizing flows. The most similar to our work is [91], which matches the posterior distribution of VAE, akin to our matching on the exogenous distribution.
5 Experiments
Experiment settings and metries We first define the prior distribution Qs over the state space of a stochastic counterfactual process 2) and train according to Eq. 8) Subsequently, based on Qs. we construct a distribution Py. concerning counterfactual events for the evaluation of the estimator (Eq. 4). Inspired by the metrics for rare event sampling effectiveness, we utilize Effective Sample Proportion (ESP) and Failure Rate (FR) with threshold m to measure the performance the importance sampling, defined as:
m())) m) m
ESP- Ey, Py. [n()))] and FR ERREEy-o
Ey.
(14)
oeheewise
where v()))= -E: I ().) (u()))) //n is the proportion of effective samples (nonzero indicator) to estimate the probability P():.), which indirectly reflects the sampling efficiency of single estimate; this is equivalent to the success rate in rare event sampling. Therefore, ESP reflects the overall efficiency of the sampling method in supporting TThm. I while FR reflects its overall effectiveness in the state space to support Cor. 1 and Eq. 8. For further discussion on metrics, see App. C.I
Specifically, the following two types of stochastic counterfactual processes are involved in subsequent experiments: i) D with state space 5 {s |s| = k} (that is, with & submodels) and prior distribution Q where the indicators for intervention and observation follow Bernoulli distributions, and their values follow endogenous distributions; ii) Di, where Q E {ATE. ETT, NDE, CUDEJ, the design of state space S and prior distribution Qº are detailed in App. C.2
During evaluation, the distribution of counterfactual events Py. is contingent upon the type of counterfactual variables: for discrete, we focus on counterfactual event V E }y= for continuous, we focus on counterfactual event V.- E (y.),, where (y=)) is a cube with side length 0.02 centered at y.. The latter can be further used to estimate the counterfactual density in continuous space, as discussed in App. D.2]
Our experiments involve the following fully spec- ified SCMs, which fall into three categories: i) Markovian diffeomorphic (42), where there is
One Run
Initial State
Final State
1.0
no hidden confounder and causal mechanisms are diffeomorphic, including SIMPSON-NLIN, TRIANGLE-NLIN and LARGEBD-NLIN; ii) Semi-Markovian continuous [107). where hidden confounders may exist and endogenous variables are continuous, including M and NAPKIN; iii) Re- gional canonical [08], where hidden confounders may exist and endogenous variables are discrete and finite, including FAIRNESS and FAIRNESS- XW. The consistent performance observed in ex- periments across these different types of SCMs will of the proposed method.
0.4
0.2
11.4
0.2
000
E.I.
I.I.
Figure 3: LL (negative Eq. 8) and ESP, FR on SIMPSON-NLIN. As LL increases, ESP in- creases while FR decreases, until convergence. serve as a reference for evaluating the robustness
8

---

Table 1: Comparison of RS, CEIS, NIS, and EXOM (ours) across 3 different SCMs, with 1,3, Among the three SCMs, SIMPSON-NLIN is Markovian diffeomorphic, NAPKIN is Semi-Markovian continuous, and FAIRNESS- XW is Regional canonical. The results are averaged over runs. Higher ESP and lower FR indicate better performance.
SIMPSON-NLIN
NAPKIN
FAIRNESS-XW
Model
ESP f
FR
ESP
FR 4
ESP 1
FR +
RS
0,008
0.971
0.006
0.967
0.156
0.065
CEIS
0.037
0.809
0.017
0.911
0.255
0.532
NIS
0,008
0.961
0.007
0.965
0.193
0.527
EXOM[GMM]
0.117
0.245
0.039
0.476
0.358
0.009
EXOM[MAF]
0.581
0.005
0.306
0.066
0.339
0.007
RS
0.000
1.000
0.000
1.000
0.038
0.281
CEIS
0.000
1.000
0.000
1.000
0.122
0.704
NIS
0.000
1.000
0.000
1.000
0.116
0.686
EXOM[GMM]
0.024
0.500
0.011
0.613
0.247
0.056
EXOM[MAF]
0.606
0.075
0.397
0.183
0.261
0.043
RS
0.000
1.000
0.000
1.000
0.030
0.368
CEIS
0.000
1.000
0.000
1.000
0.106
0.727
NIS
0.000
1.000
0.000
1.000
0.094
0.724
EXOM[GMM]
0.020
0.497
0.009
0.644
0.231
0.084
EXOM[MAF]
0.698
0.031
0.482
0.094
0.237
0.070
Markos Boundary Masked
Madeol
GMM
MAF
NICE
GMM
NICE
SOSPF
SIMPSON NLIN
LARGEBD NLIN
M
NAPKIN
Figure 4: Ablation study for Markov boundaries on 4 different settings SCMs: (a) SIMPSON-
NLIN, (b) LARGEBD-NLIN, (c) M. (d) NAPKIN. A higher ESP signifies greater sampling efficiency. In most cases, EXOM with Markov boundaries masked (orange bar) exhibits superior performance compared to when the Markov boundaries are not masked (blue bar).
Convergence We integrated Exogenous Matching with various density estimation models (GMM, MAF (73). NICE (24) SOSPF (4)) to conduct experiments the stochastic counterfactual process D (with |s| = 3) over different SCM settings. We present the results ESP and FR during the training process. As shown in Fig. (with all results detailed in App. C.5), both the ESP and FR metrics change as expected until convergence, along with the decrease in training loss [5.. 8 (equivalent the increase in L.I in 3). This validates the effectiveness of Thm. aand and empirically suggests that the conditions we assumed are generally not violated under our settings.
Comparison We compare our proposed method with three different but related sampling methods: i) Rejection Sampling (RS): ii) Cross-Entropy based Importance Sampling (CEIS, (32)): iii) Neural Importance Sampling (NIS, (68). To apply the latter setup, we extended them as detailed in Appp B.I As shown in Tab. our method outperforms other approaches in both continuous and discrete cases. Specifically, since MAF exhibits stronger representational capabilities than GMM, it is expected that EXOM using MAF as the conditional distribution model generally outperforms EXOM using GMM.
methods under the same experimental

---

Table 2: Estimation of counterfactual densities on CausalNF and counterfactual effects on NCM.
Here, "O" represents the original SCM, and "P" represents the proxy SCM. For SIMPSON-NLIN the proxy SCM is Ca and the metric used is FR ("-" indicates FR equals 1); whereas for FAIRNESS, the proxy SCM is NCM, and the metric used is the average bias w.r.t. the ground truth.
The subscript denotes the 95% CI error bound over 5 trials. For more details, see[App C.9
SIMPSON-NLIN
FAIRNESS
Method
SCM
|8|
[s]
=8
ATE
ETT
NDE
CtfDE
Q
0.89p.01
0.01g.013
0.01.
0.01g.015
0.01
RS
P
0.90p.012
0.01g.013
0.01o.021
0.01o.014
0.01p.023
O
0.00p.00s
0.02o.015
0.01p.021
0.04g.095
0.06o.168
0.04g.131 131
0.070.236
EXOM[MAF]
P
0.01p.uus
0.40g.012
0.61p.v16
0.01g.013
0.01p.024
0.01g.013
0.01p.0s4
EXOMINICHI
O
0.00p.onG
0.020.016
0.01o.0a6
0.01g.018
0.01
0.030
0.01a.020
0.01
P
p.006
0.40g.013
0.61
0.01
0.01 0.021
0.01g.012
0.01
To investigate the impact of injecting Markov boundaries, we present the improvement of ESP after injecting Markov boundaries during the training process compared to the default scenario. We conducted ablation experiments on 4 continuous SCMs and various density estimators. As shown in InFFig. 4 the use of Markov boundaries as masks significantly improves the performance of EXOM under various settings. We also observed that when masks are used, the width of the hidden layers of neural networks in conditioning affects practical performance, which will be discussed in App. C.8
Ablation
Counterfactual Estimation on Proxy SCMs We combine identifiable proxy SCMs with EXOM for counterfactual estimation. We employ CausalNF 42 for SIMPSON-NLIN (based on the stochas- tic process 2 and NCM [108] for FAIRNESS (based the stochastic process 2) For the combination of CausalNF and SIMPSON-NLIN, we estimate the counterfactual probability on cubes (y))) with side length = 0.02 centered around 1024 randomly sampled y. from Py. (this can be used further to estimate counterfactual density), measure the sampling FR, and employ a dimension- regularized 95% CI error bound. For the combination of NCM and FAIRNESS, we estimate different counterfactual queries, measure the average bias of the query results relative to the ground truth, and use 95% CI error bound to demonstrate the unbiasedness of the estimates.
The results in Tab. 2 empirically demonstrate the effectiveness of EXOM in counterfactual density and effect estimation tasks compared to RS. Specifically, when estimating the cube áy (y.) (equivalent to counterfactual density), RS almost fails in high-dimensional settings, while EXOM exhibits good sampling efficiency and lower error compared to RS. When estimating counterfactual queries, the bias and error of EXOM are similar to those of RS (partly due to the lower sampling difficulty in discrete cases). Further experimental details and analysis of the conclusions can be found in in|Ap.. C.9
Conclusion
We propose Exogenous Matching for tractable estimation of expressions in C3 in general settings. Specifically, leveraging importance sampling, we find variance upper bound for estimators poten- tially applicable to any relevant counterfactual probability (Cor. 1). and then introduce an optimization objective (Eq. 8) amenable to sampling, transforming the problem into a conditional distribution learning problem. Theoretical and empirical results demonstrate that this approach can be combined with identifiable proxy SCMs to address practical problems. We also explore injecting Markov boundaries as prior knowledge and empirically validate effectiveness in several scenarios.
Limitations i) The proposed method does not directly estimate counterfactuals through the available distribution but requires a partially specified SCM. This work explores its effectiveness in combination with identifiable proxy SCMs, which are still an active area of research; ii) It is important to emphasize that the "general case" mentioned in this paper does not cover all scenarios. For instance, our estimand, an La expression Q. requíres the assumptions of finite submodels and approximability by a finite number of counterfactual probabilities; iii Furthermore, the specific structures and tricks demonstrated in this paper (such as the injection of Markov boundaries) may not be applicable to SCMs with general parameter settings, which warrants further investigation.
10

---

Acknowledgments and Disclosure of Funding
This work was supported in part by National Key R&D Program of China (No. 20222D0120302.
References
[1] A. Agrawal, D. R. Sheldon, and J. Domke. Advances in black-hox vi: Normalizing flows, importance weighting, and optimization. Advances in Neural Information Processing Systems, 33:17358-1736, 2020.
[2] C. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos, Local causal and markov blanket induction for causal discovery and feature selection for classification part i: Algorithms and empirical evaluation. Journal of Machine Learning Research, 11(7):171 2344, 2010.
[3] C. Aliferis, Statnikov, I. Tsamardinos, Mani, and D. Koutsoukos, Local causal and markov blanket induction for causal discovery and feature selection for classification part it: Analysis and extensions. Journal of Machine Learning Research, 11(8):235- 2844244,,,, 2010.
[4] S. Balgi, J. M. Peña, and A. Daoud. Personalized public policy analysis in social sciences using causal-graphical normalizing flows. Proceedings of the AAAI Conference on Artificial Intelligence, 36(11):111 11100111111111111111111111111111 1118111881,,,1111111111111111881 Jun. 2022.
[5] E. Barcinboim, J. D. Correa, D. Ibeling, and T. Icard. On Pearl's Hierarchy and the Foundations of Causal Inference, page 507-556. Association for Computing Machinery, New York, NY, USA, 1 edition, 2022.
[6] M. Bauer and A. Mnih. Resampled priors for variational autoencoders. Chaudhuri and M. Sugiyama, editors, Proceedings of the Twenty-Secon International Conference on Artificial Inteliience and Statis- tiex, volume 89 of Proceedings of Machine Learning Research, pages 66-75, PMLR, 16-18 Apr 2019.
[7] M. Belghazi, M. Oquab, and D. Lopez-Paz. Learning about exponential amount of conditional distributions. H. Wallach, H. Larochelle, A. Beygelzimer, d'Alché-Bue, E, Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[8] S. Bongers, P. Forr'e e, J. Peters, and J. M. Mooij. Foundations of structural causal models with cycles and latent variables. The Annals of Statistics, 2016.
[9] Bothmann. T. JanBen, M. Knobbe, T. Schmale, and Schumann. Exploring phase space with Neural Importance Sampling. SciPost Phys., 8:069, 2020.
[10) Y. Cao and E. Handen-Eijnden Learning optimal flows for non-equilibrium importance sampling. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 12352 12364. Curran Associates, Inc., 2022.
[11] O. Cappé, R. Douc, A. Guillin, J.-M. Marin, and C. P. Robert. Adaptive importance sampling in general mixture classes. Statistics and Computing, 1844):447-459 Dec 2008.
[12] P. Chao, P. Blöbaum, and S. P. Kasiviswanathan. Interventional and counterfactual inference with diffusion models. CoRR, abs/2302. 000600 2023.
[13) A. Chen, R. Shi, X. Gao, R. Baptista, and R. G. Krishnan. Structured neural networks for density estimation and causal inference. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 66438-66450. Curran Associates, Inc., 2023.
[14] D. Cheng. Z. Xu, Li, Liu, J. Liu. and D. Le. Conditional instrumental variable regression with representation learning for causal inference. CoRR, abs/2310.01865. 2023.
[15] M. Cheng, X. Liao, Q Liu, B. Ma. J. Xu, and B. Zheng. Learning disentangled representations for counterfactual regression via mutual information minimization. In Proceedings the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR '22, page New York, NY, USA, 2022. Association for Computing Machinery.
[16] J.M. CORNUET, J.M MARIN, A. MIRA, and C. P. ROBERT. Adaptive multiple importance sampling. Scandinavian Journal of Statistics, 39(4):798 812, 2012.

---

[17] J. Correa, S. Lee. and E. Bareinboim. Nested counterfactual identification from arbitrary surrogate experiments. In M. Ranzato, A. Beygelzimer. Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances Neural Information Processing Systems, volume 34, pages 6856-6867. Curran Associates, Inc., 2021.
[18] C. Cortes, Y. Mansour, and M. Mohri. Learning bounds for importance weighting. In Lafferty, C. Williams, J. Shawe- Taylor. R. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 23, Curran Associates, Inc., 2010.
[19] C. Cundy and S. Ermon. Flexible approximate inference via stratified normalizing flows. In J. Peters and D. Sontag, editors, Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pages 1288-1297. PMLR, 03-06 Aug 2020.
[20] A. Dasgupta and E. A. Johnson. Rein: Reliability estimation via importance sampling with normalizing flows. Reliability Engineering & System Safety, 242: 242:109729, 2024.
[21] M. Dax, S. R. Green, J. Gair, M. Pürrer, J. Wildherger, J. H. Macke, A. Buonanno, and B. Schölkopf. Neural importance sampling for rapid and reliable gravitational-wave inference. Physical Review Letters, 130(17):171 2023.
[22] E. De Brouwer. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 35213-35225. Curran Associates, Inc., 2022.
Deep counterfactual estimation with categorical background variables. In S. Koyejo,
[23] F. De Sousa Ribeiro, T. Xia, M. Montciro, N. Pawlowski, and B. Glocker. High fidelity image counterfac- tuals with probabilistic causal models. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 3325-7425, 23-29 Jul 2023.
[24] L.. Dinh, D. Krueger, and Y. Bengio. NICE; non-linear independent components estimation. In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Workshop Track Proceedings, 2015.
[25] G. Duarte, N. Finkelstein, D. Knox, J. Mummolo, and I. Shpitser. An automated approach to causal inference in discrete settings. Journal of the American Statistical Association, 0(0):1-16, 2023.
[26] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios. Neural spline flows. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, Buc, F. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[27] Elvira, L. Martino, D. Luengo, and M. F. Bugallo. Generalized Multiple Importance Sampling. Statistical Science, 34(1):129 - 155, 2019.
[28] I. R. Fulcher, I. Shpitser. S. Marealle, and E. J. Tchetgen Tchetgen. Robust Inference on Population Indirect Causal Effects: The Generalized Front Door Criterion. Journal of the Royal Statistical Society Series Statistical Methodology, 82(1):199 214, 11 2019.
[29] C. Gao, S. Höche, J. Isaacson, C. Krause, and H. Schulz. Event generation with normalizing flows. Phys. Rev. D, 101:076002, Apr 2020.
[30] C. Gao, J. Isaacson, and C. Krause. i- flow: High-dimensional integration and sampling with normalizing flows. Machine Learning: Science and Technology. 144)0450233, oct 2020.
[31] M. Germain, K. Gregor, I. Murray, and H. Larochelle. Made: Masked autocncoder for distribution estimation. In F. Bach and D. Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 881-889, Lille, France, 07-09 Jul 2015. PMLR.
[32] S. Geyer, I. Papaioannou, and D. Strauh. Cross entropy-based importance sampling using gaussian densities revisited. Structural Safety, 76:15-27, 2019.
[33] T. Guo, H. Wang, J. Li, and H. Wang- Sampling -bsed adaptive design strategy for failure probability estimation. Reliability Engineering & System Safety, 241:109664, 2024.
[34] Han, Y. Chen, and A. Darwiche. On the complexity of counterfactual reasoning. In E. Elkind, editor. Proceedings of the Thirty Second International Joint Conference on Artificial Intelligence, IJCAI-23, pages 5676-5684 International Joint Conferences on Artificial Intelligence Organization, 8 2023. Main Track.
12

---

[35] J. Hartford, G. Lewis, K. Beyton-Brown, and M. Taddy. Deep IV: A flexible approach for counterfactual prediction. In D. Precup and Y. W. Teh, editors, Proceedings the 34th International Conference on Muchine Learning, volume of of Proceedings of Machine Learning Research, pages 1414-1423, PMLR, 06-11 Aug 2017.
[36] T. Heimel, R. Winterhalder, A. Butter, J. Isaacson, C. Krause, F. Maltoni, O. Mattelacr, and T. Plchn. Mednis-neural multi-channel importance sampling. SciPost Physics, 15(4):141, 2023.
[37] c. Hizh, S. T. John, A. T. Juuti, T. T. Saarinen, K. H. Pietiläänen, and P. Marttinen. Causal modeling of policy interventions from treatment- sequences. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors. Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Machine Learning Research, pages 13050-13088 PMLR, 23-29 Jul 2023.
[38] D. G. Horvitz and D. J. Thompson. A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association, (76260))633-688) 1952.
[39] C.-W. Huang, D. Krueger, A. Lacoste, and A. Courville. Neural autoregressive flows. J. Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2078-2087. PMLR, 10-15 Jul 2018.
[40] O. Ivanov, M. Figurnov, and Vetrov. Variational autoencoder with arbitrary conditioning. Interna- tional Conference on Learning Representations. 2019.
[41] P. Jaini, K. A. Selby, and Y. Yu. Sum-of-squares polynomial flow. In K. Chaudhuri and R. Salakhutdino; editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine I Learning Research, pages 3009-3018. PMLR, 09-15 Jun 2019.
[42] A. Javaloy, P. Sanchez-Martin. and I. Valera. Causal normalizing flows: from theory to practice. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages (8833-58864. Curran Associates, Inc., 2023.
[43] H. Jeong, J. Tian, and E. Barcinboim. Finding and listing front-door adjustment sets In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 33173-33185. Curran Associates, Inc., 2022.
[44] Y. Jung, J. Tian, and E. Bareinboim. Learning causal effects via weighted empirical risk minimization. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 1277-12709, Curran Associates, Inc., 2020.
[45] Y. Jung, J. Tian, and Bareinboim. Estimating identifiable causal effects through double machine learning. Proceedings of the AAAI Conference on Artificial Intelligence, 35(13):12133-1222 May 2021.
[46] Y. Jung, J. Tian, and E. Bareinboim. Estimating joint treatment effects by combining multiple experiments. A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 15451-18 15527, PMLR, 23-29 23-29 Jul 2023.
[47] A. Karimi, B. Schölkoff, and Valera. Algorithmic recourse: from counterfactual explanations to interventions. In Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency, FAceT 21, page 353-362, New York, NY, USA, 2021. Association for Computing Machinery.
[48] A.-H. Karimi, J. von Kügelgen, B. Schölkopf, and I. Valera. Algorithmic recourse under imperfect causal knowledge: probabilistic approach. In H. Larochelle, M. Ranzato, R. Hadsell. M. Balcan, and H. Lin, editors. Advances in Neural Information Processing Systems, volume 33, pages 265-277. Curran Associates, Inc., 2020.
[49] A. Katharopoulos and Fleuret. Not all samples are created equal: Deep learning with importance sampling. In Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2525-2534. PMLR, 10-15 Jul 2018.
[50] I. Khemakhem, R. Monti, R. Leech, and A. Hyvarinen. Causal autoregressive flows. In A. Banerjee and K. Fukumizu, editors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 5520-3522 PMLR, 13 15 Apr 2021.
13

---

[51] L. Klein. A. Krämer, and F. Noé. Equivariant flow matching. In A. Oh. T. Naumann, A. Globerson. K. Saenko, M. Hardt, and S. Levine, editors. Advances in Neural Information Processing Systems 36: Annual Conference on Neural Information Processing Systems 2023, NeurlPS 2023, New Orleans, LA. USA, December 10- 16, 2023, 2023.
[52] M. Kocaoglu, C. Snyder. A. G. Dimakis, and S. Vishwanath. CausalGAN: Learning causal implicit generative models with adversarial training. In International Conference on Learning Representations. 2018.
[53] J. Köhler, A. Krämer, and F. Noé. Smooth normalizing flows. In M. Ranzato, A. Beygelzimer, Y. N. Dauphin, P. Liang, and W. Vaughan, editors, Advances in Neural Information Processing Systems 34; Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, viriual, pages 2796-2809, 2021.
[54] A. Komanduri, Y. Wu, E. Chen, and X. Wu. Learning causally disentangled representations via the principle of independent causal mechanisms. CoRR, abs/2306.01213, 2023.
[55] A. Komanduri, Y. Wu, W. Huang, F. Chen, and X. Wu. Sem-vae: Learning identifiable causal representa- tions via structural knowledge. In 2022 IEEE International Conference on Big Data (Big Data), pages 1014-1023, 2022.
[56] M. J. Kusner, J. Loftus, C. Russell, and R. Silva. Counterfactual fairness. In I. Guyon, V. Laxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[57] A. Li and J. Pearl. Probabilitics of causation with nonbinary treatment and effect. Proceedings of the AAAI Conference Artificial Intelligence, (8(18))2046522047) Mar. 2024.
[58] H. Li, C. Zheng, Y. Cao, 7. Geng, Y. Liu, and Wu. Trustworthy policy learning under the counterfactual no-harm criterion. A. Krause, F. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Muchine Learning Research, pages 205988 PMLR, 23-29 Jul 2023.
[59] Y. Li, S. Akbar, and J. Oliva. ACFlow: A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 5831-5841. PMLR, 13-18 Jul 2020.
Flow models for arbitrary conditional likelihoods. In H. D. III and
[60] Liu, W. Gao, Z. Wang, and C. Wang. Pathflow: normalizing flow generator that finds transition paths. In J. Cussens and K. Zhang, editors. Proceedings of the Thirtt-Eighth Conference on Uncertainty in Artificial Intelligncce, volume 180 of Proceedings of Machine Learning Reseurch, pages 1232-1242. PMLR, 01-05 Aug 2022.
[61] Y. Liu, J. Wang, and Li. Edvae: Disentangled latent factors models in counterfactual reasoning for individual treatment effects estimation. Information Sciences, 652:1 19578, 2024.
[62] I. Loshchilov and F. Hutter. Decoupled weight decay regularization. International Conference on Learning Representations, 2019.
[63] J. Ma, R. Guo, S. Mishra, A. Zhang, and J. Li. Clear: Generative counterfactual explanations on graphs. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Curran Associates, Inc., 2022.
Information Processing Systems, volume 35, pages 25895-200 259007
[64] D. Margaritis. Toward provably correct feature selection in arbitrary domains. In Y. Bengio. D. Schu- urmans. Lafferty, Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009.
[65] L. Martino, V. Elvira, D. Luengo, and J. Corander. Layered adaptive importance sampling. Statistics and Computing, 27(3):599-627 623, May 2017.
[66] V. Melnychuk, D. Frauen, and S. Feuerriegel. Partial counterfactual identification of continuous outcomes with a curvature sensitivity model. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and Levine, editors. Advances in Neural Information Processing Systems, volume 36, pages 32020-32060. Curran Associates, Inc., 2023.
[67] L. I. Midgley, V. Stimper, G. N. C. Simm, B. Schölkopf. and J. M. Hernández. Lobato. Flow annealed importance sampling bootstrap. In The Eleventh International Conference on Learning Representations. 2023.
14

---

[68] T. Müller, B. Mewilliams, F. Rousselle, M. Gross, and J. Novák. Neural importance sampling, ACM Trans. Gruph., 38(5), oct 2019.
[69] A. Nasr- Eassr-Esahhny, M. Alizadeh, and D. Shah. Counterfactual identifiability of bijective causal models. In Krause, E. Brunskill, K. Cho, B. Engelhardd, S. Sabato, Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 25733-25754. PMLR, 23-29 Jul 2023.
[70] F. Noé, S. Olsson, J. Köhler, and H. Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(6457) 147, 2019.
[71] M. Oberst and D. Sontag. Counterfactual off- poolcc evaluation with Gumbel-mat structural causal models. In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference on Muchine Learning, volume of Proceedings of Machine Learning Research, pages 09-15 Jun 2019.
PMLR.
[72) M.-S. Oh and J. O. Berger. Adaptive importance sampling in monte carlo integration. Journal of Statistical Computation and Simulation, 4133-4):143-16. 1992.
[73] G. Papamakarios, T. Pavlakou, and I. Murray. Masked autoregressive flow for density estimation. In I. Guyon, U. V. Luxburg. S. Bengio. H. Wallach, R. Fergus. S. Vishwanathan, and R. Garnett, editors. Advances in Neural Information Procesxing Systems, volume 30. Curran Associates, Inc., 2017.
[74] N. Pawlowski, D. Coelho de Castro, and B. Glocker. Deep structural causal models for tractable counterfactual inference. H. Larochelle, M. Ranzato, R. Hadsell, Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 857-869. Curran Associates, Inc., 2020.
[75] J. Pearl. Probubilistic Reasoning in Intelligent Systems. Morgan Kaufmann, San Francisco (CA), 1988.
[76] J. Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009.
[77] J. Pearl and D. Mackenzie. The Book of Why: The New Science of Cause and Effect. Basic Books, Inc.. USA, 1st edition, 2018.
[78] D. Plecko and E. Bareinboim. Causal fairness analysis. CoRR, abs/2207.11388, 2022.
2022
[79] D. Prangle and C. Viscardi. Distilling importance sampling for likelihood free inference. Journal of Computational and Graphical Statistics, 1264):146111471, 2023.
[80] D. J. Rezende, G. Papamakarios, S. Racaniere, M. Albergo, G. Kanwar, Shanahan, and K. Cranmer. Normalizing flows on tori and spheres. H. D. III and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages
PMLR, 13-18 Jul 2020.
[81] F. Rozet, F. Divo, and S. Schnake. Zuko: Normalizing Flows in PyTorch, 2024.
[82] D. B. Rubin. Bayesian inference for causal effects: The role of randomization. The Annals of Statistics, 6(1 6(1):34-58, 1978.
[83] R. Rubinstein. The cross-entropy method for combinatorial and continuous optimization. Methodology And Computing In Applied Probability, 1(2):127-190. Sep 1999.
[84] P. Sanchez and S. A. Tsaftaris. Diffusion causal models for counterfactual estimation. In First Conference on Causal Learning and Reasoning. 2022.
[851 A. Shah, R. Dwivedi, D. Shah, and Womell. On counterfactual inference with unobserved confounding. NeurlPS 2022 Workshop on Causality for Real-world Impact, 2022.
[86] I. Shpitser and J. Pearl. Complete identification methods for the causal hierarchy. Journal of Machine Learning Research,
2008
[87] A. Statnikov, N. I. Lytkin, J. Lemeire, and C. F. Aliferis. Algorithms for discovery of multiple markov boundaries. Journal of Machine Learning Research, 14(15):499-56,, 2013.
[88] B. Stienen and R. Verheyen. Phase space sampling and inference from weighted events with autoregressive flows. SciPost Phys., 10:038, 2021.
15

---

[89] V. Stimper. B. Schölkopf. and J. Miguel Hernandez-Lobato. Resampling base distributions of normalizing flows. In G. Camps-Vall Valls, F. J. R. Ruiz, and I. Valera, editors, Proceedings The 25th International Conference on Artificial Intelligence and Statisticx, volume 151 of Proceedings of Machine Learning Research, pages 4915 4936. PMLR, 28-30 Mar 2022.
[90] R. Strauss and J. B. Oliva. Arbitrary conditional distributions with energy. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems. volume 34, pages 752-763. Curran Associates, Inc., 2021.
[91] R. Strauss and J. B. Oliva. Posterior matching for arbitrary conditioning. In S. Koyejo, S. Mohamed, Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances Neural Information Processing Systems, volume 35, pages 18888-18099 Curran Associates, Inc., 2022.
[92] P. Minchez-Martin, M. Rateike, and I. Valera. Vaca:
Designing variational graph autoencoders for causal queries. Proceedings of the AAAI Conference on Artificial Intelligence. 36(7):8159 119-88899--16. Jun. 2022.
[93] J. Tan, J. Blanchet, and V. Syrgkanis. Consistency of neural causal partial identification, 2024.
[94] S. Tsirtsis, A. De, and M. Rodriguez. Counterfactual explanations in sequential decision making under uncertainty. Advances in Neural Information Processing Systems, 34:30127-30133 2021.
[95] S. Tsirtsis and M. Rodriguez. Finding counterfactually optimal action sequences in continuous state spaces. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 3220-3244
Curran Associates, Inc., 2023.
[96] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and L. Polosukhin. Attention is all you need. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS'17, 17, page 6000-6010, Red Hook, NY, USA, 2017. Curran Associates Inc.
[97] H. Wang, Fan, Z. Chen, H. Li, W. Liu, T. Liu, Q. Dai, Y. Wang, Z. Dong, and R. Tang. Optimal transport for treatment effect estimation. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 5404-5418. Curran Associates, Inc., 2023.
[98] M. Wienöhst, B. van der Zander, and M. Liškiewicz. Linear-time algorithms for front-door adjustment in causal graphs. Proceedings of the AAAI Conference on Artificial Intelligence, 38(18):00577-2054, Mar. 2024.
[99] P. Wirnsberger, A. J. Ballard, G. Papamakarios, S. Abercrombic, S. Racanière, A. Pritzel, D. Jimenez Rezende, and C. Blundell. Targeted free energy estimation via learned mappings. The Journal of Chemical Physics, 153(14), 2020.
[100] A. Wu, K. Kuang, B. Li, and F. Wu. Instrumental variable regression with confounder balancing. In K. Chaudhuri, S. Jegelka, I.. Song, C. Szepesvari, G. Niu, and S. Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 24056-24077 PMLR, 17-23 Jul 2022.
[101] H. Wu, Köhler, and F. Noć. Stochastic normalizing flows. Advances in Neural Information Processing Systems, 33:5933-5944, 2020.
[102] P. A. Wu and K. Fukumizu. SMhetaS-intact-VAN Identifying and estimating causal effects under limited overlap. In International Conference on Learning Representations, 2022.
[103] X. Wu, B. Jiang, Wu, and H. Chen. Practical markov boundary learning without strong assumptions. Proceedings of the AAAI Conference on Artificial Intelligence, 37388 103988, Jun. 2023.
[104] X. Wu, B. Jiang. K. Yu, c. Miao, and H. Chen, Accurate markov boundary discovery for causal feature selection. IEEE Transactions on Cybernetics, 50(12):4983 4996, 2020.
[105] Y. Wu, L. Zhang, and X. Wu. Counterfactual fairness: Unidemtification, bound and algorithm. In Proceedings of the Eighty International Joint Conference on Artificial Intelligence, IJCAI-19. pages 1438-14400 1444, International Joint Conferences on Artificial Intelligence Organization, 2019.
[106] Q. Xi and Bloem- Bloem-Reddy. Indeterminacy in generative models: Characterization and strong identifiabil- ity. In F. Ruiz, J. Dy, and J.-W. van de Meent, editors, Proceedings of The 26th International Conference on Artificial Intelligence and Statisties, volume 206 of Proceedings of Machine Learning Research, pages 6912 6939. PMLR, 25 27 Apr 2023.
16

---

[107] K. Xia. K.-Z. Lee, Y. Bengio, and E. Bareinboim. The causal-neural connection: Expressiveness. learnability, and inference. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan. editors, Advances in Neural Information Processing Systems, volume 34, pages 10823-10836. Curran Associates, Inc., 2021.
[108] M. Xia, Y. Pan, E. Bareinboim. Neural causal models for counterfactual identification and estimation. In The Eleventh International Conference on Learning Representations, 2023.
[109] Xu, D. Cheng, Li, Liu, L. Liu, and K. Yu. Causal inference with conditional front-door adjust- ment and identifiable variational autoencoder. In The Twelfth International Conference on Learning Representations, 2024.
[110] M. Yang. Liu, Z. Chen, X. Shen, J. Hao, and J. Wang. Causalvae: Disentangled representation learning via neural structural causal models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pottern Recognition (CVPR), pages 9093-9602, June 2021.
[111] S. Yaramakala and D. Margaritis. Speculative markov blanket discovery for optimal feature selection. In Fifth IEFE International Conference on Data Mining (ICDM 005), pages 4 pp--, 2005.
[112] M. Zešević, D. S. Dhami, and K. Kersting, Not all causal inference is the same. Transactions on Machine Learning Research, 2023.
[113] J. Zhang and E. Bareinboim. Fairness in decision-making - the causal explanation formula. Proceedings of the AAAI Conference on Artificial Intelligence, 32(1), Apr. 2018.
[114] J. Zhang, J. Tian, and E. Bareinboim. Partial counterfactual identification from observational and experimental data. In K. Chaudhuri, S. Jegelka, L. Song. C. Szepesvari, G. Niu, and S. Sabato, editors, Proceedings of the 39th International Conference on Machine Learning. volume 162 of Proceedings of Machine Learning Research, pages 26548-2666 26558. PMLR, 17-23 Jul 2022.
[115] W. Zhang, L. Liu, and Treatment effect estimation with disentangled latent factors. Proceedings of the AAAI Conference on Artificial Intelligence, 35(12):1092) 10930, May 2021.
[116] Zhou, R. Bai, Kulinski, M. Kocaoglu, and D. Inouye. Towards characterizing domain coun- terfactuals for invertible latent causal models. In The Twelfth International Conference on Learning Representations, 2024.
[117] A. Zuo, S. Wei, T. Liu, B. Han, K. Zhang, and M. Gong. Counterfactual fairness with partially known causal graph. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 1238-1252. Curran Associates, Inc., 2022.
17

---

Appendix
Table of Contents
A
Proofs
19
A.
Countertactuals
19
A.2
Importance Sampling
20
Learning and Inference
27
A.4
Markov Boundary
28
B
Related Works
30
B.
Importance Sampling
30
B.2
Proxy SCMs
31
B.3
Arbitrary Conditioning
32
Injecting Prior Knowledge
C
Experiments
34
Metrics
34
Design of Stochastic Counterfactual Processes
34
C,3
Fully-specifics SCMs
36
C.4
Reproducibility
36
Convergence
Comparison
38
Combinaitons SCMs Density Estimation Models
38
C88
Ablation
................................
40
C.9
Counterfactual Estimation Proxy SCMs
43
Discussions
44
D.1
Broader Impacts
44
Counterfactual Density
44
Limitations
45
18

---

A Proofs
A.1
Counterfactuals
In this section, we demonstrate the validity of counterfactual probability-related concepts, providing intuitions based on derivations in discrete spaces, and proving lemmas to streamline subsequent derivations. Firstly, we establish the validity of counterfactual probabilities as defined in mn.. 2
Proposition 1. The probability measure defined by Eq. 2 (on recursive SCMs) adheres to the Kolmogorov axioms.
Proof. For any V, E 2y, Pr. (J+) S 1. Since I (a)().)(u) S the proof follows from the Monotonicity of Lebesgue integration:
0 == / OdPu
(n) -
IdPu
(15)
Sbu
Suu
where Pv, = Inn (no())) dPu, thus 0 S (Y.) 1.
(ii) Pr. (fy.) ) 1. Since SCM is recursive, given any u Qu. the values of exogenous variables determined, and hence all values in Y, are determined, that is, the potential responses (.
Since (u) E fty,, according to the definition of ll (arr( (a = 1, thus:
Px. (9xx))= = /
Du (u) dPu
IdPu =1 1.
(16)
Dau
(iii) Py. (U2 - Py. (D. where (
is countable (possibly countably
infinite) pairwise disjoint counterfactual events. For flu if (u) y since for any
j # i, (0 ny )o 9, it follows that (u) & y(o), thus 1 Aooo (u) = 1 Du 2l (u) = 1
and Σjti ( (u) = 0, hence.
I (eo (u) (u) dPu
Py. (.!))
Definition
(17)
121
= So (La0cos9/())= Σ Laucess dPu
(18)
jsti
- -/a (u) + I
Σ
(no))))((1))) (u) dPr
Linearity
(19)
= -/n Inuret)
(u) dPu
Zero measure
(20)
= = 1 (noc)oo((u))0 (u) + Σ
1 ( dPy
(21)
ifi su
Σl
121 Su (u)
(22)
Σ (O.!(
(23)
for i, (u) & ! then Σiz 1 "no(sS9 (u) =
w0 1 (u) = Therefore, both sides
(u
of the equation are zero, so it holds. Hence, we have proved Pr. (Uiz) y²! = Σiz1 Pr. ().
Given that axioms (i), (ii) and (iii) hold, which are the Kolmogorov axioms, it follows that Eq. satisfies the definition of probability measure.
From a measure-theoretic perspective, the indicator function in Eq. can be transformed into the Dirac measure, and vice versa:
19

---

Lemma 1. For any Y. E Σy, and any E Qu.
lo))) (v) = / y.
(u)-
(24)
Proof. Given Y. E M and u E Lu, we can consider Y, as u being constant, thus
(
Definition of Su ():.)
(25)
=Iy.(Y = (D.)
Let y. (u)
(26)
Definition of yy.
(27)
= -Son, (dv.(
Integral w.r.t. Dirac measure
(28)
//- dd).
(29)
where the integral fy. doy remains function of Then, substituting Lem. into Eq. 2 have:
P(Y.) / (u)())(u)dRR dPu
S / dōx. dPr
s..
Lem
(30)
[ / dPx..u,
(31)
where Py, satisfying dPy., dōx.(u (u) dPu is called the joint distribution of counterfactual variables and exogenous variables, which we abbreviate as the joint distribution.
In words, the joint distribution is determined by the exogenous distribution and the Dirac distribution, while the counterfactual probability the integral over the endogenous space after marginalizing the exogenous variables from the joint distribution.
A.2 Importance Sampling
Importance sampling estimator for counterfactual estimation where exogenous variables are continuous. Here, we will provide a complete derivation. For discrete cases, the Lebesgue integral can be expressed as summation, thus,
Eq. 4 4jolyy illustrates the case
P(Y.) = / (u)
(32)
Σ ((u))n)()))(u)) ueft
(33)
P(u)
Σ uetu
Q(u)
( Qu1)
(34)
P(u) Ae (u)
Eu-Qu
(35)
[{uu)
For continuous cases, the Radon-Nikodyy theorem is required. Suppose that the probability density functions of Pu and Qu exist, denoted by p(u) and q(u), respectively. According to the definition of probability density functions, p(u) = dPu/ / du and (u) = dQu/ du. Then
P().) / In)(r:
(36)
f Lno())) (u) dRu/du dQu dQu/
(37)
20

---

p(u) I (n)
(38)
Au g(u)
((u) [q(u)
Qu
(1)
(39)
We assume that the integrand is Lebesgue integrable. It can be observed that the difference between the discrete and continuous cases lies only in the importance weights. Therefore, subsequent deriva- tions, we will proceed in the continuous case, with the discrete case requiring only the substitution of the importance weights.
The expectation term can be estimated using Monte Carlo methods, resulting in the basic form the importance sampling estimator:
P().) Eu-Qu Jo (u) oy. (u (u
where oy. (u) (s) and each sample n Pu.
Variance of the importance sampling estimator
The variance of the estimator in Eq. 4
oy (u
n² no (u!
Linearity
(40)
n V u0 [ox;
i..d
(41)
V Qu 1oy, (u)]
the first equation)
n
Vu-qu Joy, (u)] can be further expanded as:
p(u) (
V Vu-Qu loy. (u)] V VunQu
(42)
[q()
((u)
27 (n)()))(u)
p(u) qqu)
BunQu
Eunqu
(43)
((u)
p²(u), (2( (1) -P²()))
=RanQu
(44)
-Se (p(U) dPu
(n) (u) P²().)
(45)
dQuu
p(u)
-av q(u)
Ino dPu - p2())))
(46)
p(u)
=Eu-Pu
l (Z)) P²(
(47)
[q(u)
= -Eu-Po [oy. (u)] - ()))).
the second equation)
The proposal distribution Qu only affects the term En-ro (ay, (u) in the variance Eq. The only difference between and the estimator P().)) En-Qu (ay. (u) is that En-Po oy. (u) is the expectation over the exogenous distribution Pu.
Optimization of Eq. 5 and the limitations Reducing variance implies improving estimation efficiency; other words, means achieving the same level of estimation accuracy with fewer samples. One direct method to optimize variance is by sampling from the exogenous distribution Pu using Monte Carlo methods to estimate the term Eg-Pt [oy: (u)] in the variance Eq. 5) and then minimizing the estimated value of this term using stochastic optimization techniques. Alternatively, as suggested by 68 this term can be transformed into an optimization problem involving the cross-entropy x divergence between the proposal distribution and the integrand, and then optimizing the values of these mismatches estimated by sampling under the proposal distribution Qu can also indirectly optimize variance. This optimization process is referred to proposal learning.
21

---

However, a learned individual proposal distribution is typically only applicable to single estimator, meaning can only estimate the probability for a single counterfactual event V. as the indicator function 1 (u: (u) in oy: (u) operates solely for V But L3 expressions can typically contain multiple counterfactual probability terms, and this limitation implies that, in the worst case, we need to perform separate optimization for each counterfactual probability, still affecting the efficiency of the counterfactual estimation task.
To address this issue, our solution is to replace the sampling proposal distribution Qu with a conditional proposal distribution QUIY., adding extra degree of freedom conditioned on Y.. Our overall idea is to let these different proposal distributions participate in the same optimization (i.e., sharing the same optimization objective) that a single optimization yields different proposal distributions applicable to different counterfactual events.
Additionally, to limit the scope of responsibility of each proposal distribution Quly. that is, to prevent too many proposal distributions from targeting the same counterfactual event V and conse- quently neglecting some counterfactual events we further impose a one-to-one correspondence between conditions and responsibilities for inference. Intuitively, this would make the proposal distribution QUly. corresponding to y. E fty focus on the support (o )). Given that for any counterfactual event there is an approximation V. R Uy. 5(y.), we can use the proposal distributions Quly. corresponding to these y. E V. to serve the estimation of P()).).
Variance upper bound When certain conditions are met in V. for any E V., the proposal distri- bution Quly. shares a common upper bound on variance, which is exactly the common optimization objective we seek:
Theorem (Variance Upper Bound). Let σy. (u)
(nn where q(u y-) denote the density of the proposal distribution QUly. and let Y. (u) be the potential response w.r.t. If for any y+ E V., there exists k 2 1 such that 1/x S p(u)/g(uu (y..) S x holds almost surely on the support lu ()+).. then for any Quly. where y. E V.:
Vu-Quy. lay. (u)] S -uuro [log q(u Y. Y.())) +
where the constant e solely dependent on is and Pu.
Proof. For a V. E Σy. that satisfies bounded importance weights (i.e., there exists k 2 1 such that 1/x S p(u)/q(uu y=) < ε holds almost surely on the support ty ().)), the following inequality holds for any y+, y² D. holds almost surely:
(uly)=) q(uly=) qu p(u)
p(u)
2
(48)
(u y
R
In the following derivations, cases where Eq. 48 does not hold will be implicitly omitted. This is because, according to the assumption of bounded importance weights, the exogenous distribution's 's probability measure corresponding to the cases where Eq. 48 does not hold is zero (in other words, the probability that the inequality Eq. holds is 1). Therefore, omitting the cases where it does not hold does not affect the expectation w.r.t. the exogenous distribution.
For ease of implementation and computation, we focus on log-weights by constructing Ey. (u) to bridge the inequality conclusions over Eu--P [cy, (u)] and the inequality relationship with the variance, where Ey. (u) log (p(u)/q(uly, (no())(u))
En-Po (u))
(49)
p(u) q(uly.).
Eu-Pu
log
(
(50)
p(u)
= =Eu~Pu
log
((Ilx.))
k existence
(51)
p(u)
(ullr=)
p(u)
< EurPu
log
h²2
Eq. 48
(52)
q(u (u)
p(u)
c(ulx.))
p(u)
EunPu
log
x q(u| Y.())
(n)
(53)
22

---

p(u)
I EurPer log
16?
Monotonicity (54)
q(u) (.)
= =- -uEuPP [log (oogq(((Y(()))
Y.(Y)..))
EuwPu [(o)p(u) + 2 2
Linearity
(55)
Next, look for a lower bound for En [5y (u)j. Firstly, according the bounded density ratio condition, for any u E Pv(('') with p(u) > 0 and y. E V., x- V p(u)/q(u ( & K holds. Therefore, based on the concavity of the logarithmic function, the following inequality holds
(assuming A > 1, as the conclusion is trivially true when K =):
p(u)
2x log
p(u) (.)
K ++1 111
log
(56)
qq(ll.))
R22
2k
Due to the boundedness of Eu-Po K5y, (u)] in the logarithmic space (with S log and 2 log k). by employing the monotonicity and linearity of the Lebesgue integral, we obtain:
En-ro (5,((u)
(57)
2x log K
[
p(u)
x² +2+11
2 R2
Eunro
(58)
-
q(uly.)
2x
2k log Ki
p(u)
+1+1 +1
EunPu
[q(uly)) 1 18u1))) (w)
Eurou [Iano:)(u))
(59)
R2
2x
2x log Ki
(1s2 + 1) lllg
Bunro [oy. (u))
(60)
2
R² -1
Substitute the above inequalities into Eq. 5
Vu-Quy. [oy;(1)
(61)
Pu loy. (u)] (V())
Eq. 5 Non-negativity
(62)
< PunPu (u)]
(63)
x²
Bunry [5x. (1)] + x²2 +
Eq. 60
(64)
2x 2xlog
2x
n²2
BuwPu [log quu (u)))
Eq. 55
(655
2k log
162
-1
E BuwPu [log((u)= + 3x2
2x log
k
2x
22
Eunty (lo (u))) +c
Extracting constant terms
(66)
2x log
x² - > 1 for k: > 1
$ - EusPu [log (u (u))] + c.
6
2x log
The constant e is solely dependent on x and Pu. whereas En-Pu (u)) equals the negative entropy of Pu.
Relaxation to concentration inequalities via empirical distribution The bounded importance weights assumption (i.e., there exists x 2 such that 1/x S p(u)/q(u y.) < holds almost surely on the support ().))) may be too strong for general cases, especially when the exogenous distribution has an infinite support set and ().) is also improper. For instance, even for two Gaussian distributions, their density ratio cannot be guaranteed to be bounded (18), E.g. 4.1).
One remedy to mitigate this strong assumption to weaken the statement of the theorem from "almost surely" to "with probability at least 1 - 8." In this way, if 5 is small, the conclusion can still hold with high confidence. To reach such a conclusion, we typically need to use concentration inequalities and introduce other weaker assumptions. One of such weaker assumptions we introduce that for any y- E Ωy the expected value of the square of the importance weights under Pt; is finite, that is, En-ro
(u|y.)] < +oc. Another assumption is about the boundedness of the importance weights in the empirical distribution Po:: this is a discrete relaxation of the original boundedness assumption of importance weights in the exogenous distribution, which makes Thm. 1 trivially hold on the empirical distribution.
23

---

Corollary (Tm.. with Cantelli's Upper Bound). Assume that for any Ve, the expectation Denote by Po the empirical distribution formed by i.i.d. samples u Pu, = 1,...;" n drawn from Pt. and such that 1/K S (/ql))/q(uul ()) y.) for all I,...,n.
Eu-Po [p² [-))y]]]== +00 Let = supy.eR. Eu-Pu ( (u'(u)/4(ulxx-)
assume there exists constant
Then, with probability at least 1 - 6, the following inequality holds:
Vu-Quy. loy, (u)) - -B u,, llog (u (.))
1) 1)= cc
(67)
where the constant e is the same as in Thm.
Proof. First, consider the empirical distribution Pu formed by ii.d. samples u(9)
Pu. Given
the existence of the bound к for these samples, meaning that the importance weights exist on the support of Pu, it is evident that Thm. holdd for P:
Vu-Quiy- Quiy- Joy, (u)] S Lu [log (u Y. +
(68)
where V denotes the empirical variance, which calculated as:
Vu-Quyy Joy, [o)( E u loy; (u)] - (V))
Eq. 3
(69)
- oy (u p2 (2)))
(70)
where we treat P().) as constant.
Next, need establish the inequality relationship between the empirical variance and the true variance. According to Cantelli's ' inequality, we have:
V
[X]
P(X E S -X)
(71)
X X2
Reversing the inequality direction, substituting the random variable X with [[X] = 1 Σm X (where each X; i.i.d),, and letting
V (t 1) yxI obtain:
n
E E + V0-1
V [
P
E
X
11
(72)
Let X; oy oo..((())) (u and subtract P(Y.) from both sides, then have:
p Vu-Qur. loy. (u)] < V Vu-Qur. [oy. + / V PurPo Joy, (u)
(u)]
1 116.. (73)
n
where
V Pu-Pu [u)] = Eu-Po Pu [o]. (u)]
- 1 Pu (ov.
(74)
BurPo Po. [o}_(1) (u)
(75)
p p2(
sup E Eu-Po y.ey.
(76)
q² yz)
M.
(77)
Substituting Eq. 68 into this equation, we obtain the result:
( ((-1) v V Wu-Quis [oy. (u)) < EE -Bu-A [log((u Y.)) + +
P
+c 21-8,
(78)
where the constant the as in Thm.
The following obtains a tighter bound after introducing the constraint of aub-Gaussian conditions (specifically, assuming that the two-sided Bernstein's 's condition is satisfied):

---

Corollary 3 (Thm. with Bernstein's s Upper Bound), Assume that for any y, E V. the expectation Additionally,
En-Po (u)/qq [y)] +oo. Let supy.ey, Eu-Pu (p"(u)/q"(ulx,
assume that the variable oy. (u) satisfies Bernstein's condition, meaning that there exists parameter b > 0 such that for any À E (-1/b, 1/b), the following inequality holds:
(V (VumPy (oy. (w)))1/2))
Py (exp (N(r,, (u) Pu-no (oy. (u)])] S exp
(79)
1 b-b
Denote by Pu the empirical distribution formed by ii.d. samples u(() Pu. = 1..... drawn
from Pu. and assume there exists constant 1 such that 1/x p(u (/(((()))))))))))))))))))))))/ ((. yy. S K for all
i 1,...,. n. Then, with probability at least the following inequality holds:
b V Vu-Quly. loy, (u)] S -B uu [log (-))))
- log (3)
1 L + c, (80)
+
2 log
where the constant e the same as in Thm.
Proof. When Bernstein's 's condition is satisfied, the following inequality holds:
1
1 EE [V]
P E x E
[X]
<
log
2 log
6.
(81)
n
The remaining derivation is similar to that Cor. 2
It is trivially provable that / - / 2
(+) when å approaches zero. This implies that, for the same d close to zero, the remainder term in Cor. 3/s typically smaller than that in Cor. meaning that the former provides a much tighter upper bound with high probability.
guard proposal distribution to ensuring boundedness For a proposal distribution Quy-- define the bounded part as 52( {u 1/x s p(u)/q(u y.) x). If the probabilities of the unbounded part under the exogenous distribution and the proposal distribution, Pu (Ru \ \
Eu-Po [novaly) and Quly. (ou \ ( = Bu-QQuy.
I Du vabs are known, we can splice the exogenous distribution over the unbounded part to form proposal distribution that ensures the bounded importance weight condition:
Proposition 2 (Importance Weight Bound Guard), For proposal distribution Quly. construct new proposal distribution as follows:
1- I-Po(auiio (auiin()))
L-Qui (Du) q q(//),), if s
(uly.)
(82)
p(u),
otherwise
which we refer to as the guard of the proposal distribution Quy, Then, for any u E Qu there exists constant
Pu Lu 1,2(2) N = max
Quly- (Qu \. 1.2(²))
(83)
1 -u-Q. @uu \ s2( Pu au 2040
such that 1/k' (ln)(uly)))// Sn'
25

---

Proof. First, we prove that the constructed proposal distribution is valid distribution, mainly by showing that satisfies the normalization condition:
[ g (uly.) du nuu
(84)
- Pu (se (1(())
[60 1 - Quix. (ou | (.2(()) g(u y.) du S Sauivoby p(u) du
(85)
1 1 Pt (sav \ 1
(1- Qus (Nu ((l)) + Pt (Au (24)
(86)
Quy-- S \ (2(())
=1.
(87)
Next, we prove the boundedness: when u E 2 according to the definition of [2 we have 1/K S p(u)/q(ul S k. Since K' 2 it follows that 1/x' S ppuu)/gh (u y.) S R'. When u & 2.rl trivially holds that p(u) ( (uly.)=11
In the subsequent experiments, we assume that the probability of the unbounded part could be non-zero but approximates zero. Formally, Po (g(b) ~ 1 and Quyy (ol)) 1.This implies that
the error terms in Cor. 2|nd Cor. 3 would be small enough, that when using the guard proposal distribution [Prop. 2] we could obtain a relatively unbiased estimate without the need for re-weighting (where (1 -
Pu Pu Pu \. (2)))) (1-Quy. (Su \ .2(" ~ 1).
The generalization of Thm. Cor. is generalized from Thm. in a straightforward manner through monotonicity. particular, the relaxation and guarantees discussed earlier also apply to Cor. when the scope of assumptions extends from a single counterfactual variable V. to multiple
counterfactual variables of different forms, yy! ,8 E S.
Corollary 1 (Expected Variance Upper Bound). Let Qs denote arbitrary distribution defined over the state space S of a stochastic counterfactual process, and let p(S) denote an arbitrary distribution defined over the o-algebra >(6 corresponding to a set of counterfactual variables Y given a state
E S. Py. is the joint distribution induced by p(a) and Qs. If the conditions in Thm. are met for
any () and any S, then for any QU y. where y. ()(() and any S:
E By. [Vn-Quyy [oy. < -EnQs EumPo [lg qu Y (w)) + c,
E
(7)
where the constant is the same as in Thm.
Proof. For any state s E S, and any y!!) E xo(e) according Thm. if y! satisfies the bounded
density ratio condition, we have:
Vu-Qusy Joy: (uly.)l & Eu [log Y. Y.(u)) + c.
b
According to the monotonicity expectation, we take the expectation w.r.t. Pl on both sides, then:
E By00 p(o)
p(o))
V Va-Qui. (ox. (uly! (u]
(88)
< E 2g00
Euw-Pu [log (u Y. +c]
(89)
E Byla plo) En~Pa olog q Y (u))]l + E E [c]
(90)
BurPu [logq(u| Y! y!! cc.
(91)
where E Eu-Po logg(u| Y!! (u)) is independent of yo(e) and the constant e irrelevant to {)(*)
26

---

We then take the expectation w.r.t. Po on both sides of the inequality, producing:
By,nPs.
Vu-Quiy loy, (u/y.)]]
(92)
Es
E B)(0)>
V [o (u [+>0]]
(93)
x.
I EENQ8
Bu-Pu
loga(u Y (u))
+
(94)
R
Eur
log Y! (m))
+ Rx-Qq [c]
(95)
-
Rgnqs
Bu
log ((u Y! (w))
+
where e is a constant independent of any state s S.
A.3 Learning and Inference
Learning Cor. 1 provides upper bound that
directly gives an optimization objective:
Algorithm Exogenous Matching Learning
[Bu-.r (l x4"(w))]
Input: Exogenous distribution Pu, prior dis- tribution for states Qs- conditional distribu- tion model QUIY_ with parameters oy., mini- batch size and learning rate 17. while Eq. &not converged do for i + and i S do
arg min,
E.-qs
This equivalent to learning the conditional
distribution Qujx, using maximum likelihood. To obtain an estimate of this log-likelihodd, can directly use Monte Carlo methods, which simply involve two sampling processes and one computation: i) sampling state from the prior state distribution Qs; ii) sampling u from the ex- ogenous distribution Pu; iii) computing the log density (or probability for discrete exogenous variables) log quul Y (u)).
1:
2:
3:
8~ Qs. u ~ Pu
4:
0y. - fy. - Vo,, logg (u))
5:
end for
6:
end while
Figure 5: The learning algorithm for Exogenous Matching.
Assuming that we can compute the gradient of g(u y-) w.r.t the parameters fy, of Quly. as modeled in the paper, then combining with the minibatch gradient descent algorithm, we can write the complete algorithm for the training process as Fig.
Inference If for any y. E V. the support of the proposal distribution Quyy covers the sup- port of Pt. then constitutes a valid proposal, and unbiased estimate be obtained using [gq. 4 More can be achieved through multiple importance sampling, any estimate obtained from the proposal Quly. is unbiased for such y- E Ve. Consequently, the expectation over these estimates remains unbiased, i.c.,
Algorithm 2 Multiple Proposals Inference
Input: Prior distribution Qy conditional pro- posal distribution QU sample size n.
1: S -- 0
2: for + and S n
3:
y. ~ Qr. (
(6)
OuIy!
4:
a0) - (p(tu) (y.)) (n))))(
P().) Eu-QQuy. [d))xx][(u]]
(96)
5: S So
=y. [Ruu-Qay. [ox])=][ui]]]
9
6: end for
7: Output: P().) + (Eiro()) /n
which corresponds to the estimator used Fq.. and provides more robust results when an ap- propriate prior QY is chosen. Our inference algorithm involves two sampling processes and one computation: i) sampling y. from the prior distribution Qy ii) sampling from the proposal distribution Quly. iii) computing the importance weights and checking if Y. (u) E y. to obtain (y) This process is detailed in Fig. 6 In our experiments, because the tasks are density estimation (for continuous) and point probability queries (for discrete), we actually still use single proposal distribution at y., which a special case where
Figure 6: The inference algorithm for multiple importance sampling.
Qy is Dirac distribution.
27

---

A.4 Markov Boundary
In order express counterfactuals that span multiple submodels using single model, the twin SCM concept is often used to articulate counterfactuals for a particular setting within a unified framework. Definition 3 (Twin SCM), For an SCM M = (U.V F. Pu) and a set of its submodels M we refer to .M.- (U. F., Pu) as Twin Structural Causal Model (Twin SCM), where V. = = Soul nous variables in the i-th submodel, and fxx denotes the corresponding causal mechanisms.
[Mijkj)..Maaee) Ma)re))... ................................
and F. = U Fipe]- Visx.l represents the endoge-
In general, the twin SCM amalgamates multiple distinct submodels into one SCM, thereby rendering conclusions drawn on single SCM applicable to the twin SCM. It is trivial to demonstrate that if the original SCM is recursive, then the twin SCM remains recursive. Consequently, the causal graph derived from the twin SCM is also an ADMG, and the augmented graph is a DAG.
Our proof is based on d-separation [75 76. which is a criterion for quickly determining conditional independence on a Directed Acyelic Graph (DAG). Suppose a distribution is denoted by P and its G. We say that a path p is blocked by a set of nodes Z if and only if: i) p contains a chain A + B + C or fork A f- B - C' and the intermediate node B is in Z; ii) p contains a collider A + B < C, and neither the intermediate node B nor its descendants are in Z. If all paths between two sets of nodes X and Y are blocked, then X and Y are said to be d-separated given Z, denoted as X le Y Z, implying X Lp Y Z, i.e., X and Y are assume faithfulness, which means that if the
corresponding DAG is
conditionally independent given Z. In this work, conditions for conditional independence are met, then d-separation holds. Combining acyclicity allows for mutual inference between the d-separation and independence.
First, let us restate the definition of the Markov boundary.
Definition 2 (Counterfactual Markov Boundary). For exogenous variable U, U and set of counterfactual variables Y, along with their joint distribution P. If U, is independent of Y \ X given X under P. i.. U, _ (Y)) \ X) X, then X is termed a Markov blanket of U, on Y. The collection of all Markov blankets of Uy on Y is denoted as B,(Y), If X e B, (Y) is a Markov blanket of U, on Y, and there exists no X' ç X such that X' E BY)), then X termed a (counterfactual) Markov boundary of U3 on Y, denoted as (Y).
This definition relies on conditional independence. With the aid of the augmented graph Gd and the d-separation based on the twin SCM, it is possible to derive Thm. 2 and Thm. from the properties on the graph.
Theorem 2 (Counterfactual Markov Boundary Independence). If = U Yijx,) and each Y ifx.] corresponds to a different submodel Mx,, then for each U, E U, there exists a Markov boundary B,(Y.) = U B,(Y, on Y., where By (Y ((x)) is Markov boundary on Mx..
Proof. We aim to demonstrate that U._1 (Y(i).))) forms a Markov boundary for Y., where B,( (()x.)) denotes the Markov boundary over Y [px.).
( U (Y
(Yxx.)) (xx.)) is Markov blanket of Uy: For any submodel Mj, according to the faithfulness assumption, since BY is a Markov blanket over Yix.). it blocks all paths hetween U, and Yssx.]] B, Y (ixx.) in the augmented graph of the twin SCM. We proceed to prove that, given (Y (x.)) does not alter the blocked paths of B, (Yx..)
). For a blocked path p E P. unblocking can occur only at the intermediate node and its descendants in the colliders. However, the nodes in BY (x))) .) are not Y (xsx.): thus demonstrating that its blocking status remains unchanged. Therefore, U (Y.x.)) is a Markov blanket of Uy.
(ii) Any proper subset of U (, (Yxx.)) is not a Markov blanket of Uy: Suppose that a proper subset D C U (Y (xx))) is a Markov blanket. Let D, D n Y (xe be subset in some submodel of this Markov blanket. Obviously, since any proper subset of B, ( (xx)) is not a Markov blanket in Y six::. Dj is not a Markov blanket on Y This implies that there exists a subset S, of Y ((x.) which is not independent of Uj given Dj. We now show that S, remains dependent on U/ given D. By the faithfulness assumption, as S, not independent of U, given Dj. there still exists an unblocked path from S, to U3 in the augmented graph of the twin SCM. Continuing with the condition
28

---

Algorithm 3 Finding Counterfactual Markov Boundary
Input: An exogenous variable Ug E U, a set of counterfactual variables Y., an augmented graph go of SCM M (V, U, F, Pu).
1: B,..)) Y.)+
2: for Y [xx.] C Y. do
3:
( +
4:
for E Y ([x.] do
5:
if D-SSEERATTION (U,), Yxx (x,)9%) then
6:
B (Y t (xxi)) U (Y)
Thm.
7:
end if
8:
end for
9: (Y...) + (Y,Y)) (xx)
Thm. 2
10:
end for Output: counterfactual Markov boundary BY...).
Figure 7: Inferring counterfactual Markov boundaries through d-separation, Thms. and 3 with faithfulness is assumed. D- D-SEPARATION(XX, Y, Z, returns true if and only if x L Y Z.
on D \ D,, we discuss all the paths from Sj to Uj. For an intermediate node of collider on path P, since it is unblocked, its descendants must be conditioned, thus continuing to condition on D Dj does not affect the colliders. For a non-collider on path p E P. for p to be blocked, there must exist an intermediate node of a non-collider in D Dj. however, (D Dj) Y(xx]] 0. so this condition cannot be met. Hence, continuing to condition on D D, does not affect any unblocked paths. According d-separation, given D, Sy and Uy remain dependent. Therefore, D is not Markov blanket, which contradicts the assumption. This means that no proper subset of U (Y.l.) is Markov blanket of Ug-
Theorem 3 (Counterfactual Markov Boundary on Graph). For an exogenous variable U, E U and counterfactual variable set Y, from the submodel Mx. the counterfactual variable Yx (Yx)) and only if Yx &gg U, Y [Yx}. i.e., when given Yx \ (Yx), Yx and Uj are not d-separated on gg, where O% is the augmented graph induced from the submodel Mx.
Proof. (i) Xx E By(Yx) implies Yx M U, Yx ((x): Assuming Yx 11 U, Yx (Yx) according to the faithfulness assumption, all paths between U, and Yx are blocked given Yx \ [Yx)- If a path consists solely of colliders, then, since the intermediate node are given, the path remains unblocked, contradicting the assumption. This implies that on every path, there exists at least one non-collider intermediate node, denoted this set by R. Since each path has least intermediate node, given R, all paths between Yx and U, are blocked. Let R = R n BN (Yx)) and R = R. R For R the Markov boundary property holds, blocking all paths between U, and Yx given (Yx)). Thus, for every path between Yx and Uj. it either passes through nodes in R and is blocked given B (Yx)x) Yx. or it passes through nodes in R We need to prove that the paths that pass through R are also blocked given B, (Yx) x Yx. Taking any R E R., for any path p between R and Uy. if p does pass through Yx. then p blocked given (Y Yx, which means that any path from Y. passing through R and then through p to reach U, is also blocked by (Yx) \ Yx. If p passes through Yx, it must pass through another R' R and eventually reach U" through a path that does not include Yx. which is also blocked by ( xx) Yx. Therefore, any path involving Y is blocked by B, (Yx) ) Yx. which means that any path from U, to Yx is also blocked. For colliders on ancestors of Yx not conditioning on Yx adds blocking to the paths, not removing any existing blocks; for Y₂ being a non-collider, any path involving Y, is blocked by (Y Yx, thus remaining blocked even if Yg is not conditioned on. Therefore, Y2 does not influence the blocking status of any relevant paths. Consequently, given (Yx) \ Xx, all the paths between any other nodes in Yx (including Yx) and U, are blocked. This implies that B, (Yx)) ) However, it is proper subset B, which contradicts B, (Yx) being a Markov boundary. Therefore, Yx e By (X)) implies Yx Kgg Uy Yx (Yx).
non-collider
Yx is also a Markov blanket.
29

---

(ii) Yx La U3 Yx (Xx) implies Yx E By(Yx): It can be proven that the contrapositive holds. If Yx not the Markov boundary, then it is independent U, given the Markov boundary Yx E BY(Yx) x). Clearly, this conclusion follows directly from the definition of the Markov blanket. Therefore, Yx V og U, Yx (V²) holds trivially.
Our final algorithm for efficiently computing Markov boundaries is obtained directly from d- separation, Thms. and 3 as shown in Fig. 7 However, a limitation of these theorems is the reliance on the Taithfulness assumption, which may not always hold, even in the case of fully specified SCM, where numerical violations can potentially occur.
B Related Works
B.1 Importance Sampling
To estimate Eq. 2 a straightforward approach is to start from the definition and employ rejection sampling, as implemented in [T08]:
P().) - Ea-P. [L)o()=)(u]_
(97)
This method operates effectively within a countable, finite exogenous space. However, as discussed in the main text, exhibits low efficiency general settings due the infrequency with which the indicator function in the counterfactual probability expression is effective (i.e., non-zero). One solution to this problem is importance sampling, where selection of the proposal distribution is critical issue. In this section, we discuss related works in detail.
Cross-entropy based importance sampling is straightforward to prove that when the density of the proposal distribution (u)) x f(u)]f(u) (where, in the context of counterfactual estimation, f()) is (D2)) (u)). ). the corresponding importance sampling estimator possesses the minimum variance. This proposal distribution is known as the optimal proposal distribution, denoted as Qū The idea behind a class of methods is to minimize the cross-entropy (or equivalently, the KI. divergence) between the conditional proposal distribution Qu and the optimal proposal distribution Qū- The expression for this in the context of counterfactual estimation is:
Ha 4) -E-ee [log q(u)] x Eu ( log q(u)].
(98)
In the task of sampling rare events, due to the difficulty of directly sampling from Pu to obtain effective samples for which f(u) is valid, it common to estimate the cross-entropy by sampling from a proposal distribution Qu. This method, known as importance sampling, is considered adaptive because continuously adapts based on the samples already drawn from the proposal distribution. For instance, (32) modeled Qu as a Gaussian Mixture Model (GMM) and optimized an equivalent form of cross-entropy to improve sampling efficiency. When applied to counterfactual estimation, the corresponding optimization objective is:
au log s(u) q(u)
arg min -E Eu-Qu
(99)
Neural importance sampling In addition to GMM, normalizing flows provide a more expressive capability as density estimators. Therefore, as enumerated in the main text, numerous recent methods have leveraged normalizing flows as proposal distributions for importance sampling. Many of these methods still rely on optimizing cross-entropy (and KL divergence), as introduced carlier. Another optimization method discussed in [68] is based on the x² divergence between the optimal proposal distribution Qú and the proposal distribution Qu:
Dxa("" (g)))= = /
(u) - q())²² du
(100)
q(u)
Pn
(q)) du
Sa, 2 (2//, g" ('uuduu / q(u)du)
(101)
q(u)
(q" (()) du
(102)
q(u)
Du
30

---

p((u)laaDas( du-1
(103)
q(1)
Ba
p(u)=
Eu-Po
Inooos( - g(u)
(104)
can be seen that the only difference between Eq. 104 and Eq. 5 lies in a constant term. Thus, directly optimizing this equation is equivalent to optimizing the variance itself. Naturally, considering the difficulty of sampling from Ptr as discussed in the paper, an adaptive approach is employed. This involves sampling from the proposal distribution Qu. with the corresponding optimization objective being:
p (u) (u)
arg min Qu~Qu
(105)
In 688 the proposal distribution Qu is modeled as NICE 24
Extension and limitation The above method, without extensions, is only applicable to estimate single counterfactual probability P(V.). . We extend these methods similarly to our proposed approach by introducing a conditional proposal distribution QuIy., and employing the same estimator in Eq. 9 during inference as baselines. During the variance optimization phase, by introducing prior distributions Qs and Pu and letting y. y (u') be sampled from stochastic counterfactual processes (u' ~ Pu and s ~ Qs), we reformulate the optimization expressions for cross-entropy based importance sampling (CEIS) as:
p(u)) (u) (uly~)
log log (ol ]]
arg min Ey, Bu-Quix
(106)
and neural importance sampling (NIS) as:
p²(u) (uly.. 10000>)((1)
arg min By. Ba-Quiy.
(107)
where 8(y.) denotes a small region around y.. The above optimization objective tasks Quly- with sampling from (r where any V. can be decomposed into of y! Y., such that V.- - ((y("), ()), ). and (V+) = S0(66(
An unavoidable drawback of the extended method is the presence of indicator functions in the optimization term. the initial proposal distribution Qu is similar to Pa in performance, the entire optimization process requires a prolonged warm-up period before effective samples start to emerge from the proposal distribution Qu. For example, if the entire space is the unit plane [0, 1|2 and (o(() is a small region {(x, y) x 0.001, < 0.001). then if the chosen initial proposal distribution is uniform, it has only probability 10-6 of sampling an effective sample, resulting in low efficiency for optimizing the loss functions.
In contrast, the optimization term in Exogenous Matching according to Thm. 1 does not include the indicator function, ensuring that every sample is effective during the optimization process, thus significantly enhancing learning efficiency.
B.2
Proxy SCMs
In this section, we provide a detailed description of the construction of the proxy SCMs used in the experiments and their identifiability.
Nerual Causal Model The Neural Causal Model (NCM, 107 108 ) is a recently proposed general approach to mimicking the original SCMs to construct proxy SCMs. For identifiability, its primary assumption is the causal graph of the original SCM. By replacing the causal mechanisms in the original SCM with neural networks that have strong expressive power, it can effectively proxy the original SCM.
Specifically, an NCM is 4-tuple (U, V, F.F, which is almost identical to the SCM definition, except that each fv, E F is a neural network with sufficient expressive power, and Py any well-defined distribution over U.
31

---

Their method for inducing submodels is the same as for a regular SCM, i.e., replacing causal mechanisms under perfect interventions, which they term neural mutilation. For counterfactual estimation, they use the rejection sampling introduced earlier because the endogenous variables in their experimental scenarios are discrete and finite, making rejection sampling efficient enough to sample effective samples.
Regarding identifiability, this work assumes causal graphs and has been shown to operate effectively when endogenous variables are discrete and finite. Specifically, they construct a class of models called G-constrained NCMs, which have the same causal graph 9 as the original SCM. Then, they prove that any NCM in the G-constrained NCM family has dual identifiability with the causal graph G. i.e., the NCMs in the G-constrained NCM family (given some observational interventional distributions) output the same result when answering a Q E C3, if and only if the graph algorithms for counterfactual identification output identifiable results given input Q and observational or interventional distributions.
Furthermore, based on this duality, they develop a sound and complete identification algorithm, and since this identification algorithm is based on the learning process of NCMs, it only requires gradient optimization. Certainly, although this work guarantees thcoretical results for identifiability in discrete and finite settings [108] well partial identifiability in Lipschitz continuous settings 93 further improvements are necessary to apply its estimation methods in broader scenarios.
Causal Normalizing Flow method to construct proxies of SCMs based on flow models. For identifiability, it assumes the causal graph or causal ordering of the true SCM, and the true SCM is Markovian and diffeomorphic.
Causal Normalizing Flow (CausaINF, 42) is recently summarized
Specifically, a CausalNF is 4-tuple (U, (T-1,P)) where T = [Vv, (w:; pay,) I V/2 E v) is sequence-preservinn Autoregressive Normalizing Flow (ANF), and Pu is the base distribution of this flow model.
Their approach derive submodels differs from conventional SCMs because the modeled CausaINTs have multiple architectures, and substituting causal mechanisms is only applicable in the recursive case. Their solution is alter exogenous distributions and leverage the Markov assumption to ensure that the density of the intervened exogenous distribution is:
(u) = II VEX
5(fm: = fvi (u;; pav))))
II p(u;). V
(108)
ó denotes the Dirac distribution, where a density only exists when the condition inside its parentheses is satisfied. In practice, the distribution of the intervention is determined by the following steps:
= (u): iii) Replace VI, resulting in v2 {wr E VI, 02 E x} U x; iv) Compute u2 = T(v2): v) Replace U2, yielding U3 {ur u, E u2 vs & X {m: M2 E u1, 2/2 & x): vi) Compute vx - Here, step (v) corresponds to replacing relevant parts of the density of the intervened exogenous distribution with the Dirac distribution, as shown in Eq. 108 The last four steps among these six steps are employed for counterfactual reasoning, but this form of counterfactual reasoning is not for generalized cases. Specifically, the counterfactual query corresponding to these four steps is P(Vx x V1). Although the submodel cannot be directly obtained, the Vx obtained from the above steps is indeed the potential response Vx(u), so method is still applicable to this model.
i) Sample from Pe:: ii) Compute
As for identifiability, this work requires numerous assumptions, including Markovianity, diffeomor- phism, causal order, and uses results from [106] to assist in proving the identifiability of the learned model up invertible transformations.
B.3 Arbitrary Conditioning
There is a challenge of modeling an exponentially large number of conditional distributions with single model. Several existing works have focused on this problem, which is termed arbitrary conditioning. Arbitrary conditioning is an unsupervised learning task that, given a dataset over V, aims answer any query of the form P(X Y = y), where Y C V and y By are arbitrary.
The initial step to address the problem of arbitrary conditioning is to encode the condition y into a fixed-length vector. Subsequently, this vector is Iearned through various modeling methods, and the
32

---

training process incorporates various random masks to simulate different forms of Y. Specifically, related works introduce a binary mask b to indicate whether a variable has been observed. In this context, the conditional probability can be expressed as P(VI-bl Vb, b). Consequently, the learning objective can equivalently be viewed as maximizing the following log-likelihone term:
arg max Bb-Pa Ps [E [log q Vb. b)ll.
(109)
where Po represents the prior distribution of the random mask, while Py denotes the observational distribution. The model corresponds to the probabilistic model employed for this task. Specific methods employed in related works are enumerated the main text, involving a variety of generative models that are suited to different task scenarios. For instance, some methods permit direct training with partially observed data, some can only output samples without density estimation, and others are incapable of performing efficient marginal probability computation.
The similarity between this work and these methods lies in the encoding scheme and the training objectives. For example, the conditional encoding method in this paper also employs binary masks, introducing both the observation mask (X) and the intervention mask w wY). Furthermore, the training objective, as defined in Eq. 8 closely resembles [Eq. 109
The distinction between our work and these methods lies in the available data and the target tasks. The problem of arbitrary conditioning is based on learning from the observed distribution and is typically applied to imputation tasks. In contrast, our conditional distribution modeling can be seen as an extension of arbitrary conditioning, generalizing the observed distribution counterfactual distributions. This allows the conditions to originate from multiple hypothetical worlds and is applied in counterfactual reasoning.
The work most similar to ours in terms of overall structure and optimization objectives is Posterior Matching [91] It addresses the problem of arbitrary conditioning by matching the encoder of a given VAE. The optimization formulation is as follows:
arg min Eb-Pa
[Ev~Pv [Bu~Poiv q(ogq(u \vb b)]]]
(110)
where Pujv represents the posterior distribution in a VAE, also referred to as the encoder. By comparing the optimization objective of our method, Eq. 8 with Eq. 110 we observe a notable similarity between the two. We can establish a one-to-one correspondence between the terms in the optimization formulations: the mask b can be analogized to the state in the stochastic counterfactual process; the decoder of the VAE (mapping from latent variables to observable variables) corresponds to the SCM's potential responses (mapping from exogenous to endogenous variables); the encoder of the VAE (mapping from observable variables to latent variables) corresponds to the abduction process in the SCM (mapping from endogenous to exogenous variables). Consequently, the posterior distribution aligns with the abduction process in the SCM, represented by the product of the indicator function and the exogenous distribution.
An intuitive explanation is that posterior matching attempts align the conditional distribution with the posterior distribution, whereas our work can be interpreted as attempting to align the conditional distribution with the optimal proposal distribution. The interpretation of the optimal proposal distribution is detailed in App. B. where in the counterfactual estimation task, the optimal proposal distribution is proportional to the product of the indicator function and the exogenous distribution.
B.4 Injecting Prior Knowledge
This work involves injecting prior structure into the learning process. Specifically, such methods of injecting prior structure are achieved by directly influencing the internal connectivity of the model.
Related work includes MADE (3), Zuko [8]), and StrNN [3.. MADE interprets from the perspec- tive of conditional probability, altering the connectivity between inputs and outputs in autoregressive models by masking. thereby enabling the model to output specific conditional distributions p(z; yy) only when there exists connectivity from inputs in y to output z;; however, this method only re- spects autoregressive properties. The method in Zuko is primarily used for implementing masked autoregressive normalizing flows, but is not limited to autoregressive dependencies. Based on the interpretation of Jacobian matrices, it alters the connectivity between inputs and outputs in ML.s through masking, ensuring that the Jacobian matrix of MLP outputs is non-zero only when the
33

---

mask is one, thus representing the dependency between input and output variables. However, the masking algorithm in Zuko only approximately preserves the maximum number of connections under dependency constraints in variable order. StrNN extensively investigates the issue of variable dependency and connectivity, formalizing it as an integer programming problem, and providing exact and greedy algorithms that outperform the previous two methods in terms of performance.
The tasks typically addressed by the above methods include representing autoregression (autoregres- sive models), constraining the dependency of Bayesian networks (Bayesian network inference), and constraining the dependency of causal graphs (causal inference). Different from these tasks, this paper applies them to representing the dependency of Markov boundary (Markov boundary feature selection). The masking algorithm adopted in this paper is Zuko, as it can be directly combined with the normalizing flows it implements and used directly in subsequent experiments.
Experiments
C.1 Metries
We formally introduce all the metrics used in the subsequent experiments, which aim to reflect the quality of counterfactual estimation:
Effective Sample Proportion (ESP) As mentioned earlier, effective samples here refer to samples that make the indicator function equal to one (i.e., the constraints of potential responses are satisfied). For specific V.. we define
7))) 1 Σ (x)
(111)
which reflects the proportion of effective samples sampled in counterfactual estimation for V. In rare event sampling, this metric essentially represents the success rate, indicating the proportion of samples in which the rare event occurs. ESP is about the expectation this function over prior distribution Py, (related to the stochastic counterfactual process):
ESP Ey.~Py. [n()).))
(112)
Specifically, ESP reflects the sampling efficiency when the importance weights are ignored. It focuses on the probability that the sampled u falls within the support Pu ()².), for all V
Failure Rate (FR) This work claims that the proposed method holds the potential for counterfactual estimation in any form following single training session. To verify the reliability of this claim, we introduce the FR metric, which reflects the efficacy of the proposal distribution post single training across the entire state space of the relevant stochastic counterfactual processes, defined as follows:
( m
FR -y.-Ps.
(113)
otherwise
where < m < is a constant representing the acceptance threshold. An v)))) below this threshold is considered failure. The closer the FR is to zero, the more it indicates that for the vast majority of D.. there are sufficiently many samples of u falling within the support Lu,,,, making the counterfactual estimation considered tractable.
C.2 Design of Stochastic Counterfactual Processes
This paper introduces stochastic counterfactual process Def. 1. allowing the proposal distribution after a single training to potentially estimate counterfactual probabilities multiple or even infinite forms in tractable manner. This facilitates efficient estimation of a L3 expression in real-worle scenarios (assuming that the L3 expression can be approximated by a finite set of counterfactual probabilities). Therefore, if the proposal distribution after single training can tractably estimate these counterfactual probabilities, it implies that estimation of the La expression only requires one training. Specifically, employing stochastic counterfactual processes, this paper conducts stochastic optimization on various forms of counterfactuals and demonstrates its feasibility (Cor. Tp. In this section, we elaborate on its specific design in experiments.
34

---

Table The state space of 2 (i.e., all involved counterfactual variables), as well as the expression of Q. Q E (ATE, ETT, NDE, CtfDE).
Counterfactual Variables (Y.)
Expression
Yx=1 Yx=0
ATE
P(Yx=1=1) = P( = 1)
xXx=1,XX X (Y X}
ETT
( P(X-1)
{Xx=1. ox w E Sw Yx=o
NDE
En = -
(Yx=1,== Wx=0h w E Sw Yx=0 X
Ew P(Yx=I,Ew 1 P(Yx
CtDE
P(X=1)
Construction of ~ If Q E L3 is an arbitrary P(V.) such that V, S Sty is ensured, and precisely involves & distinct submodels, then can design D answer any such Q. Specifically, the state space of _ S = {s | = k}.
Here, design the prior distribution Q to cover S:
blu ~ Bernoulli(p)),
0< <1,
(114)
b00h = [b,} V2 E V},
(115)
x (V, b1] 1],
(116)
(1)
Pxi,,
p(x(1) 0,
(117)
2
Bernoulli(p2),
0 0< P2 <1
(118)
(0)
[b (1)
V, E V V x
(119)
2j
y
(V, (1) = 1},
(120)
Y. = (Y(0)
(1) x(0
i S k},
(121)
where Bernoullip:) is a Bernoulli distribution with parameter P1, meaning has a probability P1 yield 1; Px the endogenous distribution corresponding X and for any fxx p(x) 0.
This distribution, due to the constrained Bernuulii distribution parameters (0 < P1. 1), implies that any X, Y C V may potentially be sampled; and for intervention values, assuming positivity (p(x(()= 0), any possible intervention value can be sampled. Therefore, for any of these submodels, any counterfactual variable is covered, and hence any variable set is also covered.
In subsequent experiments, we set P1 0.2 and P2 0.75, which are utilized during both the training and testing phases. The rationale behind employing 2 for experiments lies in the large size of S, which provides more credible evidence to verify whether the conclusion of Cor. ] satisfies the claim and whether the designs bridge the gap between theory and practice.
Construction 2 apply our method real-world tasks, we construct stochastic counterfac- tual processes 2 for several common causal query tasks, where Q E (ATE, ETT, NDE, CUDDF).
As shown in Tab. 3 the state spaces for these queries are finite or easily sampled, allowing the design of the prior distribution Q to simple uniform distribution over these available states. We assume that the endogenous variable space is discrete and finite, which means that we only need to estimate counterfactual probability for each of them. By computing according to their expressions, we can obtain the estimate for Q
35

---

C.3 spely-specified SCMs
In this section, we briefly introduce all fully specified SCMs used in our experiments.
Markovian diffeomorphic SCMs These SCMs are sourced from 42 and are primarily used in the experiments for CausaINF, including: CHAIN-LIN-3, CHAIN-NLIN-3, CHAIN-LIN-4, CHAIN-LIN-S, COLLIDER-LIN FORK-LIN, FORK-NLIN, LARGEBD-NLIN, SIMPSON-NLIN, SIMPSON-SYMER
TRIANGLE-LIN and TRIANGLE-NLIN, Markovianity denotes the ab- sence of hidden confounders in these SCMs, meaning that each exogenous variable is involved in the causal mechanisms of only one endogenous variable. The diffeomorphism a parametric assumption, positing that the transformations between exogenous and endogenous variables satisfy the properties of a diffeomorphism. For detailed construction of these SCMs, see the Appendix of 42
Semi-Markovian continuous SCMs These SCMs are described in [07]. We have reformulated their functional forms to adapt them to continuous scenarios. Specifically, this involves replacing their causal mechanisms with combinations of linear and nonlinear functions, all of which originate from the Markovian diffeomorphic SCMs already present above. These SCMs include: BACK- DOOR, FRONT-DOOR, M, and NAPKIN. Here, Semi-Markovianiie indicates that there may be latent confounders within these SCMs, yet they still satisfy recursiveness. Continuity refers to the assumptions regarding their parameters.
Regional canonical SCMs These SCMs are originated from [008. where all endogenous vari- ables are discrete and finite (binary). with complexity arises from the causal mechanisms being type of stochastic mapping. These SCMs include: FAIRNESS, FAIRNESS- XW,, FAIRNESS-XY FAIRNESS- YW, where subscripts indicate the presence of hidden confounders between variables. Please refer to the Appendix of [108] for the specific constructions of these SCMs.
C.4 Reproducibility
In this section, present details pertinent to reproducibility, encompassing the complete architecture of the models utilized, the computational resources employed, and the exccution and random seeds of the reported experiments.
The entire model architecture involves three components: the sampler, the conditioning model, and the density estimation model.
Architecture details
The sampler is used to sample from the prior distribution Qs of the stochastic counterfactual process and u from the exogenous distribution Pt. then compute the potential outcome { (u). In our experiments, we selected a batch size of 256 and generated training set of size 16.384 through the sampler. The validation set is also generated using this sampler. Specifically, we utilize the same sampler to construct Py.. from which random y. is sampled. For discrete distributions, counterfactual events are assigned as V. {y.). ). and for continuous distributions, counterfactual events are assigned as ), E (y.), where dc is cube centered at y. with side length 1.
brief overview the conditioning model is shown in Fig. 2] where both g and h are neural networks. When masking with Markov boundaries, h is a neural network that allows dynamic influence on internal connectivity, as depicted in App. B.4] On the other hand, g is permutation invariant neural network used to aggregate hidden encodings from multiple different submodels. We conduct ablation studies related to the construction of these two neural networks in App. C.8
For the density estimation model, we chose normalizing flows and GMM. Specific models for normalizing flows include MAF (73,, NSF (26.. NCSF [80]. NICE (24) NAF [39), and SOSPF 41. The implementations of these models are based on the Zuko library 8I, with some internal
implementation details modified to fit our work. Regarding the hyperparameters of these models, the number of components in GMM is fixed at 10, the number of transformations for all flow models set to 5, and all the neural networks involved contain 2 hidden layers with 64 neurons each (in some experiments, 256 neurons were used, which will be specifically indicated and discussed subsequent sections).
36

---

Hardware Each experiment was conducted on GPUs, and the primary computational overhead arising from the sampling process. All models, including the pre-trained SCM agents, were trained and tested on an NVIDIA RTX 4090 and Intel(R) Xeon(R) Gold 6430 platform. The training time for EXOM primarily influenced by the number of variables the SCM, the number of submodels, and the chosen density estimation model. For example, in the case of SIMPSON-NLIN with |s] 5, one epoch (16,384 samples, batch size of 256) including the time for the sampler to prepare data takes approximately 8 seconds. In all EXOM experimental setups, we fixed the training to 200 epochs. All experiments, including the unreported ones, consumed approximately 1000 GPU hours total.
Execution To reflect the error, each type of experiment was conducted five times with different random seeds, and the range of errors exhibited reflected in the experimental results. During training, EXOM does not require sampling and is trained for a fixed 200 epochs. The sampling size for CEIS is set to 10 and trained for 200 epochs, whereas NIS has a sampling size of 10 and is trained for 50 epochs (we found that its training speed is excessively slow due to the time-consuming nature of sampling from the flow model). During validation and testing, the validation set queries 1024 counterfactual probabilities from the sampler, with sampling sizes for EXOM, CEIS, and NIS set to 10 and RS set to 10 For all EXOM, CHIS, and NIS training, we used Adam [62] with an initial learning rate of 0.001 the optimizer and ReduceOnPPlaeauu with patience of and factor of 0.5 as the learning rate scheduler. To improve training robustness, the conditional inputs y. and both the input and the output exogenous samples probability densities and importance weights were processed in logarithmic form, with the terms where the indicator function is zero replaced by inf, and the log-sum-exp technique used to estimate the probabilities.
standardized. During inference, the
C.5 Convergence
Here, we present additional results on convergence experiments. Our experiments span across combinations of 4 different types of SCMs (SIMPSON-NIIIN, LARGEBD-NLIN, FAIRNESS-XW NAPKIN) and 4 different density estimation models (GMM, MAF, NICE, SOSPF). All models were trained and validated based on the stochastic counterfactual process D which provides a sufficiently large state space increase difficulty, where |s| 3.
As shown in Fig. 8, for all combinations, LL looellikl,, which is the negative form of the training objective Eq. 8) minimizing [Eq. is equivalent to maximizing LL) increases, ESP improves. In particular, there a clear inflection point in the experiments with three continuous SCMs, after which the ESP rapidly improves until it stops increasing because of the convergence of LL. This series of experiments indicates that can indeed optimize the upper bound (indirectly equivalent to LL) to optimize the variance (as discussed in App.. C.I. ESP and variance are strongly related). Thus, empirically, these experiments support the correctness of Thm. and also empirically indicate that our design potentially satisfies the bounded importance weight condition discussed in Appp. A.2
Fig. 9 reflects the relationship between LL and FR. It is observed that for all combinations, LL increases, FR decreases (with one exception in the SOSPF and LARGEBD-NLIN combination) until it stops increasing due to LL convergence. This phenomenon empirically supports the conclusion in
One Run GMM
Initial State
Final State
GMM
MAF
NICE
SUSPF
NICE
SOSPF
SIMPSON-NIIN
B. LARGERD-NI IN
FAIRNESS-XW
NAPKIN
Figure 8: The relationship between ESP and LL on 4 SCMs with distinct types and 4 different density estimation models. Within these combinations, ESP always increases with the improvement of LL.
37

---

One Run
Initial State
Final State
GMM
MAP
NICE
SOSPF
GMM
MAP
NICE
SIMPSONNNLIN
FAIRNESS-XW
NAPKIN
Figure 9: The relationship between FR and LL on 4 SCMs with distinct types and 4 different density estimation models. Within these combinations, FR always decreases with the improvement of LL.
Cor. and manifests the claim of simultaneously optimizing the variance of counterfactual estimators across the entire state space.
C.6 Comparison
We provide more comprehensive and detailed information regarding the comparison experiment in the main text. These comparisons are conducted on three different types of SCMs (SIMPSON-NLIN, NAPKIN, FAIRNESS-XW) EXOM employs two density estimation models (MAF, GMM), while NIS utilizes NICE, and CEIS uses GMM. All models are trained and validated based on a stochastic 1,3,5. As illustrated in Tab.1 overall, our proposed method, EXOM, outperforms other importance sampling methods, while CEIS and NIS show comparable performance in the aforementioned experiments. The advantage of EXOM is particularly pronounced when dealing with multiple distinct submodels, significantly surpassing the other methods. All three importance sampling methods outperform RS, demonstrating that variance optimization indeed enhances estimation efficiency.
counterfactual process 2 where
C.7 Combinaitons of SCMs and Density Estimation Models
In this section, we detail the performance of EXOM when using various density estimation models (GMM, MAF, NSF, NCSF, NICE, NAF and SOSPF) in different SCMs enumerated in App. C.3 These experiments are based on D with [s] 3.
In Figs. 10a and Tla deeper colors indicate a higher number of effective samples, which implies greater sampling efficiency. For density estimation models, it is evident that MAF and NSF generally perform well, while other density estimation models are suitable for specific scenarios. Regarding SCM, observe that FORK-LIN and TRIANGLE-NIAN perform worse than other SCMs, whereas some SCMs such as CHAIN-LIN-S, LARGEBD-NLIN, and NAPKIN show compatibility only with certain density estimation models. The performance improved significantly after increasing the hidden layer width from 64 to 256.
In Figs. 10h and TTb lighter colors indicate fewer failure cases, meaning that the method works effectively within the state space 2 It can be noted that for any of the aforementioned SCMs. there is always density estimation model that matches and keeps the FR at low level. Further- more, compared with Figs. 10a and Tla there is a strong negative correlation between FR and ESP, consistent with the theoretical results. In regions where FR is particularly high, in addition to the poor performance of GMM due to its limited expressive capacity, there are some notable aspects. For example, FORK-NLIN and TRIANGLE-NLIN are the most challenging settings to learn. Furthermore, as diameter of the SCM increases, the performance EXOM gradually decreases, which is reflected in the performance under the CHAIN-LIN-N series settings. Finally, comparing Figs. 10h and Tbb. some models with high parameter number requirements (such as NSF, NCSF, NAP, SOSPF) show significant improvements when the width of the hidden layers is increased.
38

---

CHAINIIN
0.05
0:08
CHAINUIING INS.
0.79
0.64
CHAINNLINN
0.04
CHAINLINA
0.08
0.71
0.13
0.07
0.57
0.56
0.68
CHAINEIN-
0.80
0.13
0.52
0.59
0.12
0.12
0.12
CRAININNS
0.00
0.48
0.00
0.00
0.57
0.05
0.03
CHAIN TIN
0.99
0.21
0.99
0.99
0.21
0.56
0.50
0.87
0.87
0.81
0.57
0.8T
COLLIDER.LI
0.16
0.00
0.12
0.02
0.02
PORKLIN
0.99
0.14
0.07
PORKLIN
0.43
0.05
0:06
0.04
FORK NENN
0,00
0.17
0
0.01
0.15
0.11
0.13
FORK NT IN
0.96
0,23
0.87
0.95
0.22
0.22
0.21
LARGEND SUIN
0.04
051
0.00
0.00
0::49
B.00
0000
LAROEROANINN
0.61
0.12
1:00
1.00
0.13
1.00
1.00
0.21
0.20
SIMPSON SYAPPOD
0.25
STUPSON SYMPAON
0.10
0.30
THANGLE-LIN-
0.00
0.56
0.59
0.36
0.34
0:57
060
TRAMKELLLIN
0.91
0.11
0.00
0.15
0.31
0:09
0.09
TRIANFLE NLIN.
0.00
0.11
0.14
0.05
0,08
0.09
0,10
TESANGLE
0.98
0.44
0.46
0.55
0.52
0:44
0:44
0.00
0.74
0.05
0.56
0.36
0.05
0.36
0.26
0.67
0.40
0,00
0.17
0.04
0:00
0.03
M. 0,01
0.49
0.33
0.15
0.32
0.30
0.42
M.
0.42
0,06
0.07
0.10
0.10
0:05
0,00
NAPKIN
0.01
0.40
0.13
0.02
0.27
0.09
0.10
NAPKIN
0.61
0.38
00:66
0.67
0.25
00:44
0:41
PARRNESS
PARMESSS
0:01
PARNEESS XW.
0.26
0.29
0.25
FAIRNESS NW.
0,05
0.04
0.04
0.05
0.05
0.06
PAINESS- XY
0,26
0.25
0.41
0.41
0.24
0.31
0.30
PAIRNESS XV.
0.03
0,03
01.4
0.04
0.03
0:05
0.05
PARNESS-YW-
0.18
0.17
0.29
0.29
0.17
0.21
0.21
PARCNESS JW-
0.04
0.04
0.03
0.03
0.03
0.06
0.09
WAF
WESF
NAP
(a) ESP
(b) FR
Figure 10: Performance on 20 SCMs and density estimation models using EXOM. The experiment is based on the stochastic counterfactual process 2) where |s| = 3. The width of the hidden layer is set to 64.
0:00
CHAIN IN
0.46
0.84
CHAIN NLIK
0.02
0.03
0,01
0.02
CHAINGIIY
0.41
0.21
0.24
0.59
0.53
0.72
0.21
CHAIN-IN-
0.30
0.07
0.04
0.12
0.11
0:04
0.07
SI.
0,00
0:5
0.51
0.25
0.52
0.57
CHAINAIN
0.95
0.14
0.13
0.21
0.69
0.12
0.12
0.95
0,90
0.75
0.03
0.00
LIR
0.84
0:85
26
0.84
PORR-IN
0.17
0.00
0.05
0.05
0:00
0.03
YORK NLNN
0,02
0.36
0.37
0.25
0.28
0.28
0.90
PORK MLNN
0.65
0,13
0.11
0.15
0.18
0.12
0.12
ARGERD SUIN
0.02
0.64
0.73
0.55
0.56
0.64
LARGIRD
0.74
0,08
0.02
0.07
0.09
0.04
0.05
0:00
SIMPSON SYMPROD
0.76
0.80
0.67
STUPRON SYMPROD
0,41
0.06
0.05
0:00
0.05
TRINNGE
0.22
0.75
0.75
0.69
0.28
0.75
0.77
TREANGE EN
0.57
0,05
0.03
0.07
0.36
0:04
0.04
TRIANGLE MLNN.
0.01
0.24
0.34
0.21
0.29
TEDANGLE-ANI
0.2
0.19
0.14
0.29
0.53
0.16
0.15
0.59
0.83
0.12
0:01
FRONEDOOTR
0.78
8:84
0:23
088
0.8
FRONT ROOK
0.20
0.01
0.00
0.02
0.01
0.00
0.00
M. 0.04
0.64
0.63
0.49
0:48
0.51
0.63
M.
0.19
0.00
0.01
0.03
0.03
0:01
0.02
NAPOIN
0.02
0.58
0.59
0.40
0.43
0.47
0.53
NAPKIN
0.52
0.30
0.05
0.14
0.17
0:08
0.10
FAIRNESS
0:01
PARINESS XW.
0.59
0.59
0.30
0.40
RAINNESS NWV
0.03
0,00
0.00
0.04
0:01
0.00
PAIRNESS XY
0.27
0.29
0.38
0.60
0,27
0.39
0.41
RAIRNTSS NV
0.02
0,02
0.00
0.00
0.02
0:00
0.00
PARNESS-YW.
0.21
0.53
0.53
0:20
0.35
FARNNESS TW:
0.02
0.00
0.00
0.00
0.01
0:00
0.00
NCSF
GMN
MÀ7
NESP
NiC
(a) ESP
(b) PR
Figure 11: (Fig. 10 continued). Performance on 20 SCMs and 7 density estimation models using EXOM. The experiment is based on the stochastic counterfactual process 2) where |s| = 3. The width of the hidden layer is set to 256.
The impact of hidden layer width on performance Experiments show that increasing the width of hidden layers significantly enhances performance. This improvement is not only due to the increased expressive power of the neural network itself but also because the experiments mentioned above utilize the Markov boundary masking trick. According to the description in App. B.4 we directly inject the Markov boundary masking into the weight matrix, masking certain weights (i.e., specific neurons are always inactive for specific tasks). This implies that when the graph structure of SCM is sparse, the Markov boundary is typically smaller, which leads to a substantial increase in the number of masked weights, thereby significantly reducing the neural network's expressive power. Therefore, increasing the width of hidden layers helps mitigate the degradation in expressive power caused by injecting the Markov boundary.
In practice, the recommended width is closely related to the number of parameters required by the conditional distribution models. For instance, models such MAF, NICE, and GMM require fewer
39

---

parameters, so setting the hidden layer width to 64 is sufficient. In contrast, models like NSF, NCSF, NAF, and SOSPF, which require larger number of parameters, will need a wider hidden layer 256 is recommended since, with a smaller Markov boundary, each parameter connection's is hidden neurons are insufficient. Optimal width configuration should be fine-tuned according to the settings.
C.8 Ablation
Our ablation experiments were conducted on the architecture of the two neural networks g and h in the conditioning model.
Markoy boundary masking Let us first introduce how the mechanism of Markov boundaries oper- ates specific density estimation models. Firstly, we calculate the Markov boundary corresponding to each U, according to Theorem 3 and generate masks m, using the method described later. The subsequent processing is then model-specific.
1. For GMM, h primarily outputs each component's 's mean, diagonal matrix, and anti anti-diagonal matrix. Here, we apply the mask m; only to the mean corresponding to the i-th exogenous variable in each component, ensuring that if the counterfactual variable is not the Markov boundary, the inference of the mean is independent of it.
For Autoregressive models (including MAF, NSF, NCSF, UNF, and SOSPF), each layer trans- formation follows an autoregressive sequence, where all previous layers' hidden encodings may influence the subsequent layer. Consequently, each U.'s hidden encoding appears in each layer. Therefore, we only need to directly apply the mask to the inference of hidden encodings, ensuring that the inferred next layer hidden encodings are only relevant to the counterfactual variables when they are in the Markov boundary.
3. For Coupling models (including NICE), each layer transformation involves only a subset of hidden encodings from the previous layer; hence, each layer's output only corresponds to some U.. Thus, we create masks only for the U, that serve as the output for each layer transformation, ensuring that the inferred hidden encodings are only relevant to the counterfactual variables in the Markov boundary.
Of course, the above initial design will raise some issues in practice:
There are some inevitable discrepancies between theory and practice, which we have summa- rized in Appp D.3
Another issue worth discussing is how to handle the counterfactual Markov boundaries in submodels. In perfect interventions, all connections between intervened variables and their parent variables are cut (All Cut), which also a way to represent causal graphs in submodels. However, theoretically, directly cutting connections would result in intervened variables never participating in the conditional reasoning of their parent exogenous variables. Hence, we have also devised two relaxed and compromised methods: one is to retain all connections (No Cut), pretending that the submodel has not been intervened upon, thus not reflecting the intervention through the mask of the Markov boundary but instead relying on the intervention information included in the conditional encoding; the second to only retain connections with parent exogenous variables (Endo. Cut), while relationships between endogenous variables still cut.
3. As discussed in App. C.7] the introduction of the Markov boundary theoretically increases the model's focus on specific variables. However, this comes at the cost of expressiveness, especially in models with a larger number of parameters. Therefore, in the ablation experiments, we use hidden layer width of 256 for NICE and SOSPF, while for GMM and MAF, the hidden layer width is set to 64. This does not affect the faimess of the ablation, as the goal is to compare the performance difference with and without the Markov boundary. Under the same hidden layer width, the presence of a Markov boundary actually means fewer neurons are involved in reasoning.
Fig. 12) and Fig. 13 respectivell illustrate the effects of Markov boundary masks on ESP and FR under different SCMs and density estimation models, where error bars and outliers are computed based on the quartile algorithm. The first phenomenon that can be observed is that introducing Markov boundary masks Markovian SCMs (SIMPSON-
LARGEBD-NLIN, TRIANGLE- NLIN) results in significant performance improvement compared to not introducing them. There one exception, namely LARGEBD-NLIN and SIMPSON-NLIN, which will be discussed in App. D.3 When applied to semi-Markovian SCMs (M, NAPKIN), the scenario changes, with the improvements becoming less pronounced and the effects post-application not being very stable. Across all intervention designs (All Cut, No Cut, Endo. Cut), these methods perform similarly
40

---

withoua Markey Boundary Masked
Markes Masced (Endis
Markon Boundary Masked (All Cuo
Marker Bounday Masked (No Cuil
ESP
GMM
MAF
NICE
SUSPP
GMM
MAP
NICE
SOSPE
SIMPSON.NLIN
LARGEBD-NLIN
TRIANGLE
M
NAPKIN
Figure 12: Ablation study for Markov boundaries on 4 different settings of SCMs: (a) SIMPSON- NLIN, (b) LARGEBD-NLIN, (c) M, (d) NAPKIN. Higher ESP indicates higher sampling efficiency.
Wilhout Markey Boundary Masked Boundary Masked Coo
Markes Bozrdary Masked (Eadis. Cull
Marken Boundary Masied
MAP
SOSPE
GMM
SIMPSON-NLIN Mes
LARGEED-NNIN
TRIANGL FNJIN
143
M
NAPKIN
Figure 13: Ablation study for Markov boundaries on 4 different settings of SCMs: (a) SIMPSON-
NLIN, (b) LARGEBD-NLIN, (c) M, (d) NAPKIN. Lower FR indicates broader coverage.
without a clear outperforming approach. However, we find that when the All Cut method performs weaker than that without Markov boundary masks, the latter two serve as compromises and indeed provide better results than All Cut. Given these observations, Endo. Cut was chosen in our previous experiments.
Aggregation function The next problem to address is the choice g. We have posited the use of cormutation-invaiiann functions in this paper because of the permutation invariance of sets. Empirical evidence in our experiments supports this assertion.
We compare 4 different choices of i) Concatenation: directly concatenate different encodings into a fixed-length vector, with aero-nadding in empty positions. This not a cermutation-invaritt function and is therefore included as a limiting case for comparison. ii) Summation: the simplest cermutation-invariann method, which directly sums all encodings. iii) Weighted Summation: an extension of Summation that provides additional expressive power. Specifically, encodings are first mapped through a network to obtain element- wise weights (via softmax), which are then multiplied by the encodings themselves yield the result. iv) Attention: offering stronger expressive
41

---

Summation
Sunmatio
ESP
GMM
MAF
NICE
SUSPE
GMM
MAP
NICE
SOSPE
SIMPSON-NLIN
LARGEBD-NLIN
TRIANGLE
M
NAPKIN
Figure 14: Ablation study for Markov boundaries on 4 different settings of SCMs: (a) SIMPSON-
NLIN, (b) LARGEBD-NLIN, (c) M, (d) NAPKIN. Higher ESP indicates higher sampling efficiency.
ZZ
Concaeattion
Weishaed Sarmaation
Sunnation
Amtencia
GMM
MAP
NICE
SOSPE
GMM
MAP
NICE
SOSPE
H-a
14-2
NIMPSON-NIIN
b. LARGERD-NI IN
B
TRIANGLE-NIN
NAPKIN
Figure 15: Ablation study for Markov boundaries 4 different settings of SCMs: (a) SIMPSON-
NLIN, (b) LARGEBD-NLIN, (c) M, (d) NAPKIN. Lower FR indicates broader coverage.
power, Attention differs from Weighted Summation in that it uses encodings before input h rather than after output h, and new network a maps them to attention weights. Attention encoding is calculated through queries and keys softmax(h(2)a (g(x))), where we set the values in the attention mechanism to be constantly 1.
We no longer explore aggregation methods with stronger expressive power, as we found in our experiments that they are likely converge prematurely local optima. This causes the variance in the state space to increase as the optimization process progresses, which in turn prevents Bq. 8 from being uniformly applied across the entire state space.
Fig. Fig. 15 respecccvy illustrate the impact of selection on ESP and FR under different SCM and density estimation models. Overall, permutation-invariant methods (Summation, Weighted Sum- mation, Attention) consistently outperform permutation- variant methods (Concatenation), aligning with theoretical expectations and intuition. Among all methods, Attention consistently outperforms others, thus serving as the choice for all preceding experiments.
42

---

Table 4: Counterfactual density kemel (ō(y.)) estimation using different sampling methods on Ca extended for Tab. 2] Here, "O" represents the original SCM, and "P" represents the proxy
SCM. Regularization are applied to the dimensions. The main metric used FR ("-" indicates FR equals 1), with subscripts representing the regularized error bound.
SIMPSON-NLIN
TRIANGLE-NLIN LARGEBD-LIN
Method
SCM
[s]
3
[s| =
1
[s| =3
[s] =
O
0.89
0.84g 015
1.00p.000
0.99p.019
RS
P
0.900.012
0.85g.013
O
0.000.005
0.020.015
0.01o.021
0.11g.002
0.12a.013
0.00c.007
EXOM[MAF]
P
0.01 0.005
0. 0.400.012
0.61
0.1 11..008
0.50g.014
0.04g.007
O
0.00g.000
0.02p.016
0.010.046
0.08g.008
0. 0.
0.010.007
EXOM[NICE]
P
0.01o.005
0.40
0.610.018
0. 55.007
0.55g.018
0.08p.007
C.9
Counterfactual Estimation on Proxy SCMs
We combine identifiable proxy SCM with Exogenous Matching for counterfactual estimation. Specif- ically, we employ two distinct proxy SCM methods tailored for these scenarios to conduct counter- factual estimation two types of SCMs.
Counterfactual estimation with CausalNF We conduct experiments on 3 different Markovian diffeomorphic SCMs. According to the previous arguments, CausalNF belongs to IMM Pv.0), indi- cating that they are La-identifiable. In our experiments, we mainly measure the error counterfactual probability queries (( (y.) )), comparing the SCM represented by NF to the original SCM in 1024 different queries. These results are reflected in Tab. 4] with averages presented. For the original SCM, train EXOM five times and calculate the error bound for the same query. For CausalNF, we first train proxy SCMs using observational data, then train EXOM times for each model, calculate the error bound for the same queries, and finally take the average error bound across proxy SCMs.
Since in the continuous case it is challenging for RS to sample instances that satisfy the counterfactual event ,(y.), obtaining ground truth in experiments becomes difficult. Therefore, we primarily use FR as the main evaluation metric, while the error bound is employed to reflect the bias. When measuring error bounds, it is imperative to regularize the dimensionality of the counterfactual variables. Failure to do so may result in significant discrepancies in the final outcomes (estimates in high- digh-dimensional scenarios tend to be considerably lower than those in -iw-nimnnioonal ones, leading to a severe underestimation of errors in high-dimensional cases). Therefore, we calculate the error bounds as follows:
<( ( = 22 (ta((()) 3. (2.)
(122)
where 1. (>))) represents the estimate of P().) in the i-th experiment, regularized by the exponent 1 | which denotes the expectation in single dimension assuming independence across dimensions. The term std refers to the standard deviation. L = 5 the number of trials.
The results in Tab. 4 indicate that EXOM exhibits a relatively small error, comparable that of RS. Furthermore, when EXOM is applied to the proxy SCM, it yields an even smaller error (in part because all failure cases are excluded from consideration).
Counterfactual estimation with NCM We conducted experiments on Regional canonical SCMs. According to their proof, NCM exhibits duality with respect to the identifiability of the counterfactual query @ and the graphical causal identifiability. In our experiments, we focus primarily on the unbiasedness of EXOM combined with NCM. The results are shown in Tabs and
As RS is straightforward to sample in the discrete case, the ground truth can be computed. We used RS with 10 samples as ground truth (in the experiment, we compared against RS with 10 samples) and then computed the average bias. Additionally, since these queries only involve single dimension (regarding solely to Y), there is no need for dimensional regularization of the error bound.
In estimating these causal queries, EXOM yielded relatively small errors and showed even smaller errors when applied to proxy SCMs. This indicates that integration of EXOM with proxy SCMs
43

---

Table 5: Counterfactual effect estimation using different sampling methods on NCM, extended for
Tab. 2] Here, "O" represents the original SCM, and "P" represents the proxy SCM. The main metric
used the average bias w.r.t. the ground truth, with the subscript denotes the 95% CI error bound
over 5 trials
FAIRNESS
FAIRNESS-XW
Method
SCM
ATE
ETT
NDE
CUDE
ATB
ETT
NDB
CUDE
O
0.01g.013
0.01o.01
0.01g.m1
0.01g.020
0.01o.013
0.01p.0zr
0.01g.013
0.01g.223
RS
P
0.01g.013
0.01p.021
0.01a.014
0.01g.02s
0.01o.012
0.02a.osi
0.01p.013
0.02o.081
O
0.04g.005
0.06c.168
0.04a.181
0.060.161
0.05o.140
0.07a.213
0.050.149
0.07a.244
EXOMIMAF]
0.01 P
0.01
0.01
0.0
0.00p.011
0.02g.02s
0.00p.011
0.01
O
0.01g.oms
0.0
0.01g.020
0.01g.cov
0.010.000
0.02a.xs
0.02
0.04g.000
EXOMINICE]
P
0.01g.017
0.01
0.01
0.01g.ozz
0.010.012
0.02p.029
0.01p.014
0.01g.001
Table 6: (Tb.. continued) Counterfactual effect estimation using different sampling methods on NCM, extended for Tab. 2] Here, "O" represents the original SCM, and "P" represents the proxy
SCM. The main metric used is the average bias w.r.t. the ground truth, with the subscript denotes the
95% CI error bound over 5 trials
FAIRNESS-XY
FAIRNESS-Y YW
Method
SCM
ATE
ETT
NDE
CHDE
ATR
ETT
NDE
CHDE
O
0.01a.01s
0.01c.022
0.01a.015
0.01o.02s
0.006.012
0.01a.02s
0.01o.015
0.01a.027
RS
P
0.01g.012
0.01g.021
0.010.014
0.01p.019
0.01o.022
0.02a.oss
0.01p.017
0.01a.029
O
0.05p.1s0
0.09
0.04g.100
0.13p.335
0.04o 106
p.10g.
0.05
0.13 0.130.306
EXOM[MAF]
P
0.00p.009
0.01g.020
0.00o.013
0.0
0.01o.026
0.02s
0.01
0.01
O
0.01g.019
0.03p.111
0.05g.006
0.05g
0.030 129
0.06g.121
0.06o.150
0.1 0.1l
EXOMINICE]
P
0.01g.033
0.01p
0.00g.011
0.01p.054
0.0
0.01
0.01p
0.01g
is practical to address real-world problems. In particular, the errors in ETT and CtfDE within the original SCM were substantially higher than those in RS, primarily because their expressions include the denominator P(X), ), and inaccuracies in P(X) markedly influence their results. This matter will be summarized in App. D.3
D Discussions
D.1 Broader Impacts
Our work aims to enhance the efficiency of counterfactual reasoning and make it applicable to theoretically broader contexts. On the positive side, current counterfactual reasoning has a positive impact on society in terms of counterfactual fairness and counterfactual decision making. However, on the flip side, assumptions and premises are not applicable to real-world situations, counterfactual reasoning can still lead to potential biases, inequalities, and misuse. Therefore, careful consideration of the correctness of the assumptions and the appropriateness of usage scenarios required before applying counterfactual reasoning.
D.2 Counterfactual Density
For continuous distributions, is the density function rather than the probability measure that garners widespread attention and study within the relevant field. According to the definition of probability density, we can find its connection with the Lebesgue differentiation theorem, which directly provides method for density estimation:
Theorem 4 (Counterfactual Density Estimation). Assuming the counterfactual variables Y. are absolutely continuous on (Ωx,,, Σy., Py.), with probability density existing at y. E Ωy Let 5(y.) be function that maps y. cube (or ball) centered at y., and 5(y.) + y. denotes the cube's side length (or the ball's diameter) tending to O, then the following equation holds:
P(Y. E 8(y-))
P(Y, E 5(y.))
p(y=)
lim 5(y.)ty.
10(y+))
~
10(y))
(123)
where 18(y.)) ) is the volume (i.c., the Lebesgue measure) of 5(y.).
44

---

Proof. Directly inferred from the definition of the probability density function (Radon-Nikodym derivative) and the Lebesgue differentiation theorem.
This method is analogous to kernel density estimation (with the Parzen window as cube). However, the precision of the estimates provided by this method is highly dependent on (ō(y.)1, further highlighting the difficulty of satisfying the indicator function in the counterfactual probability measure. We indirectly investigated the performance of EXOM on this issue in our related experiments on continuous counterfactual estimation, where we chose Y. E (y²) as the counterfactual event for the estimation. The results of the related experiments confirmed the potential of EXOM in addressing this problem.
Of course, kernel-based density estimation has many limitations. We hope that Thm. 4 will serve as a lemma for future work, leading to the derivation of more suitable density estimation methods for continuous counterfactual distributions.
D.3 Limitations
1. Dependency counterfactual generation. In practical tasks, we typically only have access to observational or experimental data. In some settings, we also assume the acquisition of causal graphs or causal orderings through causal discovery and prior knowledge. However, the proposed method not directly applicable to such real-world settings, as indicated by the assumption [nc.. 3] which requires the ability for counterfactual generation. In practical scenarios, one approach is to pre-train a neural proxy SCM to perform the counterfactual generation task. Boundedness condition on importance weights. In the context learning and optimization, we introduce additional designs: the bounded importance weight condition. We assumed conditional distribution model to satisfy the boundedness of importance weight as much as possible. However, in practical applications, this depends on the specific forms of the exogenous distribution and the conditional distribution model, necessitating case-by-case analysis.
3. Faithfulness. For injection of the Markov boundary, we also introduce the faithfulness assumption, which assists us proving theorems and providing fast algorithms. However, this assumption might be violated, leading to more variables being included in the actual Markov blanket than our calculations predict. Additionally, modifications to autoregressive and coupling models do not fully align with the ideal effect. In fact, due to potential interactions within the hidden encodings of autoregressive and coupling models, even if each layer adheres to the Markov boundary dependencies, the aggregation of multiple layers may violate these dependencies. Recursiveness. Another potential assumption in this work is the recursive SCM, which might be violated in more general settings, such as systems with eycles. All the theoretical foundations of this work are based on recursive SCMs, thus requiring these theories to be extended to more general cases.
5. Variance of L3 expressions. In estimating counterfactual expressions Q E La. although the final result is unbiased, computing each term P()).) individually introduces significant variance, as demonstrated in the experiments on counterfactual query estimation with real SCMs (see Tabs 5 and 6 Therefore, future work should focus on reducing the estimation variance of specific queries.
Potential sensitivity to the specific forms of SCM in practice. With respect the parameter form of the SCM, although we impose no theoretical constraints on the form of SCM, experiments indicate that in some settings the form of SCM still significantly influences outcomes. For example, combinatorial experiments and ablation studies reveal unexpected scenarios, particularly with the combinations of long-chain SCMs (CHAIN-LIN-4, CHAIN-L.N--5, LARGEBD-NIIIN) As the causal chain lengthens, the performance of these density estimation models gradually weakens, suggesting that long causal chains complicate counterfactual inference. This phenomenon indicates that the form of SCM in practice affects the specific effectiveness of our methods, highlighting the necessity selecting appropriate density estimation models and making corresponding adjustments.
45

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All claims made in the abstract and introduction are reflected in the methodol- ogy (Sec. 3) and experiments (Sec. 5).
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. It is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
Justification: We briefly outline the limitations in Sec. 6land provide detailed discussion on the limitations in App. D.3]
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-texx system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms
and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faireess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We present all the assumptions in Sec. and complete the proofs in App. A Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
46

---

All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented
by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide a detailed description of the specific model architectures used in our experiments in App. C.4
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (c.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results. 5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We have included the complete code in the attachment, encompassing both training and validation processes.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/ public/guidess http:/uuudeessooeebmbiiiiie SoodeSuumisiionoool....
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not
47

---

reproducible, they including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
httpsccccpbplccccc ceciii/guiiiie http://///uieeeeeeeuububbisii CodeSubmissss wwwiibbmssssiiiii ssinniiiccy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why. At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended the paper) is recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide detailed experimental details in App. C
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: the experimental figures and tables, we present error bars, with 95% CI and quartile outliers.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula,
call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma eтror bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: In App. C.4 we describe the computation time and computational resources.
48

---

Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper). 9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and have confirmed that we not in violation of it.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., there a special consid-
eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We discussed the Broader Impacts in App. D.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinfornation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA]
Justification: The data and model presented in this paper do not pose risks misuse. Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring
49

---

that users adhere to usage guidelines or restrictions to access the model or implementing
safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors
should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We introduce the sources of code and data in App. C Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines: This paper does not release new assets.
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/mmodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects. Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
50

---

15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
51


Reviews:
[{"review_id": "LWxBWZsbKk", "reviewer": "Reviewer_8MDB", "summary": "Based on the importance sampling methods, the authors propose an exogenous matching approach to estimate counterfactual probability in general settings. They derive the variance upper bound of counterfactual estimators and transform it into the conditional learning problem. They also employ the Markov boundaries information in the inference to improve the learning performances further. Extensive experiments validate the superiority and practicality of their method.", "questions": "Please see my questions in Weaknesses.", "limitations": "Not applicable.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- This paper is clearly and well written.\n\n- The authors give a theoretical analysis of their estimator, its log-variance upper bound in the general settings, and the counterfactual Markov boundary, and they also perform extensive experiments to demonstrate effectiveness in several cases: with two types of stochastic counterfactual processes, with three categories of fully specified SCMs, etc.", "weaknesses": "- This paper makes it clear in lines 140-144 about the assumptions needed for the proposed method. Regarding assumption ii), I think it is not mild, and I am wondering if the proposed method would be sensitive to the specified distribution $P_{\\textbf{U}}$ of $\\textbf{U}$. The authors might have performed such experiments, but it is not quite clear. \n\n- Are the Markov boundaries learned from the observational data via the d-separation, or they are given prior? The authors claimed that such Markov boundaries are structural prior knowledge in line 223, whereas they gave Theorem 3 to demonstrate how to obtain them. \n\n- It is suggested to offer the whole procedure or pseudo code of their proposed algorithm somewhere.\n\n- In line 239, “augmentied” might be a typo error.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer’s positive evaluation and valuable feedback. Below are our responses to the questions raised in the Weaknesses section:\n\n> **Weakness 1**: Regarding assumption ii), if the proposed method would be sensitive to the specified distribution $P_\\mathbf{U}$ of $\\mathbf{U}$?\n\n**Response**: The proposed method is theoretically insensitive to the exogenous distribution $P_\\mathbf{U}$, as we do not impose any restrictions on the parameters or specific form of $P_\\mathbf{U}$. However, we acknowledge that in our experiments, we only considered specific exogenous distributions. This is because, in practice, neural proxy SCMs, as deep generative models, typically use a form-specific and relatively simple distribution for latent variables to facilitate inference.\n\n> **Weakness 2**: Are the Markov boundaries learned from the observational data via the d-separation, or they are given prior?\n\n**Response**: We would like to emphasize that the counterfactual Markov boundary (Def. 2) for counterfactual estimation tasks serves as prior knowledge, which is not necessary and is known before estimation. The Markov boundary here refers to the counterfactual Markov boundary of exogenous variable, which can be directly derived from other prior knowledge (i.e., the augmented graph) as described in Thm. 2 and Thm. 3. Specifically, if the augmented graph is provided, we can use graph algorithms (Thm. 3) to effectively derive the Markov boundary of any exogenous variable. For learning the counterfactual Markov boundary from data, some methods cited in the manuscript (e.g., 3, 94, 111) cannot be directly applied unless the data includes observations of the exogenous variables. We thank the reviewer for the question and will refine the text to clarify this information.\n\n> **Weakness 3, 4**: It is suggested to offer the whole procedure or pseudo code of their proposed algorithm somewhere; In line 239, “augmentied” might be a typo error.\n\nWe will carefully correct the typo errors in the manuscript and, in response to the reviewer's feedback, we will include pseudo code describing the proposed algorithm in the updated version.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "X4z3BYCgKl", "reviewer": "Reviewer_aexG", "summary": "This paper presents Exogenous Matching (EXOM), a new importance sampling method for estimating counterfactual probabilities in Structural Causal Models (SCMs). EXOM transforms variance minimization into a conditional distribution learning problem, providing an upper bound on counterfactual estimator variance as per Theorem 1. It outperforms existing methods across various SCM settings and integrates well with identifiable neural proxy SCMs for practical applications. By incorporating prior knowledge through Markov boundaries, EXOM further enhances performance, demonstrating its potential as an efficient tool for counterfactual estimation in diverse scenarios.", "questions": "1. In Table 1, the EXOM method with MAP shows significantly better performance on the SIMPSON-NLIN and NAPKIN datasets compared to EXOM with GMM, whereas the performance on the FAIRNESS-XW dataset is similar for both methods. Could the authors provide further explanation for this discrepancy? Why does the GMM approach underperform in these specific cases, and what factors contribute to the similar performance on FAIRNESS-XW?\n  \n2. In the ablation study investigating the impact of injecting Markov boundaries, could the authors please include the performance results of EXOM without Markov boundaries? This comparison would further illustrate the benefit of incorporating Markov boundaries.", "limitations": "The authors discuss the limitations of their work in Section 6 and Section D.3. Notably, this work does not have a negative societal impact.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. EXOM provides a tractable and efficient approach for counterfactual estimation in general settings, including scenarios with discrete or continuous exogenous variables and various observations and interventions. This flexibility makes it applicable to a wide range of causal inference problems.\n  \n2. The method is built on solid theoretical grounds, with the authors deriving an optimizable variance upper bound for counterfactual estimators.\n  \n3. The authors incorporate structural prior knowledge, specifically Markov boundaries, into the neural networks used for parameter optimization. They empirically validate the effectiveness of this approach across various scenarios.\n  \n4. EXOM consistently outperforms other importance sampling methods in various SCM settings, as demonstrated by the experimental results. Its compatibility with identifiable neural proxy SCMs further enhances its practical applicability.", "weaknesses": "1. Theorem 1 relies on the assumption that the density ratio $q(\\mathbf{u}|\\mathbf{y}_ {\\ast})/q(\\mathbf{u}|\\mathbf{y}_ {\\ast}^\\prime)\\leq \\kappa$ holds for all $\\mathbf{u} \\in \\Omega_{\\mathbf{u}}$ and $y_{\\ast}$, $y_ {\\ast} ^\\prime \\in \\Omega_{\\mathbf{Y}_{\\ast}}$. This assumption may be overly stringent, as probability measures with infinite support sets might easily violate it. Could the authors elaborate on this assumption and provide examples of distributions that satisfy it?\n  \n2. In the Sampling and Optimization section, the distribution of the exogenous variable $\\mathbf{U}$ is assumed to be known. However, in practical scenarios, $\\mathbf{U}$ is often unknown, necessitating additional efforts to estimate $\\mathbf{P_U}$ [A]​. Could the authors provide further clarification on this assumption and discuss potential methods for estimating $\\mathbf{P_U}$​?\n  \n3. I understand the authors only consider models that provide identifiability results. However, it is encouraged to include neural proxy SCM methods based on VAE and DDPM as experimental baselines. While these may lack identifiability guarantees, comparing against them would further illustrate the superiority of the proposed method in relation to current state-of-the-art techniques.\n  \n4. While the method shows good performance on the tested SCMs, it's unclear how well it scales to larger, more complex causal models. The experiments are conducted on relatively small SCMs, and scalability to high-dimensional or densely connected causal graphs isn't thoroughly addressed.\n  \n\n[A] Ren, Shaogang, and Xiaoning Qian. \"Causal Bayesian Optimization via Exogenous Distribution Learning.\" *arXiv preprint arXiv:2402.02277* (2024).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer's patience in reading our manuscript and the attention to issues concerning assumptions, generalization, scalability, and experimental aspects. Below are our responses to these concerns:\n\n> **Weakness 1**: Could the authors elaborate on the assumption about density ratio in Theorem 1 and provide examples of distributions that satisfy it?\n\n**Response**: We agree with this concern; this is indeed not a weak assumption, especially for probability measures with infinite support sets. This is because the proof of Thm. 1 relies on the linearity of Lebesgue integration, which requires Lebesgue integrable functions to be bounded. Some alternative formulations with more relaxed assumptions are as follows:\n\n1. **By Sufficient Condition**: If the importance weights satisfy $\\kappa^{-\\frac{1}{2}} \\le p(\\mathbf{u}) / q(\\mathbf{u} \\mid \\mathbf{y}_\\*) \\le \\kappa^\\frac{1}{2}$ with probability 1, then, almost surely, $q(\\mathbf{u} \\mid \\mathbf{y}'\\_\\*) / q(\\mathbf{u} \\mid \\mathbf{y}\\_\\*) \\le \\kappa$, and hence Thm. 1 holds.\n\n2. **By Concentration Inequality**: If the importance weights satisfy $\\kappa^{-\\frac{1}{2}} \\le p(\\mathbf{u}) / q(\\mathbf{u} \\mid \\mathbf{y}\\_\\*) \\le \\kappa^\\frac{1}{2}$ with probability $\\zeta < 1$, then by introducing the assumption that the second moment of the weights $\\mathbb{E}\\_w^2$ is bounded, and using Thm. 1 from [1], we can derive an inequality of the form $\\mathbb{V} \\le \\widehat{\\mathbb{V}} + f(n, \\zeta, \\kappa, \\delta, \\mathbb{E}\\_w^2)$ with probability $(1-\\delta) \\cdot \\zeta^n$, where $\\widehat{\\mathbb{V}}$ is computed from $n$ i.i.d. samples $\\mathbf{u}^{(1)}, \\dots, \\mathbf{u}^{(n)}$ with bounded weights. In this case, for $\\widehat{\\mathbb{V}}$, almost surely $q(\\mathbf{u}^{(i)} \\mid \\mathbf{y}'\\_\\*) / q(\\mathbf{u}^{(i)} \\mid \\mathbf{y}\\_\\*) \\le \\kappa$, so Thm. 1 holds for $\\widehat{\\mathbb{V}}$, and consequently, it also holds for $\\mathbb{V}$ by the inequality relationship.\n\n3. **By Approximation**: Similar to [2], consider only the approximation in bounded regions. If the importance weights satisfy $\\kappa^{-\\frac{1}{2}} \\le p(\\mathbf{u}) / q(\\mathbf{u} \\mid \\mathbf{y}\\_\\*) \\le \\kappa^\\frac{1}{2}$ in $\\Omega\\_{\\mathbf{U}}^+$, let $\\mathcal{Y}^\\*\\_+ = \\\\{\\mathbf{Y}^*(\\mathbf{u}) \\mid \\mathbf{u} \\in \\Omega\\_{\\mathbf{U}}^+ \\\\} \\cap \\mathcal{Y}\\_\\*$, and if $P(\\mathcal{Y}\\_\\*) \\approx P(\\mathcal{Y}^\\*\\_+)$, then the theorem holds for $\\mathcal{Y}^\\*\\_+$ and approximately for $\\mathcal{Y}\\_\\*$.\n\n> **Weakness 2**: Could the authors provide further clarification on the assumption that $P_\\mathbf{U}$ is known and discuss potential methods for estimating $\\mathbf{U}$?\n\n**Response**: We agree with the viewpoint that obtaining the exogenous distribution $P^*_\\mathbf{U}$ of the true SCM $\\mathcal{M}^*$ is challenging. Literature [A] addresses this by introducing additional parameter assumptions for the causal mechanism $\\mathcal{F}$ to recover the exogenous distribution. However, for counterfactual tasks, recovering the true exogenous distribution may not be necessary, since identifiable counterfactual results can still be provided by neural proxy SCMs. As deep generative models, the exogenous (or, latent) distribution $P_\\mathbf{U}$ is typically either pre-specified or trainable, so the assumption that $P_\\mathbf{U}$ is known is weak in this context.\n\n> **Weakness 3**: It is encouraged to include neural proxy SCM methods based on VAE and DDPM to further illustrate the superiority.\n\n**Response**: We thank the reviewer for the suggestion. Although in our experiments we only selected models that provide identifiability results as baselines, theoretically, the proposed method is not limited to the type of SCMs. We will consider applying it to neural proxy SCMs with VAE or DDPM as the backbone in future work for larger, more complex tasks.\n\n> **Weakness 4**: While the method shows good performance on the tested SCMs, it's unclear how well it scales to larger, more complex causal models.\n\n**Response**: In principle, we do not impose any restrictions on the size, complexity, or dimensionality of variables for SCM. However, in practice, it is well-known that computational complexity is a challenge when using normalizing flows for high-dimensional data. Therefore, scalability is an aspect worth investigating in future research.\n\n> **Question 1**: Could the authors provide further explanation for the discrepancy of performance on the SIMPSON-NLIN and NAPKIN datasets? Why does the GMM approach underperform in these specific cases, and what factors contribute to the similar performance on FAIRNESS-XW?\n\n**Response**: The FAIRNESS dataset is discrete and finite, whereas the SIMPSON-NLIN dataset is continuous and diffeomorphic. Thus, the sampling space of the former is smaller, which is why it did not show a significant difference in performance. The performance of GMM is as expected since its expressive capacity is weaker than that of MAF. We will include a breif discussion of these results in the main text to facilitate reader understanding.\n\n> **Question 2**: Could the authors please include the performance results of EXOM without Markov boundaries in the ablation study?\n\n**Response**: We have included results for \"EXOM without Markov boundaries\" (i.e., the blue bars) in both Fig. 4 and App. C.8 (Fig. 9), and compared them with \"EXOM with Markov boundaries\" (i.e., the orange bars). These results demonstrate the advantages of incorporating Markov boundaries. We will further improve the display of the legends in the figures to enhance readability.\n\n----\n[1] [Cortes, C., Mansour, Y., & Mohri, M. (2010). Learning Bounds for Importance Weighting.](https://papers.nips.cc/paper_files/paper/2010/hash/59c33016884a62116be975a9bb8257e3-Abstract.html)\\\n[2] [Tengchao Yu, Linjun Lu, & Jinglai Li. (2019). A weight-bounded importance sampling method for variance reduction.](https://arxiv.org/abs/1811.09436)", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7uawmz8Auh", "reviewer": "Reviewer_6eu3", "summary": "This paper introduces an importance sampling method for efficient estimation of counterfactual expressions within general settings. It transforms the variance minimization problem into a conditional distribution learning issue, allowing integration with existing modeling approaches. The paper also explores the impact of incorporating structural prior knowledge, i.e. Markov boundaries, and applies the method to identifiable proxy SCMs, proving the unbiasedness of estimates and illustrating the method's practical applicability.", "questions": "I would like to know whether Theorem 2 and 3 are novel. Have the papers cited (like 81, 3, 94, 111)  and other papers proposed methods to obtain Markov boundaries? If so, what is the improvement of the method proposed in this paper?", "limitations": "The authors adequately discussed the limitations in the paper.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The paper is well-structured, with many subsections and bullet points summarizing paragraphs.\n2. The approach proposed in this paper has clear intuition and is easy to implement.", "weaknesses": "1. Contributions are not disentangled well. All three points involve experimental or empirical findings.\n2. Some results of the ablation study are abnormal. First, the results show that the approach proposed is not robust. Under setting SIMPSON-NLIN and M, the inclusion of Markov Boundary Mask significantly improves ESP. However, under setting LARGEBD-NLIN and NAPKIN, the Markov Boundary Mask harms the performance, especially with backbone SOSPF. Second, the ESPs under setting LARGEBD-NLIN with backbone SOSPF are almost 0, even when $|s|=1$, which is hard to explain if including Markov Boundary Mask is effective. Third, the variance under setting LARGEBD-NLIN with backbone NICE is extremely large. \n3. Insufficient explanation or legend for figures, making it difficult for readers to understand. For example, in Figure 1, $\\mathcal{M}$ with a subscripted hammer is not explained. In Figure 3, the legend does not indicate what different colors mean.\n4. Hard to follow. Some terminologies need explanation or reference. For example, in line 234, *faithfulness* is not defined.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for highlighting the issues and providing valuable feedback to help us improve the manuscript. Below are our responses to the identified weaknesses and questions:\n\n> **Weakness 1**: Contributions are not disentangled well. All three points involve experimental or empirical findings.\n\n**Response**: We thank the reviewer for the suggestion. We will reorganize the presentation of our contributions by placing the experimental and empirical conclusions in the third point, while the first two points will separately address the two main contributions of our work: i) the introduction of EXOM and its theoretical guarantees, and ii) the injection of Markov boundaries into neural networks for conditioning.\n\n> **Weakness 2**: Some results of the ablation study are abnormal. First, the results show that the approach proposed is not robust. Second, it is hard to explain if including Markov Boundary Mask is effective under LARGEBD-NLIN with backbone SOSPF.Third, the variance under setting LARGEBD-NLIN with backbone NICE is extremely large.\n\n**Response**: We have found that the abnormality observed is not due to the lack of robustness in the method but rather is caused by the width of the hidden layers in the neural network. Specifically, the injection of Markov boundaries relies on masking connections in the neural network, which means that with the same width, the injection of Markov boundaries results in fewer weights being used during inference. In experiments where the performance was normal, the injection of Markov boundaries still exhibited superiority even with fewer weights. Moreover, as the Markov boundaries become smaller (as seen with NAPKIN and LARGEBD-NLIN), the masks become sparser, resulting in even fewer weights being involved in inference. When the number of hyperparameters to be inferred for the conditional model is very high (e.g., SOSPF), the neural network's representational capacity is significantly insufficient, leading to markedly poor performance.\n\nIn the previous ablation studies, to ensure fairness, we used a hidden layer width of 64 across all settings, but this was inadequate for SOSPF. We have now increased the width to 256, and the result has reversed, with the experimental results on SOSPF showing consistency with results from other settings, as detailed in the Author Rebuttal PDF (Fig. 1).\n\nIn the follow-up, we will also update the relevant results in the appendix. We will include results in App. C.7 and App. C.8 with a hidden layer width of 256, then add a discussion about the recommended hyperparameter settings and their theoretical rationale in App. C.4.\n\n> **Weakness 3**: Insufficient explanation or legend for figures, making it difficult for readers to understand. For example, in Figure 1, $\\mathcal{M}$ with a subscripted hammer is not explained. In Figure 3, the legend does not indicate what different colors mean.\n\n**Response**: In Fig. 1, $\\mathcal{M}$ with a subscripted hammer represents the SCM under intervention, i.e., the submodel. In Fig. 3, the different colors indicate different submodels. We will address fidelity after d-separation, where fidelity refers to the independence implied in the probability distribution corresponding to d-separation in the DAG.\n\n> **Weakness 4**: Hard to follow. Some terminologies need explanation or reference. For example, in line 234, faithfulness is not defined.\n\n**Response**: We apologize for the confusion caused by the omission of the conceptual explanation in the main text. We have explained the concept of faithfulness in the appendix (line 1115, 1116). We will revise the manuscript to include the discussion of faithfulness in the main text, following the section on d-separation. This is because faithfulness and d-separation are highly related; it signifies that every conditional independence present in the distribution is entailed by the d-separation in the DAG.\n\nRegarding the issue of lacking explanations for legends and terminology, we will thoroughly review the entire manuscript to identify any figures, concepts, or terms that may be missing explanations. We will ensure that these are properly addressed and clarified in the updated version.\n\n> **Question 1**: Whether Theorem 2 and 3 are novel? Have the papers cited (like 81, 3, 94, 111) and other papers proposed methods to obtain Markov boundaries? If so, what is the improvement of the method proposed in this paper?\n\n**Response**: In fact, the cited works and the proposed method are used for completely different tasks. The cited works are used to learn Markov boundaries among observable variables from data, while the proposed method is specifically designed to derive the counterfactual Markov boundaries (Def. 2) of exogenous variables from an augmented graph. The cited works cannot be directly applied to learn the counterfactual Markov boundary unless the exogenous variables are also observed in the data.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IfRfFLofWP", "reviewer": "Reviewer_rY6N", "summary": "The manuscript introduces exogenous matching, an importance sampling method for efficient estimation of counterfactual expressions in various settings. This method transforms the variance minimization problem into a conditional distribution learning problem, allowing integration with existing modeling approaches. The authors validate their theoretical findings through experiments with different Structural Causal Models (SCMs), showing competetive performance in a range of counterfactual estimation tasks. They also examine the impact of structural prior knowledge and demonstrate the method's unbiased estimates and practical applicability in identifiable proxy SCMs.\n\nUpdate: revising score upward following author rebuttal.", "questions": "When we say the causal model is “not fully specified”, does that just refer to the structural equations or to the graphical structure as well? In general, I was not always certain just how much causal information is used as input to this method.\n\nI’m a bit confused by Eq. 5. If this is meant to be a variance estimator, then presumably the RHS should be something like $\\mathbb{E}[X^2] – (\\mathbb{E}[X])^2$, where $X$ denotes the likelihood ratio $p(u) / q(u)$, correct? This is almost but not quite what we find here. Looking at the appendix, I don’t see why Eq. 48 follows from Eq. 47. Why do we get to drop the square from the first term?\n\nWhat does the constant $c$ denote in Eq. 6? Is it just the entropy of $P$? A brief word on this would help with intuition.\n\nShould there definitely be a negative both in front of the expectation *and* within it? I assume the first summand should just be $-\\mathbb{E} [ \\log Q(u \\mid Y_* (u)) ] $? This would look more like the classic cross entropy formula, as indicated in the following line. Eq. 7 also suggests this.\n\nOn “augmented graphs” – does this “reverse projection” of the ADMG always work? Different DAGs can have the same ADMG, for instance if two latent variables have all the same endogenous children. Perhaps there’s an unstated minimality assumption at work here?\n\nWhat is $m$ in Eq. 13?\n\nIt is not clear to me from Sect. 5 what the sample size and data dimensionality are for these tasks? In general, this section appears rushed. The performance metrics are also somewhat surprising. If data is simulated, then we presumably have ground truth with respect to counterfactual probabilities. If so, then why not just compute the mean square error of the proposed estimator, perhaps as a function of sample size?", "limitations": "Limitations are adequately addressed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The topic is timely and important, as counterfactual estimation has become an increasingly popular subject in statistics and machine learning. The proposal builds on recent results in neural causal models, specifically with normalizing flows. The ability to incorporate prior knowledge in the form of Markov boundaries is especially welcome, since computing counterfactuals is often intractable without such constraints. The theoretical results appear sound (though I confess I did not go closely through the proofs) and the empirical results are compelling.", "weaknesses": "The manuscript is not always clear, probably because a great deal of material has been moved to the appendix to accommodate page count. The result is a somewhat disjointed text that would likely be better served by a journal publication than a conference paper. That said, I am generally supportive of this submission and would be willing to revise my score upward if my questions are adequately addressed (see below).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the patience and valuable feedback. We acknowledge that, due to page limitations, we had to move some material to the appendix, which may have affected the clarity of the manuscript. We understand that this may have made the paper appear somewhat disorganized. To address this issue, we will reorganize the content to ensure a clearer and more coherent presentation. We will introduce more transitional text to improve the flow of the manuscript and make the key points more understandable.\n\nBelow are our responses to the questions:\n\n> **Question 1**: When we say the causal model is “not fully specified”, does that just refer to the structural equations or to the graphical structure as well? How much causal information is used as input to this method?\n\n**Response**: “Fully specified” means that $\\mathcal{F}$ is completely known, which means that we can write out the explicit expression of any $f \\in \\mathcal{F}$ (i.e., the structural equation). This also implies that the graphical structure is known, as the graph structure can be derived from $\\mathcal{F}$. On the other hand, “not fully specified but $\\mathcal{F}$ is evaluable” means that given an input $\\mathbf{u} \\in \\Omega_\\mathbf{U}$, we can obtain the output $\\mathcal{F}(\\mathbf{u})$. In other words, $\\mathcal{F}$ is allowed to be regarded as a black box, such as a neural network.\n\n> **Question 2**: Confusion by Eq. 5 and inadequate derivations in the appendix.\n\n**Response**: We apologize for the confusion caused by the omission of the relevant steps in the appendix. Specifically, the subscript of the expectation operator has changed:\n\n$\\mathbb{E}_Q\\left[\\left(\\frac{P(\\mathbf{u})}{Q(\\mathbf{u})}\\right)^2 f^2(\\mathbf{u})\\right]=\\sum\\_{\\mathbf{u}\\in\\Omega\\_{\\mathbf{U}}} Q(\\mathbf{u}) \\left(\\frac{P(\\mathbf{u})}{Q(\\mathbf{u})}\\right)^2 f^2(\\mathbf{u})$ $= \\sum\\_{\\mathbf{u} \\in \\Omega\\_{\\mathbf{U}}} \\frac{P^2(\\mathbf{u})}{Q(\\mathbf{u})} f^2(\\mathbf{u})$ $= \\sum\\_{\\mathbf{u} \\in \\Omega\\_{\\mathbf{U}}} P^2(\\mathbf{u}) \\left(\\frac{f^2(\\mathbf{u})}{Q(\\mathbf{u})}\\right)$ $= \\sum\\_{\\mathbf{u} \\in \\Omega\\_{\\mathbf{U}}} P(\\mathbf{u}) \\left(\\frac{P(\\mathbf{u})}{Q(\\mathbf{u})}\\right) f^2(\\mathbf{u})$ $= \\mathbb{E}_P\\left[\\frac{P(\\mathbf{u})}{Q(\\mathbf{u})} f^2(\\mathbf{u})\\right].$\n\nSimilar steps and the same conclusion also apply to continuous distributions.\n\n> **Question 3**: What does the constant $c$ denote in Eq. 6? Is it just the entropy of $P$? A brief word on this would help with intuition.\n\nAs described in the appendix (line 1018), $c$ is composed of the entropy of $P$ and $\\log\\kappa$. We will include an brief explanation of this in the main text.\n\n> **Question 4**: Typo error of $-\\mathbb{E}\\left[\\log q(\\mathbf{u}\\mid\\mathbf{Y}_*(\\mathbf{u}))\\right]$.\n\nWe thank the reviewer for the careful review and correction. This was indeed a typo, and we will fixed it in the updated version.\n\n> **Question 5**: On “augmented graphs” – does this “reverse projection” of the ADMG always work? Perhaps there’s an unstated minimality assumption at work here?\n\nIn recursive SCMs, there is not always a \"reverse projection\" between ADMGs and their corresponding augmented DAGs, as the reviewer mentioned, different DAGs may correspond to the same ADMG. We have tried to understand why the reviewer might have this question. We believe it may be due to some ambiguity in our presentation in the \"Markov boundary\" section, which might have led readers to think that \"ADMG -> DAG -> MB\" is implied. However, we only emphasize \"DAG -> MB\". We thank the reviewer for this question and will clarify this information.\n\n> **Question 6**: What is $m$ in Eq. 13?\n\nIn Appendix C.1, we explained $m$, but we forgot to explain it in the main text. We will add the missing explanation of \\$m$ in the main text, and the correct notation is $\\le m$ instead of $\\ge m$. Here, $m$ is a threshold value; if the proportion of valid samples in a test is less than $m$, it is considered a failure.\n\n> **Question 7**: It is not clear to me from Sect. 5 what the sample size and data dimensionality are for these tasks? In general, this section appears rushed. The performance metrics are also somewhat surprising. If data is simulated, then we presumably have ground truth with respect to counterfactual probabilities. If so, then why not just compute the mean square error of the proposed estimator, perhaps as a function of sample size?\n\nThe experimental settings in Section 5 are detailed in App. C.4. For sample sizes, RS was set to $10^6$, while the other methods used $10^3$. For variable dimensions, this is related to $|s|$ and the SCM settings; for example, LARGEBD-NLIN with 9 variables results in 45 dimensions when $|s|=5$. As for the choice of metrics, our use of ESP and LL is motivated by the intuition that indicator functions in counterfactual probability estimators are more challenging to satisfy (that is, equals $1$). Therefore, we refer to metrics from rare event sampling, which are computationally convenient, dimension-independent, and naturally range from 0 to 1. Another reason for using this metric is that as dimensions increase, counterfactual probabilities decrease exponentially, making it nearly impossible for RS to sample effective samples in high-dimensional cases. RS is a tractable method for estimating ground truth in general cases, so the ground truths are almost unattainable in such situations.\\\nConsidering the reviewer's feedback, we note that metrics based on the deviation from ground truth are suitable for counterfactual effect estimation, given the low dimensionality in the experimental setting. Therefore, we have updated the section \"Counterfactual Estimation on Proxy SCMs\" to use the bias with respect to the true values ($|\\hat{\\mathbb{E}} - \\mathbb{E}|$) as a metric of unbiasedness. The updated results can be found in the Author Rebuttal PDF (Tab. 1).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all the reviewers for their valuable feedbacks, which will help us improve the manuscript. Here, we summarize and address some common concerns, then list the changes to be made in the next updated manuscript.

### **Discussion on Assumptions**

- **Assumptions Required for Exogenous Matching**

    This refers to the necessary inputs for the Exogenous Matching algorithm. As indicated in lines 140-144 of the manuscript and as mentioned in 8MDB, these assumptions consist of two parts:

    1. An exogenous distribution $P_\mathbf{U}$ that is both sampleable and has computable density or probability.
    
        We do not impose any restrictions on its parameters or form. Please refer to our response to reviewer 8MDB for more details.
        
        Moreover, in practice, assuming the exogenous distribution is known is a weak assumption for using neural proxy SCMs for counterfactual estimation. See our response to reviewer aexG for further clarification.

    2. An evaluable causal mechanism $\mathcal{F}$.
        
        Here, "evaluable" means that given an input $\mathbf{u} \in \Omega_\mathbf{U}$, we can obtain the output $\mathcal{F}(\mathbf{u})$ without needing to know the specific expressions of all structural equations $f \in \mathcal{F}$. This leads to the concept of "fully specified," which is elaborated in our response to reviewer rY6N.

    These assumptions are weak for proxy SCMs; therefore, we have emphasized in the limitations section that our method does not directly estimate from data in practice but requires first training a proxy SCM, as demonstrated in our experiment "Counterfactual Estimation on Proxy SCMs."

- **Assumptions in Theorem 1**

    We acknowledge that the original description in the main text is not weak. To address this, we have introduced 3 more relaxed formulations. Please refer to our response to reviewer aexG for details. We will briefly incorporate these relaxations into the main text and provide their proofs in the appendix.


### **Clarification for Injecting Markov Boundaries**

In response to questions from reviewers 8MDB, 6eu3, and rY6N, we need to clarify some fundamental facts about "Injecting Markov Boundaries". Firstly, the "Markov Boundaries" here actually refer to counterfactual Markov Boundaries (Def. 2). Additionally, methods for learning Markov Boundaries from observational data cannot be directly applied to learning counterfactual Markov Boundaries. Therefore, we need Thm. 2 and Thm. 3 to derive counterfactual Markov Boundaries from the augmented graph. We will supplement this explanation in the updated version.

### **Discussion on Experiments**

- **Ablation**

    According to the feedback from reviewer 6eu3, the ablation experiments related to the Markov boundary exhibited anomalies, which could lead readers to perceive our method as not robust. We identified that these anomalies stem from the width of the hidden layers, which is a practical engineering issue and can be reasonably explained. For a detailed explanation, please refer to our response to reviewer 6eu3.

    We have updated the experiments accordingly, and the corresponding part of Fig. 4 in the main text will be revised to Fig. 1 in the Author Rebuttal PDF. We will also continue to refine the relevant experiments in the appendix, adding results with increased hidden layer widths, and include a discussion on hyperparameter selection and its impact.

- **Counterfactual Estimation on Proxy SCMs**

    Based on reviewer rY6N's suggestion, a more intuitive metric should be the deviation from the ground truth. We have elaborated on the intuition and rationale for choosing ESP and LL as metrics, which is due to the difficulty of obtaining the ground truth in high-dimensional settings in our experimental setup. For more details, please refer to our response to reviewer rY6N.

    Of course, for low-dimensional cases, such as in counterfactual effect estimation experiments, the ground truth is easily accessible. Therefore, we have updated the metrics for these experiments, and Tab. 2 in the main text will be revised to Tab. 1 in the Author Rebuttal PDF.

### **Discussion on Presentation**

Due to page limitations, we did indeed move some content to the appendix and inadvertently omitted some key information from the main text. We will carefully review and correct any typos, complete the missing explanations and results in the main text, and reorganize key sections with additional transitions to improve readability.

### **List of Changes**

Based on the reviewer’s suggestions, we will carry out the following changes:

- **Regarding the presentation in the main text:**
   1. Correct typo errors, reorganize the description of contributions, and ensure that legends clearly describe the elements in the figures. For example, eliminate the ambiguity caused by the missing legends in Fig. 1 and 2.
   2. Assess the inclusion of concepts, methods, and experimental results in the main text, and consider transferring relevant content from the appendix to the main text where appropriate, such as the explanation of the concept of faithfulness.
   3. Reorganize the presentation of some key sections, using smoother transitions or more intuitive ways to connect different parts to improve readability.

- **Regarding further clarification of assumptions:**
   1. Add additional discussion about the assumptions required for Exogenous Matching.
   2. Clarify the motivation and significance of (counterfactual) Markov boundaries, distinguish the proposed method from data-driven learning approaches.
   3. Add some relaxations to the assumptions in Thm. 1.

- **Regarding modifications to experimental figures:**
   1. Replace the corresponding part of Fig. 4 in the main text with Fig. 1 from the PDF.
   2. Replace Tab. 2 in the main text with Tab. 1 from the PDF.

Meta-review: This paper aims to estimate counterfactual probabilities in Structural Causal Models, and the authors proposed EXOM, an importance sampling approach to solve this task. All reviewers are inclined towards a consensus of accepting the paper. During the rebuttal phase, the authors have added supporting arguments and results to justify the proposed method.

This is a particularly long paper for a conference: the paper ends on page 54. I encourage the authors to submit a cleaner, more concise version of the camera-ready version. In fact, the reviewers all raise concerns about the paper's clarity and comment that the paper is not easy to follow. The authors should also include the discussion about the assumptions of Theorem 1 in the camera-ready version.
Decision: Accept (poster)