Paper deep dive
General sample size analysis for probabilities of causation: a delta method approach
Tianyuan Cheng, Ruirui Mao, Judea Pearl, Ang Li
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/21/2026, 12:33:17 AM
Summary
The paper proposes a general sample size framework for estimating bounds of Probabilities of Causation (PoCs), such as PNS, PN, and PS, using the delta method. It addresses the gap in determining required experimental and observational sample sizes to achieve a desired margin of error, demonstrating through simulations that the approach provides stable and less conservative estimates than existing methods.
Entities (12)
Relation Signals (10)
Ruirui Mao → affiliatedwith → Cornell University
confidence 99% · Ruirui Mao Department of Computer Science, Cornell University
Judea Pearl → affiliatedwith → University of California, Los Angeles
confidence 99% · Judea Pearl Department of Computer Science, University of California Los Angeles
Ang Li → affiliatedwith → Florida State University
confidence 99% · Ang Li Department of Computer Science, Florida State University
Tianyuan Cheng → affiliatedwith → Florida State University
confidence 99% · Tianyuan Cheng Department of Statistics, Florida State University
Probabilities of Causation → includes → Probability of Necessity and Sufficiency
confidence 95% · Probabilities of causation (PoCs), such as the probability of necessity and sufficiency (PNS)...
Probabilities of Causation → includes → Probability of Necessity
confidence 95% · Pearl (1999) defined three binary PoCs, including PNS, PN, and PS.
Probabilities of Causation → includes → Probability of Sufficiency
confidence 95% · Pearl (1999) defined three binary PoCs, including PNS, PN, and PS.
Structural Causal Model → defines → Probabilities of Causation
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Probabilities of causation (PoCs), such as the probability of necessity and sufficiency (PNS), are important tools for decision making but are generally not point identifiable. Existing work has derived bounds for these quantities using combinations of experimental and observational data. However, there is very limited research on sample size analysis, namely, how many experimental and observational samples are required to achieve a desired margin of error. In this paper, we propose a general sample size framework based on the delta method. Our approach applies to settings in which the target bounds of PoCs can be expressed as finite minima or maxima of linear combinations of experimental and observational probabilities. Through simulation studies, we demonstrate that the proposed sample size calculations lead to stable estimation of these bounds.
Tags
Links
- Source: https://arxiv.org/abs/2602.17070v1
- Canonical: https://arxiv.org/abs/2602.17070v1
Trouble viewing inline? Open PDF directly →
Full Text
50,888 characters extracted from source content.
Expand or collapse full text
General sample size analysis for probabilities of causation: a delta method approach Tianyuan Cheng Department of Statistics, Florida State University Ruirui Mao Department of Computer Science, Cornell University Judea Pearl Department of Computer Science, University of California Los Angeles Ang Li Department of Computer Science, Florida State University Abstract Probabilities of causation (PoCs), such as the probability of necessity and sufficiency (PNS), are important tools for decision making but are generally not point identifiable. Existing work has derived bounds for these quantities using combinations of experimental and observational data. However, there is very limited research on sample size analysis, namely, how many experimental and observational samples are required to achieve a desired margin of error. In this paper, we propose a general sample size framework based on the delta method. Our approach applies to settings in which the target bounds of PoCs can be expressed as finite minima or maxima of linear combinations of experimental and observational probabilities. Through simulation studies, we demonstrate that the proposed sample size calculations lead to stable estimation of these bounds. 1 Introduction Probabilities of causation (PoCs) are used in many real-world applications, such as marketing, law, social science, and health science, especially when decisions depend on whether an action caused an outcome. For example, Li and Pearl (2022) introduced a “benefit function” that is a linear combination of PoCs and reflects the payoff or cost of selecting an individual with certain features, with the goal of finding those most likely to show a target behavior. Stott et al. (2004) used it in climate event assignment to quantify how much human influence changes the risk of an extreme event. Mueller and Pearl (2023) argued that PoCs can be used for personalized decision making, and Li et al. (2020) found that PoCs can be helpful for improving the accuracy of some machine learning methods. Without extra assumptions, PoCs are generally not identifiable, so one often works with bounds instead of point values. Using the structural causal model (SCM), Pearl (1999) defined three binary PoCs, including PNS, PN, and PS. Tian and Pearl (2000) derived bounds for these quantities using both experimental and observational information, and later Li and Pearl (2019, 2024) provided formal proofs. Several papers studied how to tighten these bounds. For example, Mueller et al. (2021) used covariates and causal structure to narrow the bounds for PNS, and Dawid et al. (2017) used covariates to narrow the bounds for PN. Most of these works implicitly assume that the experimental and observational samples are large enough to estimate the needed probabilities well. However, there is little work on how to choose sample sizes so that the estimated bounds have a desired level of precision. This gap limits the use of the theory in real applications. Li et al. (2022) discussed this issue, but focused on a special case for the PNS bound. To our knowledge, there is still no unified framework that links a target error level to required experimental and observational sample sizes for general PoC bounds. In this paper, we study adequate sample sizes for estimating bounds of PoCs from a general perspective. Our starting point is that many sharp bounds can be written as a finite minimum or maximum of explicit functions of a finite set of probabilities, including experimental probabilities such as P(yx)P(y_x) and observational probabilities such as P(x,y)P(x,y). For PNS, the bound components are typically linear in these probabilities; for PN and PS, the bound components often take a ratio form. Under mild regularity conditions (e.g., denominators bounded away from zero), the bound components are smooth transformations of the underlying probability vector. The resulting lower and upper bounds can still be non-smooth because they are formed by finite minima or maxima. The fact that the components are smooth but the bounds can be non-smooth leads to different asymptotic behavior in the two cases, which our results accommodate. Our main contributions are: • We propose a general sample size framework for estimating PoC bound endpoints (i.e., lower and upper bounds) with a pre-specified margin of error, based on multivariate delta-method variance approximations for smooth endpoints, and a directional delta method implemented via numerical methods for non-smooth endpoints, following Fang and Santos (2019). • Our asymptotic results apply whenever the bound endpoints can be expressed as finite minima or maxima of explicit bound components. This covers common forms in the PoC literature and also extends to other bounded causal quantities, including linear combinations of PoCs. • We provide simulation studies showing that the proposed sample sizes are stable and sufficient in practice, and that they are far less conservative than existing results in the literature. 2 Preliminaries In this section, we briefly review the definitions of the three aspects of binary causation following Tian and Pearl (2000). Our analysis is based on the counterfactual framework within structural causal models (SCMs) as introduced in Pearl (2009). We denote by Yx=yY_x=y the counterfactual statement that variable Y would take value y if X were set to x. Throughout the paper, we use yxy_x to represent the event Yx=yY_x=y, yx′y_x for Yx′=yY_x =y, yx′y_x for Yx=y′Y_x=y , and yx′y_x for Yx′=y′Y_x =y . Experimental information is summarized through causal quantities such as P(yx)P(y_x), while observational information is summarized by joint distributions such as P(x,y)P(x,y). Unless otherwise stated, X denotes the treatment variable and Y denotes the outcome variable. For the following three probabilities of causation, we all assume that X and Y are two binary variables in a causal model M. Let x and y stand for the propositions X=trueX=true and Y=trueY=true, respectively, and x′x and y′y for their complements. Then we have: Definition 1 (Probability of Necessity (PN)). The probability of necessity is defined as: PN≜ \; P(Yx′=false∣X=true,Y=true) P(Y_x =false X=true,Y=true) (1) ≜ \; P(yx′∣x,y). P(y_x x,y). Definition 2 (Probability of Sufficiency (PS)). The probability of sufficiency is defined as: PS≜ \; P(Yx=true∣X=false,Y=false) P(Y_x=true X=false,Y=false) (2) ≜ \; P(yx∣x′,y′). P(y_x x ,y ). Definition 3 (Probability of Necessity and Sufficiency (PNS)). The probability of necessity and sufficiency is defined as: PNS≜ \; P(Yx=true,Yx′=false) P(Y_x=true,\,Y_x =false) (3) ≜ \; P(yx,yx′). P(y_x,y_x ). We also introduce the sharp bounds on the PoCs derived from experimental and observational data in Tian and Pearl (2000), where the bounds are obtained via Balke’s program in Balke (1995). The bounds are: max0,P(yx)−P(yx′),P(y)−P(yx′),P(yx)−P(y)≤PNS≤minP(yx),P(yx′),P(x,y)+P(x′,y′),P(yx)−P(yx′)+P(x,y′)+P(x′,y). \ array[]l0,\\ P(y_x)-P(y_x ),\\ P(y)-P(y_x ),\\ P(y_x)-P(y) array \ \;≤\; \ array[]lP(y_x),\\ P(y_x ),\\ P(x,y)+P(x ,y ),\\ P(y_x)-P(y_x )+P(x,y )+P(x ,y) array \. (4) max0,P(y)−P(yx′)P(x,y)≤PN≤min1,P(yx′)−P(x′,y′)P(x,y). \0,\; P(y)-P(y_x )P(x,y) \ ≤ \1,\; P(y _x )-P(x ,y )P(x,y) \. (5) max0,P(yx)−P(y)P(x′,y′)≤PS≤min1,P(yx)−P(x,y)P(x′,y′). \0,\; P(y_x)-P(y)P(x ,y ) \ ≤ \1,\; P(y_x)-P(x,y)P(x ,y ) \. (6) 3 Main Results Let θ be a parameter vector that collects all experimental and observational probabilities we need. For example, in the setting of the upper bound of (4) we can take θ=(P(yx),P(yx′),P(x,y),P(x′,y′),P(x,y′),P(x′,y))⊤.θ=(P(y_x),P(y_x ),P(x,y),P(x ,y ),P(x,y ),P(x ,y)) . (7) We first claim that in most scenarios, the bounds of PoCs can be written in the form of a complex linear combination of θ. Lemma 1 (Piecewise-linear and fractional forms of sharp bounds for PoCs). Consider a structural causal model with finite discrete variables, and focus on the case that X and Y are both binary, i.e., X has two treatment levels x,x′x,x and Y has two outcome levels y,y′y,y . Let θ∈ℝdθ ^d collect the observational and experimental probabilities used to constrain the model. For a probability of causation Q∈PNS,PN,PSQ∈\PNS,PN,PS\ as defined in Section 22, let UQ(θ)U_Q(θ) and LQ(θ)L_Q(θ) be its sharp upper and lower bounds over all SCMs. Then there exist finite positive integers J,K<∞J,K<∞, vectors aj,ck∈ℝda_j,c_k ^d and constants bj,dk∈ℝb_j,d_k such that UQ(θ)=1hQ(θ)min1≤j≤Jaj⊤θ+bj,LQ(θ)=1hQ(θ)max1≤k≤Kck⊤θ+dk,U_Q(θ)= 1h_Q(θ) _1≤ j≤ J\a_j θ+b_j\,\,\,\,L_Q(θ)= 1h_Q(θ) _1≤ k≤ K\c_k θ+d_k\, (8) where the denominator hQ(θ)h_Q(θ) is defined as: 1. For PNS: hPNS(θ)=1h_PNS(θ)=1 (reducing to a piecewise-linear form). 2. For PN: hPN(θ)=P(x,y)h_PN(θ)=P(x,y), which is a specific component of θ. 3. For PS: hPS(θ)=P(x′,y′)h_PS(θ)=P(x ,y ), which is a specific component of θ. UQ(θ)U_Q(θ) and LQ(θ)L_Q(θ) are continuous and piecewise-defined functions of θ. They are differentiable at any θ where the optimizer is unique. Proof. According to Tian and Pearl (2000), for finite discrete SCMs, sharp bounds for PNS (and the numerators of PN, PS) can be obtained by optimizing a linear functional over the set of latent-type distributions consistent with θ, which is a linear programming (LP) problem. Then, by LP duality and the polyhedral structure of the dual feasible set (see Boyd and Vandenberghe (2004)), the optimal value is a piecewise-linear function of θ and has a finite min or max representation of linear functions a⊤θ+ba θ+b. ∎ Lemma 11 gives a unified form for the sharp bounds: for PNS they are piecewise-linear in θ while for PN and PS they are piecewise linear-fractional due to normalization by P(x,y)P(x,y) and P(x′,y′)P(x ,y ). In either case, the bound endpoints are piecewise smooth and are differentiable whenever the active term is unique (optimizer is unique), and thus the inference for the plug-in estimators UQ^=UQ(θ^) U_Q=U_Q( θ) and LQ^=LQ(θ^) L_Q=L_Q( θ) reduces to studying piecewise-defined functions of θ θ. Next, we establish asymptotic normality and confidence intervals in the regular case where the optimizer is unique. We first state the setup: let m be the experimental sample size and n be the observational sample size. Assume r is the ratio of experimental sample size to observational sample size, i.e., m/n→r∈(0,∞)m/n→ r∈(0,∞). Let θ=(θe⊤,θo⊤)⊤∈ℝdθ=( _e , _o ) ^d collect the experimental and observational probabilities used in Lemma 11, with true value θ0=(θe0⊤,θo0⊤)⊤∈ℝd _0=( _e0 , _o0 ) ^d. Let θ^=(θ^e⊤,θ^o⊤)⊤ θ=( θ_e , θ_o ) be the plug-in estimator based on the two independent samples. By the multivariate central limit theorem (CLT) for multinomial proportions, m(θ^e−θe0)⇒N(0,Σe),n(θ^o−θo0)⇒N(0,Σo), m\,( θ_e- _e0)\ \ N(0, _e), n\,( θ_o- _o0)\ \ N(0, _o), (9) and the two limits are independent. Equivalently, n(θ^−θ0)⇒N(0,Ωr),Ωr=(Σe/r00Σo). n\,( θ- _0)\ \ N(0, _r),\,\, _r= pmatrix _e/r&0\\ 0& _o pmatrix. (10) Theorem 1. Assume hQ(θ0)>0h_Q( _0)>0. Let j∗=argmin1≤j≤Jaj⊤θ0+bj,k∗=argmax1≤k≤Kck⊤θ0+dk,j^*= _1≤ j≤ J\a_j _0+b_j\, k^*= _1≤ k≤ K\c_k _0+d_k\, and assume both optimizers are unique with positive gaps: ΔU:=minj≠j∗(aj⊤θ0+bj)−(aj∗⊤θ0+bj∗)>0, _U:= _j≠ j^*(a_j _0+b_j)-(a_j^* _0+b_j^*)>0, ΔL:=(ck∗⊤θ0+dk∗)−maxk≠k∗(ck⊤θ0+dk)>0. _L:=(c_k^* _0+d_k^*)- _k≠ k^*(c_k _0+d_k)>0. Suppose n(θ^−θ0)⇒N(0,Ωr) n\,( θ- _0)\ \ N(0, _r). Then UQU_Q and LQL_Q are differentiable at θ0 _0 and n(U^Q−UQ(θ0))⇒N(0,∇UQ(θ0)⊤Ωr∇UQ(θ0)), n\,( U_Q-U_Q( _0))\ \ N\! (0,\ ∇ U_Q( _0) _r\,∇ U_Q( _0) ), (11) n(L^Q−LQ(θ0))⇒N(0,∇LQ(θ0)⊤Ωr∇LQ(θ0)). n\,( L_Q-L_Q( _0))\ \ N\! (0,\ ∇ L_Q( _0) _r\,∇ L_Q( _0) ). ∎ Proof. We prove the result for the upper bound UQU_Q only and the proof for lower bound is similar. Define g(θ):=min1≤j≤Jaj⊤θ+bj,UQ(θ)=g(θ)hQ(θ).g(θ):= _1≤ j≤ J\a_j θ+b_j\, U_Q(θ)= g(θ)h_Q(θ). By assumption, the minimizer at θ0 _0 is unique and has a positive gap: ΔU=minj≠j∗(aj⊤θ0+bj)−(aj∗⊤θ0+bj∗)>0. _U= _j≠ j^* (a_j _0+b_j )- (a_j^* _0+b_j^* )>0. Since each aj⊤θ+bja_j θ+b_j is continuous and the index set is finite, the gap implies that there exist a ϵ>0ε>0 such that for all ‖θ−θ0‖≤ϵ\|θ- _0\|≤ε, the minimizer is j∗j^*. Hence g(θ)g(θ) is differentiable at θ0 _0. By Danskin’s theorem (Danskin (2012)), ∇g(θ0)=aj∗.∇ g( _0)=a_j^*. Because hQh_Q is differentiable at θ0 _0 and hQ(θ0)>0h_Q( _0)>0, it follows that UQU_Q is differentiable at θ0 _0 with ∇UQ(θ0)=hQ(θ0)∇g(θ0)−g(θ0)∇hQ(θ0)hQ(θ0)2=hQ(θ0)aj∗−(aj∗⊤θ0+bj∗)∇hQ(θ0)hQ(θ0)2.∇ U_Q( _0)= h_Q( _0)∇ g( _0)-g( _0)∇ h_Q( _0)h_Q( _0)^2= h_Q( _0)a_j^*-(a_j^* _0+b_j^*)∇ h_Q( _0)h_Q( _0)^2. Now apply the multivariate delta method to the differentiable map θ↦UQ(θ)θ U_Q(θ) at θ0 _0: n(UQ(θ^)−UQ(θ0))⇒N(0,∇UQ(θ0)⊤Ωτ∇UQ(θ0)). n (U_Q( θ)-U_Q( _0) )\ \ N (0,\ ∇ U_Q( _0) \, _τ\,∇ U_Q( _0) ). (12) ∎ Using Theorem 1 we can derive a corresponding confidence interval as follows: Corollary 1. Assume the conditions of Theorem 11. Let Ω^r _r be a consistent estimator of Ωr _r. Define the plug-in variance estimators V^U:=∇UQ(θ^)⊤Ω^r∇UQ(θ^),V^L:=∇LQ(θ^)⊤Ω^r∇LQ(θ^), V_U:=∇ U_Q( θ) _r∇ U_Q( θ), V_L:=∇ L_Q( θ) _r∇ L_Q( θ), and the corresponding standard errors SE^(U^Q):=V^U/n,SE^(L^Q):=V^L/n. SE( U_Q):= V_U/n, SE( L_Q):= V_L/n. Then, asymptotic (1−α)(1-α) confidence intervals for UQ(θ0)U_Q( _0) and LQ(θ0)L_Q( _0) are CIU=[U^Q−z1−α/2SE^(U^Q),U^Q+z1−α/2SE^(U^Q)],CI_U= [ U_Q-z_1-α/2\, SE( U_Q),\, U_Q+z_1-α/2\, SE( U_Q) ], (13) CIL=[L^Q−z1−α/2SE^(L^Q),L^Q+z1−α/2SE^(L^Q)],CI_L= [ L_Q-z_1-α/2\, SE( L_Q),\, L_Q+z_1-α/2\, SE( L_Q) ], (14) where z1−α/2z_1-α/2 can be found on z-table of standard normal distribution.∎ Proof. By Theorem 1, we have n(U^Q−UQ(θ0))⇒N(0,VU),VU:=∇UQ(θ0)⊤Ωr∇UQ(θ0), n ( U_Q-U_Q( _0) ) N(0,V_U), V_U:=∇ U_Q( _0) _r∇ U_Q( _0), and similarly n(L^Q−LQ(θ0))⇒N(0,VL),VL:=∇LQ(θ0)⊤Ωr∇LQ(θ0). n ( L_Q-L_Q( _0) ) N(0,V_L), V_L:=∇ L_Q( _0) _r∇ L_Q( _0). Let Ω^r _r be a consistent estimator of Ωr _r. Under the uniqueness conditions in Theorem 1, UQU_Q and LQL_Q are differentiable in a neighborhood of θ0 _0, hence ∇UQ(⋅)∇ U_Q(·) and ∇LQ(⋅)∇ L_Q(·) are continuous at θ0 _0. Since θ^→θ0 θ p _0, we have ∇UQ(θ^)→∇UQ(θ0),∇LQ(θ^)→∇LQ(θ0).∇ U_Q( θ) p∇ U_Q( _0), ∇ L_Q( θ) p∇ L_Q( _0). Therefore, by Slutsky’s theorem (Van der Vaart (2000)), V^U:=∇UQ(θ^)⊤Ω^r∇UQ(θ^)→VU,V^L:=∇LQ(θ^)⊤Ω^r∇LQ(θ^)→VL. V_U:=∇ U_Q( θ) _r∇ U_Q( θ) pV_U, V_L:=∇ L_Q( θ) _r∇ L_Q( θ) pV_L. Define SE^(U^Q)=V^U/n SE( U_Q)= V_U/n and SE^(L^Q)=V^L/n SE( L_Q)= V_L/n. Then U^Q−UQ(θ0)SE^(U^Q)⇒N(0,1),L^Q−LQ(θ0)SE^(L^Q)⇒N(0,1), U_Q-U_Q( _0) SE( U_Q) N(0,1), L_Q-L_Q( _0) SE( L_Q) N(0,1), (15) which gives an asymptotic (1−α)(1-α) confidence intervals. ∎ When the unique-optimizer assumption in Theorem 11 fails, means that when we have multiple active constraints, then the bound function may be non-differentiable at θ0 _0. In this case we use a directional delta method (see Dümbgen (1993), Fang and Santos (2019)) to derive the limiting distribution. The following theorem states the result. Theorem 2. Let UQ(θ)U_Q(θ) and LQ(θ)L_Q(θ) be defined as in equation (8). Assume hQ(θ0)>0h_Q( _0)>0 and hQh_Q is differentiable at θ0 _0. Define the active sets as follows: J0=argminj(aj⊤θ0+bj),K0=argmaxk(ck⊤θ0+dk).J_0= _j(\ a_j _0+b_j), K_0= _k(c_k _0+d_k). If n(θ^−θ0)⇒Z∼N(0,Ωr) n ( θ- _0 ) Z N (0, _r ), then n(U^Q−UQ(θ0))⇒minj∈J0gU,j⊤Z,n(L^Q−LQ(θ0))⇒maxk∈K0gL,k⊤Z, n ( U_Q-U_Q( _0) ) _j∈ J_0g_U,j Z, n ( L_Q-L_Q( _0) ) _k∈ K_0g_L,k Z, (16) where the generalized gradients at θ0 _0 are gU,j=ajhQ(θ0)−minℓ(aℓ⊤θ0+bℓ)hQ(θ0)2∇hQ(θ0),j∈J0.g_U,j= a_jh_Q( _0)- _ (a_ _0+b_ )h_Q( _0)^2\,∇ h_Q( _0), j∈ J_0. gL,k=ckhQ(θ0)−maxℓ(cℓ⊤θ0+dℓ)hQ(θ0)2∇hQ(θ0),k∈K0.g_L,k= c_kh_Q( _0)- _ (c_ _0+d_ )h_Q( _0)^2\,∇ h_Q( _0), k∈ K_0. ∎ Remark 1. If |J0|=1|J_0|=1 and |K0|=1|K_0|=1, then the maps UQ(θ)U_Q(θ) and LQ(θ)L_Q(θ) are locally differentiable at θ0 _0 and the asymptotic results reduce to Theorem 1. If |J0|>1|J_0|>1 or |K0|>1|K_0|>1, then in general, the limiting distribution in Theorem 2 is not Gaussian. Instead, it is the distribution of a minimum or maximum of finitely many correlated Gaussian linear forms evaluated at the Gaussian limit Z∼N(0,Ωr)Z N(0, _r). Proof. We prove the result for the upper bound UQU_Q only and the proof for lower bound is similar. For simplicity, denote m(θ)=minj(aj⊤θ+bj)m(θ)= _j(a_j θ+b_j) and r(θ)=1/hQ(θ)r(θ)=1/h_Q(θ) so that U(θ)=r(θ)m(θ)U(θ)=r(θ)m(θ). Let J0=argminj(aj⊤θ0+bj)J_0= _j(a_j _0+b_j). Take any sequence tn↓0t_n 0 and any hn→h_n→ h, and define θn=θ0+tnhn _n= _0+t_nh_n. Let δj=(aj⊤θ0+bj)−m(θ0)≥0 _j=(a_j _0+b_j)-m( _0)≥ 0. Then, m(θn)−m(θ0)tn=minj(aj⊤hn+δjtn). m( _n)-m( _0)t_n= _j (a_j h_n+ _jt_n ). If j∉J0j∉ J_0, then δj>0 _j>0. Since tn↓0t_n 0, we have δj/tn→∞ _j/t_n→∞, and thus aj⊤hn+δj/tn→∞a_j h_n+ _j/t_n→∞. Therefore, for all sufficiently large n, the minimum is attained within J0J_0, i.e. minj(aj⊤hn+δjtn)=minj∈J0(aj⊤hn+δjtn). _j (a_j h_n+ _jt_n )= _j∈ J_0 (a_j h_n+ _jt_n ). Thus, limn→∞m(θn)−m(θ0)tn=minj∈J0aj⊤h=:mθ0′(h). _n→∞ m( _n)-m( _0)t_n= _j∈ J_0a_j h=:m _ _0(h). By assumption, hQh_Q is differentiable at θ0 _0 with hQ(θ0)>0h_Q( _0)>0, so the map r(θ)=1/hQ(θ)r(θ)=1/h_Q(θ) is differentiable at θ0 _0 and rθ0′(h)=−∇hQ(θ0)⊤hhQ(θ0)2.r _ _0(h)=- ∇ h_Q( _0) hh_Q( _0)^2. Using U(θ)=r(θ)m(θ)U(θ)=r(θ)m(θ) and the decomposition U(θn)−U(θ0)tn= U( _n)-U( _0)t_n= r(θn)m(θn)−m(θ0)tn+m(θ0)r(θn)−r(θ0)tn r( _n) m( _n)-m( _0)t_n+m( _0) r( _n)-r( _0)t_n +(r(θn)−r(θ0))(m(θn)−m(θ0))tn, + (r( _n)-r( _0) ) (m( _n)-m( _0) )t_n, the last term is o(1)o(1) because r(θn)−r(θ0)=O(tn)r( _n)-r( _0)=O(t_n) and m(θn)−m(θ0)=O(tn)m( _n)-m( _0)=O(t_n). Thus, U is Hadamard directionally differentiable at θ0 _0 with Uθ0′(h)=r(θ0)mθ0′(h)+m(θ0)rθ0′(h)=1hQ(θ0)minj∈J0aj⊤h−m(θ0)hQ(θ0)2∇hQ(θ0)⊤h.U _ _0(h)=r( _0)\,m _ _0(h)+m( _0)\,r _ _0(h)= 1h_Q( _0) _j∈ J_0a_j h- m( _0)h_Q( _0)^2∇ h_Q( _0) h. (17) In summary, Uθ0′(h)=minj∈J0gU,j⊤h,gU,j=ajhQ(θ0)−m(θ0)hQ(θ0)2∇hQ(θ0),j∈J0.U _ _0(h)= _j∈ J_0g_U,j h, g_U,j= a_jh_Q( _0)- m( _0)h_Q( _0)^2∇ h_Q( _0), j∈ J_0. (18) By the directional delta method (Dümbgen (1993)), n(U^−U(θ0))⇒Uθ0′(Z)=minj∈J0gU,j⊤Z. n ( U-U( _0) ) U _ _0(Z)= _j∈ J_0g_U,j Z. (19) ∎ Theorem 2 gives the asymptotic result of U^Q,L^Q U_Q,\, L_Q, but it is not directly usable in practice because it depends on the unknown active set J0J_0 and K0K_0. A natural idea is to use the standard bootstrap, but this fails due to the discontinuity of the directional derivative map with respect to θ0 _0. To address this, we apply the numerical delta method of Fang and Santos (2019), which consistently approximates the limit law by numerically evaluating directional derivatives via finite differences without needing to identify J0J_0 or K0K_0. The corresponding confidence intervals are stated in the following corollary: Corollary 2. Assume the conditions of Theorem 2. Let Ω^r _r be a consistent estimator of Ωr _r, and let ϵn _n be a sequence satisfying ϵn→0 _n→ 0 and nϵn→∞ n\, _n→∞. Define the numerical directional derivative statistics for B simulated draws Z(b)∼iidN(0,Ω^r)Z^(b) iid N(0, _r), b=1,…,Bb=1,…,B: TU(b)=UQ(θ^+ϵnZ(b))−UQ(θ^)ϵn,TL(b)=LQ(θ^+ϵnZ(b))−LQ(θ^)ϵn.T_U^(b)= U_Q\! ( θ+ _nZ^(b) )-U_Q( θ) _n, T_L^(b)= L_Q\! ( θ+ _nZ^(b) )-L_Q( θ) _n. Let qU(τ)q_U(τ) and qL(τ)q_L(τ) denote the empirical τ-quantiles of TU(b)\T_U^(b)\ and TL(b)\T_L^(b)\ respectively. Then, asymptotic (1−α)(1-α) confidence intervals for UQ(θ0)U_Q( _0) and LQ(θ0)L_Q( _0) are CIU=[U^Q−qU(1−α/2)n,U^Q−qU(α/2)n].CI_U= [ U_Q- q_U(1-α/2) n,\ U_Q- q_U(α/2) n ]. (20) CIL=[L^Q−qL(1−α/2)n,L^Q−qL(α/2)n].CI_L= [ L_Q- q_L(1-α/2) n,\ L_Q- q_L(α/2) n ]. (21) ∎ Remark 2. The representation in Lemma 1 is not specific to PoCs. The same form holds for any causal quantity whose sharp bounds can be written as in the form of equation (8) with h(θ)>0h(θ)>0 and finite J,KJ,K. All asymptotic results in Theorems 1 and 2 continue to hold under the same setup and assumptions. 4 Simulation In this section we run experiments to show that the proposed number of samples provided by Theorems and Corollary in last section are adequate to obtain PoCs within the desired margin errors. We study the bounds for PNS, which is one of the most frequently used PoCs in causal studies. For comparison, we use the same simulation setup as Li et al. (2022), which studies the sharp PNS bound given in Section 2. It is trivial to verify that this sharp bound has the form shown in Lemma 11. 4.1 Experiment: estimating PNS bound Following Li et al. (2022), we consider two parameter settings that share the same causal graph in Figure 1 but differ in their coefficients. Here X is a binary treatment with x=1x=1, x′=0x =0, Y is a binary outcome with y=1y=1, y′=0y =0, and Z=(Z1,…,Z20)Z=(Z_1,…,Z_20) is a set of 20 mutually independent binary covariates. Detailed coefficients are provided in Supplementary materials. Figure 1: The Causal model for PNS bound experiment. X and Y are binary, and Z is a set of 20 independent binary confounders. Let UZiU_Z_i, UXU_X, and UYU_Y be independent Bernoulli exogenous variables. The endogenous covariates are defined as follows: Zi=UZi,i=1,…,20,Z_i=U_Z_i, i=1,…,20, and MX=∑i=120aiZi,MY=∑i=120biZi,M_X= _i=1^20a_iZ_i, M_Y= _i=1^20b_iZ_i, where ai\a_i\ and bi\b_i\ are drawn from Uniform(−1,1).(-1,1). Here are the structural equations and we generate (X,Y)(X,Y) following these equations with a constant C∼Uniform(−1,1)C (-1,1): X=fX(MX,UX)= X=f_X(M_X,U_X)= 1,if MX+UX>0.5,0,otherwise, Y=fY(X,MY,UY)= Y=f_Y(X,M_Y,U_Y)= 1,0<CX+MY+UY<1,or 1<CX+MY+UY<2,0,otherwise. Since all exogenous variables are binary and independent, the population experimental quantities (e.g., P(yx)P(y_x) and P(yx′)P(y_x )) and the population observational quantities (e.g., P(x,y′)P(x,y ) and P(x′,y)P(x ,y)) can be computed exactly by summing over all configurations of (UX,UY,UZ1,…,UZ20)(U_X,U_Y,U_Z_1,…,U_Z_20). For instance, we can compute P(yx)=∑uX,uY,ZP(uX)P(uY)∏i=120P(uZi)fY(1,MY(Z),uY),P(y_x)= _u_X,u_Y,u_ZP(u_X)\,P(u_Y)\, _i=1^20P(u_Z_i)f_Y\! (1,\,M_Y(u_Z),\,u_Y ), and P(x,y)=∑uX,uY,ZP(uX)P(uY)∏i=120[P(uZi)⋅fX(MX(Z),uX)⋅fY(fX,MY(Z),uY)],P(x,y)= _u_X,u_Y,u_ZP(u_X)\,P(u_Y)\, _i=1^20 [P(u_Z_i)· f_X\! (M_X(u_Z),u_X )\,· f_Y\! (f_X,\,M_Y(u_Z),\,u_Y ) ], where Z=(uZ1,…,uZ20)u_Z=(u_Z_1,…,u_Z_20) and MX(Z),MY(Z)M_X(u_Z),M_Y(u_Z) denote the corresponding realizations of MXM_X and MYM_Y. We then generate finite samples under two rules. For the experimental sample of size m, we draw (UX,UY,UZ1,…,UZ20)(U_X,U_Y,U_Z_1,…,U_Z_20) from their Bernoulli distributions, assign X∼Bernoulli(0.5)X (0.5) independently of the exogenous variables, and set Y=fY(X,MY,UY)Y=f_Y(X,M_Y,U_Y). This will give i.i.d. experimental pairs (X,Y)(X,Y). Let mabm_ab be the number of experimental samples with (X=a,Y=b)(X=a,Y=b) and ma⋅=ma0+ma1m_a·=m_a0+m_a1. We estimate P^(yx)=P^(Y=1∣X=1)=m11m1⋅,P^(yx′)=P^(Y=1∣X=0)=m01m0⋅. P(y_x)= P(Y=1 X=1)= m_11m_1·, P(y_x )= P(Y=1 X=0)= m_01m_0·. For the observational sample of size n, we again draw the exogenous variables, set X=fX(MX,UX)X=f_X(M_X,U_X) and Y=fY(X,MY,UY)Y=f_Y(X,M_Y,U_Y), and estimate the required joint probabilities by empirical frequencies. Let nabn_ab be the number of observational samples with (X=a,Y=b)(X=a,Y=b) and n=∑a∈0,1∑b∈0,1nabn= _a∈\0,1\ _b∈\0,1\n_ab. For example, P^(x,y)=n11n,P^(x,y′)=n10n,P^(x′,y)=n01n. P(x,y)= n_11n, P(x,y )= n_10n, P(x ,y)= n_01n. These estimated probabilities are then plugged into the bound formulas to get finite-sample bound estimates. Assume the ratio of experimental data to observational data is m/n→r=1m/n→ r=1. Then, Corollary 1 gives adequate experimental size m≥1921m≥ 1921, while result provided in Li et al. (2022) gives m≥6147m≥ 6147. Our method requires only about 31%31\% of the sample size suggested by the existing result. The detailed computation is in Appendix. Figure 2: Scatter plots under different sample sizes for Model 1. Left to right: n=120,481,1921n=120,481,1921. For both Model 1 and 2, we run 1000 replications, and in each replication, we generate experimental and observational samples with sizes n∈120,481,1921n∈\120,481,1921\. The results are shown in Figures 2 and 3. When n=120n=120, both bounds are highly noisy with large dispersion across replications. Increasing the sample size to n=481n=481 reduces the noise, but the variance remains relatively large. When n=1921n=1921, the adequate size we computed using Corollary 11, we can see that the estimates become very stable, with both the upper and lower bound concentrates tightly around the true bounds, which shows a good accuracy and stability. Figure 3: Scatter plots under different sample sizes for Model 2. Top to bottom: n=120,481,1921n=120,481,1921. We also examine the relationship between sample size n and the average estimation error, defined as (bound^−truebound)( bound-true\ bound) averaged over 1000 replications for each model. See Figure 4. Using the target error level 0.050.05(red horizontal line), both Model 1 and Model 2 reach the desired accuracy with sample sizes below 300300 for both the upper and lower bounds (blue/orange curves). The average error drops quickly as n increases up to around 500500, which suggests that n=500n=500 can already work well in practice. This also confirms that the theoretical sample size n=1921n=1921 is more than sufficient and it represents the adequate size that considers the worst case. Figure 4: Average error of estimations VS sample size. Left: Model 1, Right:Model 2 We further investigate the relationship between the average estimation error and the sample size n beyond Models 1 and 2. In order to approximate a more general setting, we repeatedly resample the SCM coefficients and rerun the same simulation design as in Model 11 for 2020 independent draws. For each n, we then average the estimation error over both the 1000 Monte Carlo replications, and report the resulting curve in Figure 5. It shows that both the upper and lower bound estimation errors decrease rapidly as n increases, and in this aggregated experiment the errors converge even faster than in Models 1 and 2. The target error level 0.050.05 is reached around 100100, smaller than 500500, which is consistent with our earlier observations. The theoretical sample size n=1921n=1921 remains clearly sufficient with very small errors for both bounds. In summary, these results suggest that the sample-size selection rule using our Theorems is not tied to a single parameter setting and remains stable and accurate across a range of SCM with different coefficients. Figure 5: Average error with sample size for 20 replicates 5 Discussion The bounds of PoCs are piecewise-defined, and points with multiple optimizers lie on the boundaries between pieces. Away from these boundaries, the optimizer is unique and thus Theorem 11 applies. Since the measure of boundary points are almost zero, we can think Theorem 11 holds for most of time. This is supported by our simulations, where the target error is often achieved with far fewer samples than the theoretical results (around 1/41/4 to 1/31/3), suggesting that the active piece is quite stable in a neighborhood of θ0 _0. For the non-unique optimizer case in Theorem 22, one may consider a smooth approximation (e.g., Softmin) to enforce differentiability, but this introduces approximation bias and requires hyperparameter tuning. For this reason we adopt a conservative variance-based approach. However, smoothing may still be useful for large-scale repeated evaluations or gradient-based optimization when near-ties are common and numerical stability is a concern. Most importantly, although we use PoCs (PNS, PN, PS) in binary setup as motivating examples, Theorems 11 and 22 apply more broadly. We discuss two extensions here. First, our proposed method can be directly extended to multi-level discrete X and Y, with the same CLT and (directional) delta-method arguments but at the cost of higher dimension and computational difficulty. Second, the same methods can be applied to any bounded causal quantity whose sharp bounds can be written as a finite minimum or maximum of linear (or ratio-type) functions of θ, including linear combinations of PoCs such as the “benefit function” in Li and Pearl (2019). 6 Conclusion We study the adequate sample size for estimating bounds of PoCs within a desired margin of error. In finite discrete settings, we first show that the bounds can be written as finite minima or maxima of explicit functions of observational and experimental probabilities, including both linear and ratio-type terms. Using this representation, when the bound endpoint is smooth at the true parameter, we use multivariate delta-method variance approximations to obtain asymptotic margins of error; when the endpoint is not smooth due to active-set switches in the finite min–max structure, we use a directional delta method and implement it using numerical methods following Fang and Santos (2019). Our asymptotic results apply whenever the bound endpoints has this finite min max form, which covers many bounds in related literature and also extends to other bounded causal quantities, including linear combinations of PoCs. We also provide simulation studies showing that the proposed sample sizes are stable and work well in practice, and are less conservative than existing approaches. To our knowledge, this is the first work providing a general and systematic sample size framework for estimating non-identifiable PoC bounds that handles both smooth and non-smooth endpoints. References A. A. Balke (1995) Probabilistic counterfactuals: semantics, computation, and applications. University of California, Los Angeles. Cited by: §2. S. Boyd and L. Vandenberghe (2004) Convex optimization. Cambridge university press. Cited by: §3. J. M. Danskin (2012) The theory of max-min and its application to weapons allocation problems. Vol. 5, Springer Science & Business Media. Cited by: §3. A. P. Dawid, M. Musio, and R. Murtas (2017) The probability of causation. Law, Probability and Risk 16 (4), p. 163–179. Cited by: §1. L. Dümbgen (1993) On nondifferentiable functions and the bootstrap. Probability Theory and Related Fields 95 (1), p. 125–140. Cited by: §3, §3. Z. Fang and A. Santos (2019) Inference on directionally differentiable functions. The Review of Economic Studies 86 (1), p. 377–412. Cited by: 1st item, §3, §3, §6. A. Li, S. J. Chen, J. Qin, and Z. Qin (2020) Training machine learning models with causal logic. In Companion Proceedings of the Web Conference 2020, p. 557–561. Cited by: §1. A. Li, R. Mao, and J. Pearl (2022) Probabilities of causation: adequate size of experimental and observational samples. arXiv preprint arXiv:2210.05027. Cited by: §1, §4.1, §4.1, §4. A. Li and J. Pearl (2019) Unit selection based on counterfactual logic. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, Cited by: §1, §5. A. Li and J. Pearl (2022) Unit selection with causal diagram. In Proceedings of the AAAI conference on artificial intelligence, Vol. 36, p. 5765–5772. Cited by: §1. A. Li and J. Pearl (2024) Probabilities of causation with nonbinary treatment and effect. In Proceedings of the AAAI conference on artificial intelligence, Vol. 38, p. 20465–20472. Cited by: §1. S. Mueller, A. Li, and J. Pearl (2021) Causes of effects: learning individual responses from population data. arXiv preprint arXiv:2104.13730. Cited by: §1. S. Mueller and J. Pearl (2023) Personalized decision making–a conceptual introduction. Journal of Causal Inference 11 (1), p. 20220050. Cited by: §1. J. Pearl (1999) Probabilities of causation: three counterfactual interpretations and their identification. Synthese 121 (1-2), p. 93. Cited by: §1. J. Pearl (2009) Causality. Cambridge university press. Cited by: §2. P. A. Stott, D. A. Stone, and M. R. Allen (2004) Human contribution to the european heatwave of 2003. Nature 432 (7017), p. 610–614. Cited by: §1. J. Tian and J. Pearl (2000) Probabilities of causation: bounds and identification. Annals of Mathematics and Artificial Intelligence 28 (1), p. 287–313. Cited by: §1, §2, §2, §3. A. W. Van der Vaart (2000) Asymptotic statistics. Vol. 3, Cambridge university press. Cited by: §3. Appendix A Appendix A.1 Calculation of sample size in Experiment 1, PNS bound We will prove the case of PNS upper bound in (4) only. First, focus on the term P(yx)−P(yx′)+P(x,y′)+P(x′,y)P(y_x)-P(y_x )+P(x,y )+P(x ,y). Denote A=P(yx),B=P(yx′),C=P(x,y′),D=P(x′,y)A=P(y_x),B=P(y_x ),C=P(x,y ),D=P(x ,y), and we can see the A and B come from randomized controlled trial, which can be assumed to be independent with each other, and independent with C and D; however, C and D come from a multinomial distribution, (nx,y,nx′,y=nD,nx,y′=nC,nx′,y′)∼Multinomial(n;P(x,y),P(x′,y),P(x,y′),P(x′,y′)),(n_x,y,n_x ,y=n_D,n_x,y =n_C,n_x ,y )\, \,Multinomial(n;P(x,y),P(x ,y),P(x,y ),P(x ,y )), such that C and D are negatively correlated. Combine C and D into one event to make C∪D∼Binomial(n;P(x,y′)+P(x′,y)=:PC∪D)C∪ D (n;P(x,y )+P(x ,y)=:P_C∪ D). Then we have: var(A^−B^+C^+D^)= ( A- B+ C+ D)= var(A^)+var(B^)+var(C+D) ( A)+var( B)+var (C+D) = = P(yx)(1−P(yx))mA+P(yx′(1−P(yx′)mB+PC∪D(1−PC∪D)n. P(y_x)(1-P(y_x))m_A+ P(y_x (1-P(y_x )m_B+ P_C∪ D(1-P_C∪ D)n. By AM-GM inequality, p(1−p)≤1/4p(1-p)≤ 1/4 for any p∈[0,1]p∈[0,1], thus, var(A^−B^+C^+D^)≤14mA+14mB+14n,var( A- B+ C+ D)≤ 14m_A+ 14m_B+ 14n, with constraint mA+mB=m_A+m_B=m. Assume that mA=mB=m/2m_A=m_B=m/2, then Err(A^−B^+C^+D^)= Err( A- B+ C+ D)= z1−α/2var(A^−B^+C^+D^) z_1-α/2 var( A- B+ C+ D) (22) ≤ ≤ z1−α/221mA+1mB+1n z_1-α/22 1m_A+ 1m_B+ 1n = = z1−α/21m+14n. z_1-α/2 1m+ 14n. Taking the error bound to be ϵ=0.05ε=0.05, then (22) gives that m≥(1+14r)(z1−α/2ϵ)2≈1537(1+14r).m≥(1+ 14r)( z_1-α/2ε)^2≈ 1537(1+ 14r). (23) When we take m=nm=n, i.e. r=1r=1, (22) gives adequate experimental size n=m≥1921n=m≥ 1921, which is the sample size in the worst case. A.2 Coefficients of Experiments A.2.1 Model 1 Zi=UZi,i∈1,…,20.Z_i=U_Z_i, i∈\1,…,20\. X=fX(MX,UX)=1,if MX+UX>0.5,0,otherwise.X=f_X(M_X,U_X)= cases1,&if M_X+U_X>0.5,\\ 0,&otherwise. cases Y=fY(X,MY,UY)=1,if 0<CX+MY+UY<1or 1<CX+MY+UY<2,0,otherwise.Y=f_Y(X,M_Y,U_Y)= cases1,&if 0<CX+M_Y+U_Y<1\ or\ 1<CX+M_Y+U_Y<2,\\ 0,&otherwise. cases UZ1 U_Z_1 ∼Bernoulli(0.352913861526), (352913861526), UZ2 U_Z_2 ∼Bernoulli(0.460995855543), (460995855543), UZ3 U_Z_3 ∼Bernoulli(0.331702473392), (331702473392), UZ4 U_Z_4 ∼Bernoulli(0.885505026779), (885505026779), UZ5 U_Z_5 ∼Bernoulli(0.017026872706), (017026872706), UZ6 U_Z_6 ∼Bernoulli(0.380772701708), (380772701708), UZ7 U_Z_7 ∼Bernoulli(0.028092602705), (028092602705), UZ8 U_Z_8 ∼Bernoulli(0.220819399962), (220819399962), UZ9 U_Z_9 ∼Bernoulli(0.617742227477), (617742227477), UZ10 U_Z_10 ∼Bernoulli(0.981975046713), (981975046713), UZ11 U_Z_11 ∼Bernoulli(0.142042291381), (142042291381), UZ12 U_Z_12 ∼Bernoulli(0.833602592350), (833602592350), UZ13 U_Z_13 ∼Bernoulli(0.882938907115), (882938907115), UZ14 U_Z_14 ∼Bernoulli(0.542143191999), (542143191999), UZ15 U_Z_15 ∼Bernoulli(0.085023436884), (085023436884), UZ16 U_Z_16 ∼Bernoulli(0.645357252864), (645357252864), UZ17 U_Z_17 ∼Bernoulli(0.863787135134), (863787135134), UZ18 U_Z_18 ∼Bernoulli(0.460539711624), (460539711624), UZ19 U_Z_19 ∼Bernoulli(0.314014079207), (314014079207), UZ20 U_Z_20 ∼Bernoulli(0.685879388218), (685879388218), UX U_X ∼Bernoulli(0.601680857267), (601680857267), UY U_Y ∼Bernoulli(0.497668975278). (497668975278). C=−0.77953605542.C=-0.77953605542. MX=[Z1Z2⋯Z20]βX,βX=[0.259223510143−0.658140989167−0.750258317680.1629064624260.652023463285−0.08929395865410.421469107769−0.4431296847660.802624388789−0.2257409784990.7166216317170.0650682260309−0.2206903340260.156355773665−0.50693672491−0.7070602781150.418812816935−0.08221187039860.769299853833−0.511585391002].M_X=[\,Z_1\ Z_2\ ·s\ Z_20\,]\ _X, _X= bmatrix0.259223510143\\ -0.658140989167\\ -0.75025831768\\ 0.162906462426\\ 0.652023463285\\ -0.0892939586541\\ 0.421469107769\\ -0.443129684766\\ 0.802624388789\\ -0.225740978499\\ 0.716621631717\\ 0.0650682260309\\ -0.220690334026\\ 0.156355773665\\ -0.50693672491\\ -0.707060278115\\ 0.418812816935\\ -0.0822118703986\\ 0.769299853833\\ -0.511585391002 bmatrix. MY=[Z1Z2⋯Z20]βY,βY=[−0.7928671119180.7599671361470.554377223690.503970540409−0.5271871446510.3786199880910.2692551963010.6715970435940.3960101422740.3252285766430.6578083275740.8016550239930.0907679484097−0.0713852594543−0.0691046005285−0.222582013343−0.848408031595−0.584285069026−0.3248748317990.625621583197].M_Y=[\,Z_1\ Z_2\ ·s\ Z_20\,]\ _Y, _Y= bmatrix-0.792867111918\\ 0.759967136147\\ 0.55437722369\\ 0.503970540409\\ -0.527187144651\\ 0.378619988091\\ 0.269255196301\\ 0.671597043594\\ 0.396010142274\\ 0.325228576643\\ 0.657808327574\\ 0.801655023993\\ 0.0907679484097\\ -0.0713852594543\\ -0.0691046005285\\ -0.222582013343\\ -0.848408031595\\ -0.584285069026\\ -0.324874831799\\ 0.625621583197 bmatrix. A.2.2 Model 2 Zi=UZi,i∈1,…,20.Z_i=U_Z_i, i∈\1,…,20\. X=fX(MX,UX)=1,if MX+UX>0.5,0,otherwise.X=f_X(M_X,U_X)= cases1,&if M_X+U_X>0.5,\\ 0,&otherwise. cases Y=fY(X,MY,UY)=1,if 0<CX+MY+UY<1or 1<CX+MY+UY<2,0,otherwise.Y=f_Y(X,M_Y,U_Y)= cases1,&if 0<CX+M_Y+U_Y<1\ or\ 1<CX+M_Y+U_Y<2,\\ 0,&otherwise. cases UZ1 U_Z_1 ∼Bernoulli(0.524110233482), (524110233482), UZ2 U_Z_2 ∼Bernoulli(0.689566064108), (689566064108), UZ3 U_Z_3 ∼Bernoulli(0.180145428970), (180145428970), UZ4 U_Z_4 ∼Bernoulli(0.317153536644), (317153536644), UZ5 U_Z_5 ∼Bernoulli(0.046268153873), (046268153873), UZ6 U_Z_6 ∼Bernoulli(0.340145244411), (340145244411), UZ7 U_Z_7 ∼Bernoulli(0.100912238566), (100912238566), UZ8 U_Z_8 ∼Bernoulli(0.772038172066), (772038172066), UZ9 U_Z_9 ∼Bernoulli(0.913108434869), (913108434869), UZ10 U_Z_10 ∼Bernoulli(0.364272290967), (364272290967), UZ11 U_Z_11 ∼Bernoulli(0.063667554704), (063667554704), UZ12 U_Z_12 ∼Bernoulli(0.454839320009), (454839320009), UZ13 U_Z_13 ∼Bernoulli(0.586687215140), (586687215140), UZ14 U_Z_14 ∼Bernoulli(0.018824647595), (018824647595), UZ15 U_Z_15 ∼Bernoulli(0.871017316787), (871017316787), UZ16 U_Z_16 ∼Bernoulli(0.164966968157), (164966968157), UZ17 U_Z_17 ∼Bernoulli(0.578925020078), (578925020078), UZ18 U_Z_18 ∼Bernoulli(0.983082980658), (983082980658), UZ19 U_Z_19 ∼Bernoulli(0.018033993991), (018033993991), UZ20 U_Z_20 ∼Bernoulli(0.074629121266), (074629121266), UX U_X ∼Bernoulli(0.29908139311), (29908139311), UY U_Y ∼Bernoulli(0.9226108109253). (9226108109253). C=0.975140894243.C=0.975140894243. MX=[Z1Z2⋯Z20]βX,βX=[ 0.843870221861 0.178759296447−0.372349746729−0.950904544846−0.439457721339−0.725970103834−0.791203963585−0.843183562918−0.68422616618−0.782051030131−0.434420454146−0.445019418094 0.751698021555−0.185984172192 0.191948271392 0.401334543567 0.331387702568 0.522595634402−0.928734581669 0.203436441511].M_X=[\,Z_1\ Z_2\ ·s\ Z_20\,] _X, _X= bmatrix\ \ 0.843870221861\\ \ \ 0.178759296447\\ -0.372349746729\\ -0.950904544846\\ -0.439457721339\\ -0.725970103834\\ -0.791203963585\\ -0.843183562918\\ -0.68422616618\\ -0.782051030131\\ -0.434420454146\\ -0.445019418094\\ \ \ 0.751698021555\\ -0.185984172192\\ \ \ 0.191948271392\\ \ \ 0.401334543567\\ \ \ 0.331387702568\\ \ \ 0.522595634402\\ -0.928734581669\\ \ \ 0.203436441511 bmatrix. MY=[Z1Z2⋯Z20]βY,βY=[−0.453251661832 0.424563325534 0.0924810605305 0.312680246141 0.7676961338 0.124337421843−0.435341306455 0.248957751703−0.161303883519−0.537653062121−0.222087991408 0.190167775134−0.788147770713−0.593030174012−0.308066297974 0.218776507777−0.751253645088−0.11151455376 0.785227235182−0.568046522383].M_Y=[\,Z_1\ Z_2\ ·s\ Z_20\,] _Y, _Y= bmatrix-0.453251661832\\ \ \ 0.424563325534\\ \ \ 0.0924810605305\\ \ \ 0.312680246141\\ \ \ 0.7676961338\\ \ \ 0.124337421843\\ -0.435341306455\\ \ \ 0.248957751703\\ -0.161303883519\\ -0.537653062121\\ -0.222087991408\\ \ \ 0.190167775134\\ -0.788147770713\\ -0.593030174012\\ -0.308066297974\\ \ \ 0.218776507777\\ -0.751253645088\\ -0.11151455376\\ \ \ 0.785227235182\\ -0.568046522383 bmatrix.