Partial Identification of Counterfactual Distributions
Year: 2021

Abstract:
This paper investigates the problem of bounding counterfactual queries from a 1 combination of observational data and qualitative assumptions about the underlying 2 data-generating model. These assumptions are usually represented in the form 3 of a causal diagram (Pearl, 1995). We show that all counterfactual distributions 4 (over finite observed variables) in an arbitrary causal diagram could be generated 5 by a special family of structural causal models (SCMs), compatible with the 6 same causal diagram, where unobserved (exogenous) variables are discrete, taking 7 values in a finite domain. This entails a reduction in which the space where the 8 original, arbitrary SCM lives can be mapped to a dual, more well-behaved space 9 where the exogenous variables are discrete, and more easily parametrizable. Using 10 this reduction, we translate the bounding problem in the original space into an 11 equivalent optimization program in the new space. Solving such programs leads to 12 optimal bounds over unknown counterfactuals. Finally, we develop effective Monte 13 Carlo algorithms to approximate these optimal bounds from a finite number of 14 observational data. Our algorithms are validated extensively on synthetic datasets. 15

text: This paper investigates the problem of bounding counterfactual queries from a1 combination of observational data and qualitative assumptions about the underlying2 data-generating model. These assumptions are usually represented in the form3 of a causal diagram (Pearl, 1995). We show that all counterfactual distributions4 (over finite observed variables) in an arbitrary causal diagram could be generated5 by a special family of structural causal models (SCMs), compatible with the6 same causal diagram, where unobserved (exogenous) variables are discrete, taking7 values in a finite domain. This entails a reduction in which the space where the8 original, arbitrary SCM lives can be mapped to a dual, more well-behaved space9 where the exogenous variables are discrete, and more easily parametrizable. Using10 this reduction, we translate the bounding problem in the original space into an11 equivalent optimization program in the new space. Solving such programs leads to12 optimal bounds over unknown counterfactuals. Finally, we develop effective Monte13 Carlo algorithms to approximate these optimal bounds from a finite number of14 observational data. Our algorithms are validated extensively on synthetic datasets.15
1 Introduction16
This paper studies the problem of inferring counterfactual queries from the combination of non-17 experimental data (e.g., observational studies) and qualitative assumptions about the data-generating18 process. These assumptions are represented in the form of a causal diagram [32], which is a19 directed acyclic graph where arrows indicate the potential existence of functional relationships among20 corresponding variables; some variables are unobserved. This problem arises in diverse fields such21 as artificial intelligence, statistics, cognitive science, economics, and the health and social sciences.22 For example, when investigating the gender discrimination in college admission, one may ask “what23 would the admission outcome be for a female applicant had she been a male?” Such a counterfactual24 query contains conflicting information: in the real world the applicant is female, in the hypothetical25 world she was not. Therefore, it is not immediately clear how to design effective experimental26 procedures for evaluating counterfactuals, let alone how to compute them from observations alone.27
The problem of identifying counterfactual distributions from the combination of data and a causal28 diagram has been studied in the causal inference literature. First, there exist a complete proof system29 for reasoning about counterfactual queries [19]. While such a system, in principle, is sufficient in30 evaluating any identifiable counterfactual expression, it lacks a proof guideline which determines the31 feasibility of such evaluation efficiently. There are algorithms to determine whether a counterfactual32 distribution is inferrable from all possible controlled experiments [41]. There exist also algorithms33 for identifying path-specific effects from experimental data [1] and observational data [42].34
In practice, however, the combination of quantitative knowledge and observed data does not always35 permit one to point-identify the target counterfactual queries. Partial identification methods concern36 with deriving informative bounds over the target counterfactual probability, even when the target37
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
itself is non-identifiable. Several algorithms have been developed to bound counterfactuals from the38 combination of observational and experimental data [30, 36, 3, 4, 14, 35, 23, 24, 16, 25, 49].39
In this work, we build on the approach introduced by Balke & Pearl in [3], which involves direct40 discretization of the exogenous domains, also referred to as the principal stratification [17, 34]. Con-41 sider the causal diagram of Fig. 1a, where X,Y, Z are binary variables in {0, 1}; U is an unobserved42 variable taking values in an arbitrary continuous domain. [3] showed that domains of U could be43 discretized into 16 equivalent classes without changing the original counterfactual distributions and44 the graphical structure in Fig. 1a. For instance, despite it being induced by an arbitrary distribution45 P ⇤(u) over a continuous domain of the exogenous variable U , the observational distribution P (x, y|z)46 must be reproduced by a generative model of the form P (x, y|z) = P u P (x|u, z)P (y|x, u)P (u),47 where P (u) is a discrete distribution over a finite exogenous domain {1, . . . , 16}.48
Using the finite-state representation of unobserved variables, [4] derived tight bounds on treatment49 effects under the condition of noncompliance in Fig. 1a. [11, 21] applied the parsimony of finite-state50 representation in a Bayesian framework, to obtain credible intervals for the posterior distribution of51 causal effects in noncompliance settings. Despite their optimal guarantees, these bounds are only52 applicable to the specific noncompliance setting in Fig. 1a. For the most general cases, a systematic53 procedure for bounding counterfactual queries in arbitrary causal diagrams is still missing.54
Our goal in this paper is to overcome these challenges. We investigate the expressive power of discrete55 structural causal models (SCMs) [33] where each unobserved variable is drawn from a discrete56 distribution, takes values in a finite set of states. We show that when inferring about counterfactual57 distributions (over finite observed variables) in an arbitrary causal diagram, one could restrict domains58 of unobserved variables to a finite space without loss of generality. This observation allows us to59 develop novel partial identification algorithms to bound unknown counterfactual probabilities from60 the observational data. More specifically, our contributions are as follows. (1) We introduce a61 special family of discrete SCMs, with finite unobserved domains, and show that it could represent62 all categorical counterfactual distributions in an arbitrary causal diagram. (2) Using this result, we63 translate the original partial identification task into equivalent polynomial programs. Solving such64 programs leads to informative bounds over unknown counterfactual probabilities, which are provably65 optimal. (3) We develop an effective Monte Carlo algorithm to approximate optimal counterfactual66 bounds from a finite number of observational data. Finally, our algorithms are validated extensively67 on synthetic datasets. Given space constraints, all proofs are provided in Appendices A and B.68
1.1 Preliminaries69
We introduce in this section some basic notations and definitions that will be used throughout the70 paper. We use capital letters to denote variables (X), small letters for their values (x) and ⌦X for71 their domains. For an arbitrary set X , let |X| be its cardinality. For convenience, we denote by P (x)72 probabilities P (X = x); for an arbitrary subdomain X ✓ ⌦X , P (X ) ⌘ P (X 2 X ). Finally, the73 indicator function 1X=x returns 1 if an event X = x holds true; otherwise 1X=x = 0.74
The basic semantical framework of our analysis rests on structural causal models (SCMs) [33,75 Ch. 7]. An SCM M is a tuple hV ,U ,F , P i where V is a set of endogenous variables and U is76 a set of exogenous variables. F is a set of functions where each fV 2 F decides values of an77 endogenous variable V 2 V taking as argument a combination of other variables in the system. That78 is, v fV (paV , uV ),PaV ✓ V , UV ✓ U . Exogenous variables U 2 U are mutually independent,79 values of which are drawn from the exogenous distribution P (u). Naturally, M induces a joint80 distribution P (v) over endogenous variables V , called the observational distribution. Each SCM81 is associated with a causal diagram G (e.g., Fig. 1), which is a directed acyclic graph (DAG) where82
solid nodes represent endogenous variables V , empty nodes represent exogenous variables U and83 arrows represent the arguments PaV , UV of each function fV .84
An intervention on an arbitrary subset X ✓ V , denoted by do(x), is an operation where values of85 X are set to constants x, regardless of how they are ordinarily determined. For an SCM M , let86 Mx denote a submodel of M induced by intervention do(x). For any subset Y ✓ V , the potential87 response Yx(u) is defined as the solution of Y in the submodel Mx given U = u. Drawing values88 of exogenous variables U following the probability measure P induces a counterfactual variable Yx.89 Specifically, the event Yx = y (for short, yx) can be read as “Y would be y had X been x”. For any90 subsets Y , . . . ,Z, X, . . . ,W ✓ V , the distribution over counterfactuals Yx, . . . ,Zw is defined as:91
P (yx, . . . , zw) =
Z
⌦U
1Yx(u)=y ^ · · · ^ 1Zw(u)=zdP (u). (1)
Distributions of the form P (yx) is called the interventional distribution; when the treatment set92 X = ;, P (y) coincides with the observational distribution. Throughout this paper, we assume93 that endogenous variables V are discrete and finite; while exogenous variables U could take any94 (continuous) value. The counterfactual distribution P (yx, . . . , zw) defined above is thus a categorical95 distribution. For a more detailed survey on SCMs, we refer readers to [33, Ch. 7].96
2 Discretization of Structural Causal Models97
For a DAG G with endogenous V and exogenous variables U , let P ⇤ denote the collection of all98 counterfactual distributions over variables V . Formally,99
P ⇤ = {P (yx, . . . , zw) | 8Y , . . . ,Z,X, . . . ,W ✓ V } . (2) Let M be the family of all the SCMs compatible with the causal diagram G, i.e., M =100 {8M | GM = G}1. Counterfactual distributions in G are defined as the collection {P ⇤M : 8M 2M }101 that contains all counterfactual probabilities induced by SCMs M in the candidate family M . In this102 section, we will show that counterfactual distributions in any causal diagram G could be generated by103 an alternative family of “generic” SCMs compatible with G, which we will define later.104 Definition 1 (Counterfactual-Equivalence). For a DAG G, let M ,N be two sets of SCMs compatible105 with G. M and N are said to be counterfactually equivalent (for short, ctf-equivalent) if for any106 M 2M , there exists an alternative N 2 N such that P ⇤M = P ⇤ N , and vice versa.107
Our analysis rests on a special family of SCMs where values of each exogenous variable are drawn108 from a discrete distribution over a finite set of states.109 Definition 2. An SCM M = hV ,U ,F , P i is said to be a discrete SCM if110
1. Values of every U 2 U are drawn from a discrete distribution P (u) over a domain ⌦U ; let111 ✓u denote the probability P (U = u), for any u 2 ⌦U . f112 2. Values of every V 2 V are decided by function v fV (paV , uV ) ⌘ ⇠ (paV ,uV ) V , where for113 8paV , uV , ⇠ (paV ,uV ) V is a constant in the finite domain ⌦V .114
Given a causal diagram G, our goal is to construct a family of discrete SCMs N that is counter-115 factually equivalent to the original family of SCMs M . Our construction utilizes a special type of116 clustering of nodes in the diagram, called the confounded component [45].117 Definition 3. For an DAG G, a subset C ✓ V is a c-component if any pair X,Y 2 C is connected118 in G by a bi-directed path of the form V1 $ V2 $ · · · $ Vn, n = 1, 2, . . . , where (1) V1 = X ,119 Vn = Y ; (2) {V1, . . . , Vn} ✓ V ; and (3) each Vi $ Vj is a sequence Vi Uk ! Vj and Uk 2 U .120
A c-component C in G is maximal if there exists no other c-component that contains C. We denote121 by C(G) the collection of all maximal c-components in G. Naturally, c-components in C(G) form a122 partition over endogenous variables V , which, in turn, defines a partition {[V 2CUV | 8C 2 C(G)}123 over exogenous variables U . Therefore, for every U 2 U , there must exist a unique c-component124 in C(G), denoted by CU , such that U 2 [V 2CUUV . For example, exogenous variables U1, U2 in125 Fig. 1a corresponds to c-components CU1 = {Z} and CU2 = {X,Y } respectively; while the causal126 diagram of Fig. 1b only has a single c-component {X,Y, Z}.127
1We will use the subscript M to represent the restriction to a specific SCM M . Therefore, GM represents the causal diagram associated with SCM M ; so does the collection of counterfactuals P ⇤M .
Theorem 1. For a DAG G, consider the following conditions2: (1) M is the set of all SCMs128 compatible with G; (2) N is the set of all discrete SCMs compatible with G where for every U 2 U ,129 its cardinality |⌦U | = Q V 2CU |⌦PaV 7! ⌦V |, i.e., the number of functions mapping from PaV to130 V for every variable V in the c-component CU . Then, M and N are counterfactually equivalent.131
Thm. 1 establishes the expressive power of discrete SCMs in representing counterfactual distributions132 in a causal diagram G. It implies that the counterfactual distribution P (yx, . . . , zw) in any SCM M133 could be generated using a generic model as follows, for dU = Q V 2CU |⌦PaV 7! ⌦V |,134
P (yx, . . . , zw) = X
U2U
X
u=1,...,dU
1Yx(u)=y ^ · · · ^ 1Zw(u)=z Y
U2U ✓u. (3)
Among above quantities, ✓u are parameters of the exogenous distribution P (u) over a finite domain135 {1, . . . , dU}. Counterfactual variables Yx(u) are recursively defined as follows:136
Yx(u) = {Yx(u) | 8Y 2 Y } , where Yx(u) =
( xY if Y 2X
⇠({Vx(u)|V 2PaY },uY )Y otherwise (4)
where xY is the value assigned to variable Y in constants x. As an example, consider the causal137 diagram G described in Fig. 1b where X,Y, Z are binary variables in {0, 1}. Since G has a single c-138 component {X,Y, Z}, exogenous variables U1, U2 must share the same cardinality d in the proposed139 family of discrete SCMs N . It follows from Thm. 1 the counterfactual distribution P (z, xz0 , yx0) in140 any SCM compatible with G could be written as follows:141
P (z, xz0 , yx0) = dX
u1,u2=1
1 ⇠ (u1) Z =z ^ 1 ⇠ (z0,u1,u2) X =x ^ 1 ⇠ (x0,u2) Y =y ✓u1✓u2 , (5)
where ⇠(u1)Z , ⇠ (z,u1,u2) X , ⇠ (x,u2) Y are parameters taking values in {0, 1}; ✓ui , i = 1, 2, are probabilities142 of the discrete distribution P (ui) over the finite domain {1, . . . , d}. The cardinality d = |⌦Z | ⇥143 |⌦Z 7! ⌦X |⇥ |⌦X 7! ⌦Y | = 32. The total cardinalities of domains for U1, U2 are thus 2d = 64.144
Comparison with related work One could naïvely apply the discretization procedure in [3] and145 obtain a family of discrete SCMs that are sufficient in representing distributions in an causal diagram.146 However, such parametrization is not necessarily complete. To witness, consider again the causal147 diagram in Fig. 1b with binary X,Y, Z. Applying the discretization in [3] leads to a family of discrete148 SCMs compatible with a different diagram in Fig. 1c where the cardinality of exogenous variable149 U is equal to d = 32 (see Appendix D for details). However, this parametrization fails to capture150 some critical constraints over counterfactual distributions since it does not maintain the original151 structure of the causal diagram. For instance, counterfactual variables Z and Yx in the original152 diagram of Fig. 1b are independent due to independence restrictions [33, Ch. 7.3.2]; while Z and153 Yx in Fig. 1c are generally correlated due to the presence of unobserved confounder U . Compared154 with [3], the discretization method in Thm. 1 captures all constraints over counterfactual distributions155 while requiring only a factor of |U | increase in the cardinality of exogenous domains.156
More recently, [15] proved a special case of Thm. 1 for interventional distributions in a specific157 class of causal diagrams that satisfy the running intersection property. When there is no direct arrow158 between endogenous variables, [38] showed that the observational distribution in a diagram could be159 represented using finite-state exogenous variables. Thm. 1 generalizes these results by showing that,160 for the first time, all counterfactual distributions in an arbitrary causal diagram could be generated161 using discrete exogenous variables taking values from a finite domain, without any loss of generality.162
2.1 Partial identification of Counterfactual Distributions163
To demonstrate the expressive power of discrete SCMs, we investigate the problem of partial iden-164 tification of counterfactual distributions. For an SCM M⇤ = hV ,U ,F , P i, we are interested in165 evaluating an arbitrary counterfactual probability P (yx, . . . , zw). The detailed parametrization of166 M⇤ is unknown. Instead, the learner only has access to the causal diagram G and the observa-167 tional distribution P (v) induced by M⇤. Our goal is to derive an informative bound [l, r] from the168 combination of G and P (v) that contains the actual counterfactual probability P (yx, . . . , zw).169
2For every V 2 V , ⌦PaV 7! ⌦V is the set of all functions mapping from domains ⌦PaV to ⌦V .
Let N denote the family of discrete SCMs defined in Thm. 1 which are compatible with the causal170 diagram G. We derive a bound [l, r] over P (yx, . . . , zw) from the observational data P (v) by solving171 the following optimization problem:172
[l, r] = min /max n PN (yx, . . . , zw) | 8N 2 N , PN (v) = P (v) o (6)
For instance, consider again the double-bow diagram G in Fig. 1b. The observational distribution173 P (x, y, z) in any discrete SCM in N could be written as:174
P (x, y, z) = dX
u1,u2=1
1 ⇠ (u1) Z =z ^ 1 ⇠ (z,u1,u2) X =x ^ 1 ⇠ (x,u2) Y =y ✓u1✓u2 . (7)
One could derive a bound over the counterfactual distribution P (z, xz0 , yx0) from the observational175 data P (x, y, z) by solving polynomial programs which optimize the objective Eq. (5) over parameters176 ✓u1 , ✓u2 , ⇠ (u1) Z , ⇠ (z,u1,u2) X , ⇠ (x,u2) Y , subject to the observational constraints Eq. (7).177
As a corollary, it follows immediately from Thm. 1 that the solution [l, r] of the optimization problem178 Eq. (6) is guaranteed to be a valid bound over the unknown counterfactual P (yx, . . . , zw).179 Corollary 1 (Soundness). Given a DAG G and an observational distribution P (v), let M be the set180 of all SCMs compatible with G and let Mo = {8M 2M | PM (v) = P (v)}. For the solution [l, r]181 of Eq. (6), PM (yx, . . . , zw) 2 [l, r] for any SCM M 2Mo.182
Since the underlying SCM M⇤ 2Mo, Corol. 1 implies that the derived bound [l, r] must contain the183 actual counterfactual probability P (yx, . . . , zw). Our next result shows that such a bound [l, r] is184 provably tight, i.e., it cannot be improved without additional assumptions.185 Corollary 2 (Tightness). Given a DAG G and an observational distribution P (v), let M be the set186 of all SCMs compatible with G and let Mo = {8M 2M | PM (v) = P (v)}. For the solution [l, r]187 of Eq. (6), there exist SCMs M1,M2 2Mo such that PM1(yx, . . . , zw) = l, PM2(yx, . . . , zw) = r.188
Corol. 2 confirms the tightness of the bound [l, r] obtained from Eq. (6). Suppose there exists a valid189 bound [l0, r0] strictly contained in [l, r]. One could construct from Corol. 2 an SCM M compatible190 with the causal diagram G and the observational distribution P (v), but its counterfactual probability191 P (yx, . . . , zw) lies outside [l0, r0], which is a contradiction.192
The optimization problem of Eq. (6) is reducible to equivalent polynomial programs (see Appendix E).193 Despite the soundness and tightness of derived bounds, solving such programs may take exponentially194 long in the most general case [29]. Our focus here is upon the causal inference aspect of the problem195 and like earlier discussions we do not specify which solvers are used [3, 4]. In some cases of196 interest, effective approximate planning methods for polynomial programs do exist. Investigating197 these methods is an ongoing subject of research [26, 31, 48, 28, 27].198
3 Bayesian Approach for Partial Identification199
This section describes an effective algorithm to approximate the optimal counterfactual bound in200 Eq. (6), provided with finite samples v̄ = v(n)
N n=1
drawn from the observational distribution201 P (v), and prior distributions over parameters ✓u and ⇠ (paV ,uV ) V (possibly uninformative).202
We first introduce Markov Chain Monte Carlo (MCMC) algorithms that sample the posterior distribu-203 tion P (✓ctf | v̄) over a counterfactual probability ✓ctf = P (yx, . . . , zw). More specifically, for every204 V 2 V , 8paV , uV , parameters ⇠ (paV ,uV ) V are drawn uniformly over the finite domain ⌦V . For every205 U 2 U , exogenous probabilities ✓u are drawn from a generalized Dirichlet distribution [12]. We will206 take the view of a stick-breaking construction [40] which successively breaks pieces off a unit-length207 stick with size proportional to random draws from a Beta distribution. Parameters ✓u are proportions208 of each of the pieces relative to its original size. Formally,209
8u = 1, 2, . . . , dU , ✓u = µu
u 1Y
i=1
(1 µi), µu ⇠ Beta ⇣ ↵(u)U , (u) U ⌘ , (8)
For every endogenous variable V 2 V , ⇠V =
n
⇠(paV ,uV )V | 8paV , uV
o
.
where dU = Q V 2CU |⌦PaV 7! ⌦V | and ↵ (u) U , (u) U > 0 are hyperparameters. Finally, we truncate210 this construction by setting µdU = 1. Note from Eq. (8) that all parameters ✓u for u > dU are equal211 to zero. As an example, Fig. 2 shows a graphical representation of the data-generating process over212 parameters ✓u and ⇠ (paV ,uV ) V associated with SCMs in Fig. 1b, spanning over N observations.213
Gibbs sampling is a well-known MCMC algorithm that allows one to sample posterior distributions.214 For convenience, we introduce the following notations. Let parameters ✓ = {✓u | 8U 2 U , 8u}215 and ⇠ = n ⇠(paV ,uV )V | 8V 2 V , 8paV , uV o . The set Ū = U (n) N n=1 are exogenous variables216
affecting N observations V̄ = V (n) N n=1
; we use ū to represent their realizations. Our blocked217 Gibbs sampler works by iteratively drawing values from the conditional distributions of variables as218 follows [22]. Detailed derivations of complete conditional distributions are shown in Appendix F.219
Sampling P (ū | v̄,✓, ⇠). Exogenous variables U (n), n = 1, . . . , N , are mutually independent220 given parameters ✓, ⇠. We could draw each U (n) | ✓, ⇠, V̄ corresponding to the nth observation221
independently. The complete conditional for U (n) is given by222
P ⇣ u(n) | v(n),✓, ⇠ ⌘ / Y
V 2V 1 ⇠ (pa(n)V ,u (n) V ) V =v (n)
Y
U2U ✓u. (9)
Sampling P (⇠,✓ | v̄, ū). Parameters ⇠,✓ are independent given V̄ , Ū . Therefore, we will derive223 complete conditional ⇠,✓ separately. Note that in discrete SCMs, the nth observation of variable224 V 2 V is decided by v(n) ⇠(paV ,uV )V given pa (n) V = paV , u (n) V = uV . Thus, draw values of each225 ⇠(paV ,uV )V 2 ⇠ from the complete conditional defined as:226
P ⇣ ⇠(paV ,uV )V | v̄, ū ⌘ = ( 1 ⇠ (paV ,uV ) V =v (i) if 9i, s.t. pa (i) V = paV , u (i) V = uV ,
1/|⌦V | otherwise. . (10)
Let nu = PN n=1 1u(n)=u records the number of values in u (n) that are equal to u. By the conjugacy227 of the generalized Dirichlet distribution, the complete conditional of ✓u is given by, for every U 2 U ,228
8u = 1, 2, . . . dU , ✓u = µu
u 1Y
i=1
(1 µi), µu ⇠ Beta ↵(u)U + nu, (u) U +
dUX
k=u+1
nk
! . (11)
Doing so eventually produces values drawn from the posterior distribution over ✓, ⇠, Ū | V̄ . Given229
parameters ✓, ⇠, we compute the counterfactual probability ✓ctf = P (yx, . . . , zw) following the230 three-step algorithm in [33] which consists of abduction, action, and prediction. Thus computing ✓ctf231 from each draw ✓, ⇠, Ū eventually gives us the draw from the posterior distribution P (✓ctf | v̄).232
3.1 Collapsed Gibbs Sampling233
We also describe an alternative sampler that applies to stick-breaking priors with a known Pólya234 urn characterization. Formally, consider stick-breaking priors in Eq. (8) with hyperparameters235
↵(u)U = ↵U/dU and (u) U = (dU u)↵U/dU for some real ↵U > 0. Let Ū n denote the set236 difference Ū \U (n); so does V̄ n = V̄ \ V (n). Our collapsed Gibbs sampler first iteratively draws237 values from the conditional distribution of U (n) | Ū n, V̄ , n = 1, . . . , N , as follows.238

heading: Sampling P
text: u(n) | v̄, ū n . At each iteration, draw U (n) from the conditional given by239
P ⇣ u(n) | v̄, ū n ⌘ / Y V 2V P ⇣ v(n) | pa(n)V , u (n) V , v̄ n, ū n ⌘ Y U2U P ⇣ u(n) | v̄ n, ū n ⌘ . (12)
Among quantities in the above equation, for every V 2 V ,240
P ⇣ v(n) | pa(n)V , u (n) V , v̄ n, ū n ⌘ =
( 1v(n)=v(i) if 9i 6= n, pa (i) V = pa (n) V , u (i) V = u (n) V ,
1/|⌦V | otherwise. . (13)
For every U 2 U , let ū n be a set of exogenous samples u(1), . . . , u(n 1), u(n+1), . . . , u(N) . Let241
{u⇤1, . . . , u ⇤ K} denote K unique values that samples in ū n take on.242
P ⇣ u(n) | v̄ n, ū n ⌘ =
8 >><
>>:
n⇤k + ↵U/dU ↵U +N 1
if u(n) = u⇤k, for k = 1, . . . ,K
↵U (1 K/dU )
↵U +N 1 if u(n) 62 {u⇤1, . . . , u ⇤ K}
. (14)
where n⇤k = P i 6=n 1u(i)=u⇤k records the number of values in u (i) 2 ū n that are equal to u⇤k.243
Doing so eventually produces exogenous variables drawn from the posterior distribution of Ū | V̄ .244
We then sample parameters from the posterior distribution of ✓, ⇠ | Ū , V̄ ; the complete conditional245
P (⇠,✓ | v̄, ū) are given in Eqs. (10) and (11). Finally, computing ✓ctf from each sample ✓, ⇠ gives246 us a draw from the posterior distribution P (✓ctf | v̄).247
When the cardinality dU of exogenous domains is high, the collapsed Gibbs sampler described here is248 more computational efficient than the blocked sampler, since it does not iteratively draw parameters249 ✓, ⇠ in the high-dimensional space. Instead, the collapsed sampler only draws ✓, ⇠ once after samples250 drawn from the distribution of Ū | V̄ converge. On the other hand, when the cardinality dU is251 reasonably low, the blocked Gibbs sampler is preferable since it exhibits better convergence [22].252
3.2 Credible Intervals over Counterfactual Probabilities253
Given a MCMC sampler, one could bound the counterfactual probability ✓ctf by computing credible254 intervals from the posterior distribution P (✓ctf | v̄).255 Definition 4. Fix ↵ 2 [0, 1). A 100(1 ↵)% credible interval [l↵, r↵] for ✓ctf is given by256
l↵ = sup {x | P (✓ctf  x | v̄) = ↵/2} , r↵ = inf {x | P (✓ctf  x | v̄) = 1 ↵/2} . (15)
For a 100(1 ↵)% credible interval [l↵, r↵], any counterfactual probability ✓ctf that is compatible257 with observational data v̄ lies between the interval l↵ and r↵ with probability 1 ↵. Credible258 intervals have been widely applied for computing bounds over counterfactuals provided with finite259 observations [20, 47, 37, 8, 46]. As the number of observational data N grows (to infinite), the 100%260 credible interval [l0, r0] eventually converges to the optimal asymptotic bound [l, r] in Eq. (6) [11].261
Let ✓(t) T t=1
be T samples drawn from P (✓ctf | v̄). One could compute the 100(1 ↵)% credible262 interval for ✓ctf using the following consistent estimators [39]:263
l̂↵(T ) = ✓ (d(↵/2)Te), r̂↵(T ) = ✓ (d(1 ↵/2)Te), (16)
where ✓(d(↵/2)Te), ✓(d(1 ↵/2)Te) are the d(↵/2)T eth smallest and the d(1 ↵/2)T eth smallest of264 ✓(t)
3. Our next results establish non-asymptotic deviation bounds for the empirical estimates of265 credible intervals defined in Eq. (16) for finite samples.266 Lemma 1. Fix T > 0 and 2 (0, 1). Let function f(T, ) = p 2T 1 ln(4/ ). With probability at267
least 1 , estimators l̂↵(T ), r̂↵(T ) for any ↵ 2 [0, 1) is bounded by268
l̂↵(T ) 2 ⇥ l↵ f(T, ), l↵+f(T, ) ⇤ , r̂↵(T ) 2 ⇥ r↵+f(T, ), r↵ f(T, ) ⇤ . (17)
3For any real ↵ 2 R, d↵e denotes the smallest integer n 2 Z larger than ↵, i.e., d↵e = min{n 2 Z | n ↵}.
Algorithm 1: CREDIBLEINTERVAL 1: Input: Credible level ↵, tolerance level , ✏. 2: Output: An credible interval [l↵, h↵] for ✓ctf. 3: Let T = d2✏ 2 ln(4/ )e. 4: Draw samples ✓(1), . . . , ✓(T ) from the
posterior distribution P (✓ctf | v̄). 5: Return interval h l̂↵(T ), r̂↵(T ) i (Eq. (16)). We summarize our algorithm, CREDIBLEIN-269 TERVAL, in Alg. 1. It takes a credible level270 ↵ and tolerance levels , ✏ as inputs. In par-271 ticular, CREDIBLEINTERVAL repeatedly draw272 T d2✏ 2 ln(4/ )e samples from P (✓ctf | v̄).273 It then computes estimates l̂↵(T ), ĥ↵(T ) from274 drawn samples following Eq. (16) and return275 them as the output. It follows immediately from276 Lem. 1 that such a procedure efficiently approx-277 imates a 100(1 ↵)% credible interval.278
Corollary 3. Fix 2 (0, 1) and ✏ > 0. With probability at least 1 , the interval [l̂, r̂] =279 CREDIBLEINTERVAL(↵, , ✏) for any ↵ 2 [0, 1) is bounded by l̂ 2 [l↵ ✏, l↵+✏] and r̂ 2 [r↵+✏, r↵ ✏].280
Corol. 3 implies that any counterfactual parameter ✓ctf compatible with observational data v̄ falls281 between [l̂, r̂] = CREDIBLEINTERVAL(↵, , ✏) with probability P ⇣ ✓ctf 2 [l̂, r̂] | v̄ ⌘ ⇡ 1 ↵± ✏. As282
the tolerance rate ✏! 0, [l̂, r̂] converges to a 100(1 ↵)% credible interval with high probability.283
4 Simulations and Experiments284
We demonstrate our algorithms on various simulated SCM instances and a real world patient dataset285 collected from the International Stroke Trial (IST) [10]. Overall, we found that simulation results sup-286 port our findings and the proposed bounding strategy consistently dominates state-of-art algorithms.287 When target distributions are identifiable (Experiment 1), our bounds collapse to the actual, unknown288 counterfactual probabilities. For non-identifiable settings, our algorithm obtains sharp asymptotic289 bounds when closed-form solutions already exist (Experiments 2 & 3); and improves over state-of-art290 bounds in other more general cases where the optimal strategy is unknown (Experiment 4).291
In all experiments, we evaluate our proposed bounding strategy based on credible intervals (ci). In292 particular, we draw 4 ⇥ 103 samples from the posterior distribution over the target counterfactual293 ✓ctf | V̄ . This allows us to compute 100% credible interval over ✓ctf within error ✏ = 0.05, with294 probability at least 1 = 0.95. As the baseline, we also include the actual counterfactual probability295 ✓⇤. For details on simulation setups and additional experiments, we refer readers to Appendix C.296
Experiment 1: Frontdoor Graph This experiment evaluates our sam-297 pling algorithm on interventional probabilities that are identifiable from298 the observational data. Consider the “Frontdoor” graph described in299 Fig. 3 where X,Y,W are binary variables in {0, 1}; U1, U2 2 R. In this300 case, the interventional distribution P (yx) is identifiable from P (x,w, y)301 through the frontdoor adjustment [33, Thm. 3.3.4]. We collect N = 105302 observational samples V̄ = {X(n), Y (n),W (n)}Nn=1 from a randomly303
generated SCM. Fig. 4a shows samples drawn from the posterior distribution of the target probability304 P (Yx=0 = 1) | V̄ . The analysis reveals that these samples collapse to the actual interventional305 probability P (Yx=0 = 1) = 0.5085, which confirms the identifiability of P (yx) in Fig. 3.306
Experiment 2: Instrumental Variables (IV) This experiment evaluates our bounding strategy in307 non-identifiable settings, while closed-form solutions for the optimal bounds over target probabilities308 already exist. Consider first the “IV” diagram in Fig. 1a where X,Y, Z 2 {0, 1} and U1, U2 2 R.309 The non-identifiability of P (yx) from the observational data P (x, y, z) with the instrument Z and the310 unobserved confounding between X and Y has been acknowledged in [5]. For binary X,Y, Z, [2]311 derived closed-form, sharp bounds over P (yx) (labelled as opt). We collect N = 105 observational312 samples V̄ = {X(n), Y (n), Z(n)}Nn=1 from a randomly generated SCM instance. Fig. 4b shows313 samples drawn from the posterior distribution of
P (Yx=0 = 1) | V̄ . As a baseline, we also include314 the optimal bound opt, and posterior samples obtained from the Gibbs sampler of [11], which utilizes315 the canonical partitions of exogenous domains in [2] (bp). The analysis reveals that our algorithm316 derives the valid bound over the actual probability P (Yx=0 = 1) = 0.3954; the 100% credible317 interval converges to the optimal IV bound l = 0.1468, r = 0.6617.318
Experiment 3: Probability of Necessity and Sufficiency (PNS) We now study the problem of319 evaluating the probability of necessity and sufficiency P (Yx=1 = 1, Yx=0 = 0) from the observational320 data P (x, y) in the “Bow” diagram of Fig. 1d where X,Y 2 {0, 1} and U 2 R. The sharp bound for321 P (Yx=1 = 1, Yx=0 = 0) from P (x, y) was introduced in [44] (labelled as opt). We collect N = 105322 observational samples V̄ = {X(n), Y (n)}Nn=1 from an SCM instance. Fig. 4c shows samples drawn323 from the posterior distribution of
P (Yx=1 = 1, Yx=0 = 0) | V̄ . As a baseline, we also include the324 optimal bound opt, and posterior samples obtained from the Gibbs sampler which discretizes the325 exogenous domains using canonical partitions [2] (bp). The analysis reveals that our 100% credible326 interval (ci) matches the optimal PNS bound l = 0, r = 0.6775, i.e., the proposed strategy achieves327 the sharp bound over the counterfactual probability P (Yx=1 = 1, Yx=0 = 0) = 0.1867.328
Experiment 4: International Stroke Trials (IST) IST was a large, randomized, open trial of up329 to 14 days of antithrombotic therapy after stroke onset [10]. In particular, the treatment X is a pair330 (i, j) where i = 0 stands for no aspirin allocation, 1 otherwise; j = 0 stands for no heparin allocation,331 1 for median-dosage, and 2 for high-dosage. The primary outcome Y 2 {0, . . . , 3} is the health332 of the patient 6 months after the treatment, where 0 stands for death, 1 for being dependent on the333 family, 2 for the partial recovery, and 3 for the full recovery.334
To emulate the presence of unobserved confounding, we filter the experimental data with selection335 rules f (Z)X , Z 2 {0, . . . , 9}, following a procedure in [49]. Doing so allows us to obtain N = 3⇥10 3336 synthetic observational samples V̄ = {X(n), Y (n), Z(n)}Nn=1 that are compatible with the “Double337 bow” diagram of Fig. 1b. We are interested in evaluating the treatment effect E[Yx=(1,0)] for338 only assigning aspirin X = (1, 0). Fig. 4d shows samples drawn from the posterior distribution339 of E[Yx=(1,0)] | V̄ . As a baseline, we also include a naïve generalization of the discretization340 procedure (bp) [2] (see Appendix D) and the natural bounds [36, 30] estimated at the 95% confidence341 level (nb) [49]. Posterior samples of ci and bp are drawn using our proposed collapsed sampler342 due to the high-dimensional latent space. The analysis reveals that all algorithms achieve bounds343 that contain the actual, target causal effect E[Yx=(1,0)] = 1.3418. Our bounding strategy obtains a344 100% credible interval lci = 1.2604, rci = 1.4687, which consistently improves over all the other345 algorithms (lbp = 1.1121, rbp = 1.8073, lnb = 1.1195, rnb = 1.6221).346
5 Conclusion347
This paper investigated the problem of partial identification of counterfactual distributions, which348 concerns with bounding unknown counterfactual probabilities from the combination of the obser-349 vational data and qualitative assumptions of the data-generating process, represented in the form of350 a directed acyclic causal diagram. We studied a special family of SCMs with discrete exogenous351 variables, taking values from a finite set of unobserved states, and showed that it could represent all352 counterfactual distributions (over finite observed variables) in an arbitrary causal diagram. That is,353 this new family of discrete SCMs is counterfactual equivalent to the original family of candidate354 SCMs compatible with the causal diagram. Using this result, we developed a novel algorithm to355 derive bounds over counterfactual probabilities from finite observations, which are provably tight.356
References357
[1] C. Avin, I. Shpitser, and J. Pearl. Identifiability of path-specific effects. In Proceedings of the358 Nineteenth International Joint Conference on Artificial Intelligence IJCAI-05, pages 357–363,359 Edinburgh, UK, 2005. Morgan-Kaufmann Publishers.360
[2] A. Balke and J. Pearl. Counterfactual probabilities: Computational methods, bounds, and361 applications. In R. L. de Mantaras and D. Poole, editors, Uncertainty in Artificial Intelligence362 10, pages 46–54. Morgan Kaufmann, San Mateo, CA, 1994.363
[3] A. Balke and J. Pearl. Counterfactuals and policy analysis in structural models. In P. Besnard364 and S. Hanks, editors, Uncertainty in Artificial Intelligence 11, pages 11–18. San Francisco,365 1995.366
[4] A. Balke and J. Pearl. Bounds on treatment effects from studies with imperfect compliance.367 Journal of the American Statistical Association, 92(439):1172–1176, September 1997.368
[5] E. Bareinboim and J. Pearl. Causal inference by surrogate experiments: z-identifiability.369 In N. de Freitas and K. Murphy, editors, Proceedings of the Twenty-Eighth Conference on370 Uncertainty in Artificial Intelligence, pages 113–120, Corvallis, OR, 2012. AUAI Press.371
[6] H. Bauer. Probability theory and elements of measure theory. Holt, 1972.372
[7] H. Bauer. Measure and integration theory, volume 26. Walter de Gruyter, 2011.373
[8] F. A. Bugni. Bootstrap inference in partially identified models defined by moment inequalities:374 Coverage of the identified set. Econometrica, 78(2):735–753, 2010.375
[9] C. Carathéodory. Über den variabilitätsbereich der fourier’schen konstanten von positiven har-376 monischen funktionen. Rendiconti Del Circolo Matematico di Palermo (1884-1940), 32(1):193–377 217, 1911.378
[10] A. Carolei et al. The international stroke trial (ist): a randomized trial of aspirin, subcutaneous379 heparin, both, or neither among 19435 patients with acute ischaemic stroke. The Lancet,380 349:1569–1581, 1997.381
[11] D. Chickering and J. Pearl. A clinician’s tool for analyzing non-compliance. Computing Science382 and Statistics, 29(2):424–431, 1997.383
[12] R. J. Connor and J. E. Mosimann. Concepts of independence for proportions with a generaliza-384 tion of the dirichlet distribution. Journal of the American Statistical Association, 64(325):194–385 206, 1969.386
[13] J. Eckhoff. Helly, radon, and carathéodory type theorems. In Handbook of convex geometry,387 pages 389–448. Elsevier, 1993.388
[14] R. J. Evans. Graphical methods for inequality constraints in marginalized dags. In 2012 IEEE389 International Workshop on Machine Learning for Signal Processing, pages 1–6. IEEE, 2012.390
[15] R. J. Evans et al. Margins of discrete bayesian networks. The Annals of Statistics, 46(6A):2623–391 2656, 2018.392
[16] N. Finkelstein and I. Shpitser. Deriving bounds and inequality constraints using logical relations393 among counterfactuals. In Conference on Uncertainty in Artificial Intelligence, pages 1348–394 1357. PMLR, 2020.395
[17] C. Frangakis and D. Rubin. Principal stratification in causal inference. Biometrics, 1(58):21–29,396 2002.397
[18] D. Galles and J. Pearl. An axiomatic characterization of causal counterfactuals. Foundation of398 Science, 3(1):151–182, 1998.399
[19] J. Halpern. Axiomatizing causal reasoning. In G. Cooper and S. Moral, editors, Uncertainty400 in Artificial Intelligence, pages 202–210. Morgan Kaufmann, San Francisco, CA, 1998. Also,401 Journal of Artificial Intelligence Research 12:3, 17–37, 2000.402
[20] G. W. Imbens and C. F. Manski. Confidence intervals for partially identified parameters.403 Econometrica, 72(6):1845–1857, 2004.404
[21] G. W. Imbens and D. B. Rubin. Bayesian inference for causal effects in randomized experiments405 with noncompliance. The annals of statistics, pages 305–327, 1997.406
[22] H. Ishwaran and L. F. James. Gibbs sampling methods for stick-breaking priors. Journal of the407 American Statistical Association, 96(453):161–173, 2001.408
[23] N. Kallus and A. Zhou. Confounding-robust policy improvement. In Advances in neural409 information processing systems, pages 9269–9279, 2018.410
[24] N. Kallus and A. Zhou. Confounding-robust policy evaluation in infinite-horizon reinforcement411 learning. Advances in Neural Information Processing Systems, 2020.412
[25] N. Kilbertus, M. J. Kusner, and R. Silva. A class of algorithms for general instrumental variable413 models. In Advances in Neural Information Processing Systems, 2020.414
[26] J. B. Lasserre. Global optimization with polynomials and the problem of moments. SIAM415 Journal on optimization, 11(3):796–817, 2001.416
[27] J. B. Lasserre. Moments, positive polynomials and their applications, volume 1. World Scientific,417 2009.418
[28] M. Laurent. Sums of squares, moment matrices and optimization over polynomials. In Emerging419 applications of algebraic geometry, pages 157–270. Springer, 2009.420
[29] H. R. Lewis. Computers and intractability. a guide to the theory of np-completeness, 1983.421
[30] C. Manski. Nonparametric bounds on treatment effects. American Economic Review, Papers422 and Proceedings, 80:319–323, 1990.423
[31] P. A. Parrilo. Semidefinite programming relaxations for semialgebraic problems. Mathematical424 programming, 96(2):293–320, 2003.425
[32] J. Pearl. Causal diagrams for empirical research. Biometrika, 82(4):669–710, 1995.426
[33] J. Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, New York,427 2000. 2nd edition, 2009.428
[34] J. Pearl. Principal stratification – a goal or a tool? The International Journal429 of Biostatistics, 7(1), 2011. Article 20, DOI: 10.2202/1557-4679.1322. Available at:430 <http://ftp.cs.ucla.edu/pub/stat_ser/r382.pdf>.431
[35] A. Richardson, M. G. Hudgens, P. B. Gilbert, and J. P. Fine. Nonparametric bounds and432 sensitivity analysis of treatment effects. Statistical science: a review journal of the Institute of433 Mathematical Statistics, 29(4):596, 2014.434
[36] J. Robins. The analysis of randomized and non-randomized aids treatment trials using a new435 approach to causal inference in longitudinal studies. In L. Sechrest, H. Freeman, and A. Mulley,436 editors, Health Service Research Methodology: A Focus on AIDS, pages 113–159. NCHSR,437 U.S. Public Health Service, Washington, D.C., 1989.438
[37] J. P. Romano and A. M. Shaikh. Inference for identifiable parameters in partially identified439 econometric models. Journal of Statistical Planning and Inference, 138(9):2786–2807, 2008.440
[38] D. Rosset, N. Gisin, and E. Wolfe. Universal bound on the cardinality of local hidden variables441 in networks. Quantum Information & Computation, 18(11-12):910–926, 2018.442
[39] P. K. Sen and J. M. Singer. Large sample methods in statistics: an introduction with applications,443 volume 25. CRC press, 1994.444
[40] J. Sethuraman. A constructive definition of dirichlet priors. Statistica sinica, pages 639–650,445 1994.446
[41] I. Shpitser and J. Pearl. What counterfactuals can be tested. In Proceedings of the Twenty-Third447 Conference on Uncertainty in Artificial Intelligence, pages 352–359. AUAI Press, Vancouver,448 BC, Canada, 2007. Also, Journal of Machine Learning Research, 9:1941–1979, 2008.449
[42] I. Shpitser and E. Sherman. Identification of personalized effects associated with causal450 pathways. In UAI, 2018.451
[43] J. Tian. Studies in Causal Reasoning and Learning. PhD thesis, Computer Science Department,452 University of California, Los Angeles, CA, November 2002.453
[44] J. Tian and J. Pearl. Probabilities of causation: Bounds and identification. Annals of Mathematics454 and Artificial Intelligence, 28:287–313, 2000.455
[45] J. Tian and J. Pearl. A general identification condition for causal effects. In Proceedings of the456 Eighteenth National Conference on Artificial Intelligence, pages 567–573. AAAI Press/The457 MIT Press, Menlo Park, CA, 2002.458
[46] D. Todem, J. Fine, and L. Peng. A global sensitivity test for evaluating statistical hypotheses459 with nonidentifiable models. Biometrics, 66(2):558–566, 2010.460
[47] S. Vansteelandt, E. Goetghebeur, M. G. Kenward, and G. Molenberghs. Ignorance and uncer-461 tainty regions as inferential tools in a sensitivity analysis. Statistica Sinica, pages 953–979,462 2006.463
[48] H. Waki, S. Kim, M. Kojima, and M. Muramatsu. Sums of squares and semidefinite program464 relaxations for polynomial optimization problems with structured sparsity. SIAM Journal on465 Optimization, 17(1):218–242, 2006.466
[49] J. Zhang and E. Bareinboim. Bounding causal effects on continuous outcomes. In Proceedings467 of the 35nd AAAI Conference on Artificial Intelligence, 2021.468
Checklist469
1. For all authors...470 (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s471 contributions and scope? [Yes]472 (b) Did you describe the limitations of your work? [Yes] “Throughout this paper, we473 assume that endogenous variables V are discrete and finite; while exogenous variables474 U could take any (continuous) value.”475
(c) Did you discuss any potential negative societal impacts of your work? [N/A] This work476 does not present any foreseeable societal consequence.477
(d) Have you read the ethics review guidelines and ensured that your paper conforms to478 them? [Yes]479
2. If you are including theoretical results...480 (a) Did you state the full set of assumptions of all theoretical results? [Yes] See Sec. 1.1.481 (b) Did you include complete proofs of all theoretical results? [Yes] See Appendices A482 and B.483 3. If you ran experiments...484
(a) Did you include the code, data, and instructions needed to reproduce the main experi-485 mental results (either in the supplemental material or as a URL)? [Yes] See Appendix C.486
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they487 were chosen)? [Yes] See Appendix C.488
(c) Did you report error bars (e.g., with respect to the random seed after running experi-489 ments multiple times)? [N/A]490
(d) Did you include the total amount of compute and the type of resources used (e.g., type491 of GPUs, internal cluster, or cloud provider)? [Yes] See Appendix C. “Experiments492 were performed on a computer with 32GB memory.”493
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...494 (a) If your work uses existing assets, did you cite the creators? [Yes] “IST was a large,495 randomized, open trial of up to 14 days of antithrombotic therapy after stroke onset496 [10].” See also Appendix C497
(b) Did you mention the license of the assets? [Yes] See Appendix C. The IST dataset is498 shared under “Open Data Commons Attribution License (ODC-By) v1.0”.499
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]500 501
(d) Did you discuss whether and how consent was obtained from people whose data you’re502 using/curating? [N/A]503
(e) Did you discuss whether the data you are using/curating contains personally identifiable504 information or offensive content? [N/A]505
5. If you used crowdsourcing or conducted research with human subjects...506 (a) Did you include the full text of instructions given to participants and screenshots, if507 applicable? [N/A]508 (b) Did you describe any potential participant risks, with links to Institutional Review509 Board (IRB) approvals, if applicable? [N/A]510 (c) Did you include the estimated hourly wage paid to participants and the total amount511 spent on participant compensation? [N/A]512

Review summaries:
EDIT:  I thank the authors for their clarifications and corrections to my previous line of thought, and I am prepared to revise my score up to 5 as a result.  I take the point about polynomial vs linear programming, and I can see that Lemma 4 is a non-trivial result.  The authors, however, make no effort to point to Lemma 4 as being the primary source of the novelty in Theorem 1!  Overall, the content of the paper seems too rich for a conference with a page limit like this one, it might be better off submitted to a journal with much more explanation and more examples.  ----  This paper gives a Monte Carlo method for obtaining bounds from any DAG model with hidden variables.   The authors illustrate the algorithm on the front door model, as well as the instrumental variables model; they give several more examples in the supplementary material.    I think the paper could be improved by explaining more explicitly how the Bayesian algorithm for obtaining samples helps, when presumably one could solve the linear program (6).

The paper presents a methodology to obtain bounds on counterfactual probabilities from unknown SCMs using observational data. The authors consider a setting where the causal structure (DAG) of the data-generating model is known and all the endogenous variables are discrete and finite.    In this setting, the authors show that the set of all SCMs compatible with the known DAG, G, and observational data can be mapped to a special class of SCMs, which they call discrete SCMs. These discrete SCMs have the nice property that they are compatible with DAG G and the observational distribution and all exogenous variables are discrete and finite (with known maximum value). The authors show that there must exist a discrete SCM within this family which agrees with the original (unknown) SCM on all counterfactual probabilities. Using this as the main idea, the authors attempt to bound the counterfactual probabilities by maximising/minimising the counterfactual probabilities over all discrete SCMs in this family.    The paper makes no assumptions about the missing exogenous variables in the SCM model under question. This makes for a very powerful idea which can be used to bound counterfactual queries.