Paper deep dive
A Generalized-Bayes Perspective on Counterfactual Explanations: Posterior-Based Decision-Making and Evaluation
Keita Kinjo
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/3/2026, 2:36:22 AM
Summary
This paper establishes a theoretical equivalence between distance-minimization-based counterfactual explanations (CEs) and the maximum a posteriori (MAP) estimate of a Gibbs posterior within the generalized Bayes framework. It introduces the Distance-Prior Generalized Bayes CE (DP-GBCE) and extends this to a unified framework supporting multiple decision rules, including Bayes decisions minimizing expected loss and CVaR-CE for risk aversion. The authors also propose handling model multiplicity by mixing posterior distributions of multiple models using Bayesian weights, and define metrics for evaluating both individual CEs and the posterior distribution.
Entities (8)
Relation Signals (5)
Distance-minimization-based CE → isequivalentto → MAP estimate of Gibbs posterior
confidence 95% · We show that a distance-minimization-based CE is mathematically equivalent to the maximum a posteriori (MAP) estimate of a Gibbs posterior within the generalized Bayes framework
DP-GBCE → isformulationof → Generalized Bayes CE
confidence 94% · We call this formulation the Distance-Prior Generalized Bayes CE (DP-GBCE). ... We refer to this formulation as the Generalized Bayes CE (GBCE).
CVaR-CE → istypeof → Decision Rule
confidence 93% · CVaR-CE, a risk-averse decision rule.
Generalized Bayes → extends → Conventional Bayesian Inference
confidence 92% · Generalized Bayes extends conventional Bayesian inference by allowing a posterior distribution to be constructed from a loss function even when a likelihood is unavailable
Bayesian Model Weights → usedfor → Mixing Posterior Distributions
confidence 90% · uses Bayesian model weights to mix the posterior distributions of multiple models, thereby accounting for model multiplicity
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Counterfactual explanations (CEs) enhance the interpretability of machine learning models by identifying the smallest change to an input required to obtain a desired output. Although CEs are conventionally formulated as a distance-minimization problem, the theoretical basis of this formulation has received limited attention. We show that a distance-minimization-based CE is mathematically equivalent to the maximum a posteriori (MAP) estimate of a Gibbs posterior within the generalized Bayes framework, specifically when a distance-based prior is used. We call this formulation the Distance-Prior Generalized Bayes CE (DP-GBCE). Building on this posterior perspective, we introduce two decision rules beyond MAP within a unified framework: a Bayes decision that minimizes expected decision loss and CVaR-CE, a risk-averse decision rule. We also propose an extension that uses Bayesian model weights to mix the posterior distributions of multiple models, thereby accounting for model multiplicity, where several models have comparable predictive performance. Finally, we define metrics for evaluating both individual CEs and the posterior distribution as a whole, and use experiments on simulated data and Google Trends data to quantify the trade-offs among the decision rules.
Tags
Links
- Source: https://arxiv.org/abs/2607.29077v1
- Canonical: https://arxiv.org/abs/2607.29077v1
Trouble viewing inline? Open PDF directly →
Full Text
69,639 characters extracted from source content.
Expand or collapse full text
A Generalized-Bayes Perspective on Counterfactual Explanations: Posterior-Based Decision-Making and Evaluation Keita Kinjo Faculty of Business, Kyoritsu Women’s University Abstract Counterfactual explanations (CEs) enhance the interpretability of machine learning models by identifying the smallest change to an input required to obtain a desired output. Although CEs are conventionally formulated as a distance-minimization problem, the theoretical basis of this formulation has received limited attention. We show that a distance-minimization-based CE is mathematically equivalent to the maximum a posteriori (MAP) estimate of a Gibbs posterior within the generalized Bayes framework, specifically when a distance-based prior is used. We call this formulation the Distance-Prior Generalized Bayes CE (DP-GBCE). Building on this posterior perspective, we introduce two decision rules beyond MAP within a unified framework: a Bayes decision that minimizes expected decision loss and CVaR-CE, a risk-averse decision rule. We also propose an extension that uses Bayesian model weights to mix the posterior distributions of multiple models, thereby accounting for model multiplicity, where several models have comparable predictive performance. Finally, we define metrics for evaluating both individual CEs and the posterior distribution as a whole, and use experiments on simulated data and Google Trends data to quantify the trade-offs among the decision rules. 1 Introduction Counterfactual explanations have attracted considerable attention as a means of addressing the interpretability challenge in machine learning [3, 14, 5]. Given a trained black-box model, a counterfactual explanation identifies the smallest change to a particular instance’s attributes required to obtain a desired prediction. It thereby indicates which attributes affect the prediction and what actions could change the outcome. For example, if a machine learning model denies an individual insurance coverage, a counterfactual explanation can identify the changes in attributes such as income or occupation that would result in approval. This approach is also known as algorithmic recourse [13, 6], and its connection to adversarial examples has likewise been noted [5, 14, 10]. CE methods typically minimize the cost of changing the original instance into the counterfactual, thereby encouraging similarity between the two [15, 3, 14, 5]. One rationale for this formulation draws on Lewis’s notion of the closest possible world, constructed by making only the minimal changes to the actual world needed to evaluate a causal effect [8]. Despite this philosophical foundation, theoretical discussion has been limited. Although a large required change may intuitively suggest that the corresponding variable is important, the theoretical basis for this interpretation has rarely been examined. Generalized Bayes extends conventional Bayesian inference by allowing a posterior distribution to be constructed from a loss function even when a likelihood is unavailable, as is often the case in machine learning [1]. It constructs a Gibbs posterior through exponential weighting based on the loss rather than on a likelihood. In this study, we show that the CE obtained through distance minimization corresponds to the MAP estimate of a Gibbs posterior over CEs, thereby providing a probabilistic justification for the existing cost-minimization formulation. We further propose a posterior distribution over CEs that incorporates model uncertainty, and we present decision rules other than MAP for the resulting CE posterior. In addition, we propose metrics for evaluating individual CEs and the posterior distribution as a whole, and compare their performance. By clarifying the connection between CEs and Bayesian inference – generalized Bayes in particular – this study opens up alternative approaches to challenges such as robustness, as well as applications that draw on Bayesian insights. Several previous studies have also examined CEs within a Bayesian framework (see Section 4). This study differs from the existing literature in four main respects. First, it establishes the equivalence between CE distance minimization and MAP estimation of a Gibbs posterior, thereby providing existing methods with a probabilistic foundation. Second, it offers a model-agnostic framework that requires neither an explicit generative model nor a specified likelihood and can be applied to any learner – differentiable or otherwise – as long as its loss can be evaluated. This contrasts with [12], which requires a generative probability model and gradient information. Third, it places a probability distribution (the Gibbs posterior) over counterfactual candidates and enables MAP, Bayes, and CVaR-CE decision rules to be compared within a unified framework, allowing robust decision-making under uncertainty in the success region. This contrasts with [9], which seeks a single counterfactual point. Fourth, as a natural extension for addressing model multiplicity – the coexistence of multiple models with comparable predictive performance – it yields a distributional CE that incorporates model uncertainty by mixing the posterior distributions of the individual models using Bayesian model weights (see Section 2.4). The remainder of this paper is organized as follows. Section 2 describes the proposed method, Section 3 presents the experiments, Section 4 details related work, and Section 5 provides a discussion. 2 Counterfactual Explanations via Generalized Bayes 2.1 Counterfactual Explanation Formulation Let D=(xi,yi)i=1nD=\(x_i,y_i)\_i=1^n, where xi∈⊆ℝmx_i ^m is an input and yi∈y_i is its target value. We take =ℝY=R for regression and =0,1Y=\0,1\ for binary classification. Let f:→f:X denote a prediction model trained on D. We write xbx_b for the input to be explained, y=f(xb)y=f(x_b) for its predicted outcome, and y∗∈y^* for the desired output. A counterfactual explanation for xbx_b is denoted by x~∈ x . The function d(xb,x~)d(x_b, x) represents the cost or distance associated with the change, and ℓ(f(x~),y∗) (f( x),y^*) measures failure to attain the desired output. Squared-error loss is commonly used for continuous outcomes and cross-entropy loss for binary outcomes. The CE problem is then formulated as follows [15]: x~CE∈argminx~∈[λℓ(f(x~),y∗)+d(xb,x~)] x^CE∈ x argmin [λ (f( x),y^*)+d(x_b,\ x) ] (1) The problem is therefore to find a x~ x that both brings the prediction close to the desired output y∗y^* (small loss) and remains close to the original input xbx_b (small distance). Here, λ≥0λ≥ 0 controls the weight assigned to the loss. A fixed λ is widely used for three reasons. First, it is often difficult to solve either the constrained problem of minimizing d(xb,x~)d(x_b, x) subject to ℓ(f(x~),y∗)≤ε (f( x),y^*)≤ or the bi-objective problem of jointly minimizing ℓ and d; many implementations therefore use a fixed-weight penalized objective. Second, an exact constraint may yield no feasible solution for nonlinear models or noisy predictions, whereas a fixed λ provides a soft constraint. Third, λ has a natural interpretation as the trade-off between goal attainment and feasibility, and thus expresses the degree of tolerance permitted in the explanation [3, 14, 5]. 2.2 Generalized Bayes and the Gibbs Posterior We next introduce the generalized Bayes framework and use it to recast the CE problem from a Bayesian perspective. Generalized Bayes updates a prior distribution through exponential weighting based on a task-specific loss function, without requiring the likelihood assumed in conventional Bayesian inference [1]. To review conventional Bayesian inference briefly, let D denote the observed data and θ∈Θθ∈ the unknown quantity of interest. Given a prior p(θ)p(θ) and a probability model p(y|x,θ)p(y|x,θ), the posterior is p(θ|D)∝p(D|θ)p(θ)p(θ|D) p(D|θ)p(θ) (2) Here, p(D|θ)p(D|θ) is the likelihood and, under conditional independence, can be written as ∏i=1np(yi|xi,θ) _i=1^np(y_i|x_i,θ). In generalized Bayes, a loss function L(D,θ)∈ℝ+L(D,θ) ^+ that evaluates θ against D replaces the likelihood, yielding the following generalized Bayes posterior (or Gibbs posterior): p(θ|D)∝exp(−ηL(D,θ))p(θ)p(θ|D) (-η L(D,θ))p(θ) (3) Here, η is called the learning rate and controls the concentration of the posterior distribution. Setting L(D,θ)=−logp(D|θ)L(D,θ)=- p(D|θ) and η=1η=1 recovers conventional Bayesian inference as a special case. Because the MAP estimate minimizes an objective defined by the loss and the prior, it can also be interpreted as a regularized M-estimator. Within PAC-Bayes theory, the Gibbs posterior arises as the distribution that optimally trades off expected loss against Kullback–Leibler divergence from the prior, and it can carry generalization guarantees even without a likelihood model [2]. These guarantees concern estimation of a parameter θ; we do not claim that an analogous guarantee applies to the CE posterior over x~ x introduced in Section 2.3. 2.3 Counterfactual Explanation as a Gibbs Posterior We now apply this framework to CEs and derive a Gibbs posterior over counterfactuals. Let x~ x be the object of inference and let p(x~|xb)p( x|x_b) be its prior distribution conditional on xbx_b. Using ℓ(f(x~),y∗) (f( x),y^*) as the loss yields p(x~|xb,y∗)∝exp(−ηℓ(f(x~),y∗))p(x~|xb)p( x|x_b,y^*) (-η (f( x),y^*))p( x|x_b) (4) We refer to this formulation as the Generalized Bayes CE (GBCE). Choosing the distance-based prior p(x~|xb)∝exp(−d(xb,x~))p( x|x_b) (-d(x_b, x)) gives p(x~|xb,y∗)∝exp(−ηℓ(f(x~),y∗)−d(xb,x~))p( x|x_b,y^*) \! (-η (f( x),y^*)-d(x_b, x) ) (5) We refer to this formulation as the Distance-Prior Generalized Bayes CE (DP-GBCE). Hereafter, unless otherwise noted, “Gibbs posterior” or “Gibbs” refers to the DP-GBCE. Algorithm 1 in Section 2.6 summarizes the overall computational procedure: it approximates the DP-GBCE (and the ModelUnc construction described in Section 2.4) via importance sampling, applies the decision rules, and computes the evaluation metrics. The distribution p(x~|xb)p( x|x_b) is a user-specified prior over counterfactual candidates and does not represent a data-generating process. Depending on the distance function, the Gibbs-type (or Boltzmann-type) prior exp(−d(xb,x~)) (-d(x_b, x)) includes the normal distribution (for squared Euclidean distance) and the Laplace distribution (for L1L_1 distance) as special cases. Proposition 1. Under the above setting, the MAP estimator of the generalized-Bayes posterior reduces to the following CE optimization problem: x^MAP=argmaxx~p(x~|xb,y∗)=argminx~[ηℓ(f(x~),y∗)+d(xb,x~)] x^MAP= xargmax\ p( x|x_b,y^*)= xargmin [η\, (f( x),y^*)+d(x_b, x) ] (6) Proof. p(x~|xb,y∗)∝exp(−ηℓ(f(x~),y∗)−d(xb,x~))p( x|x_b,y^*) (-η (f( x),y^*)-d(x_b, x)), and argmaxargmax is invariant under multiplication by a positive constant, so x^MAP=argmaxx~p(x~|xb,y∗)=argmaxx~exp(−ηℓ(f(x~),y∗)−d(xb,x~)) x^MAP=argmax_ x\ p( x|x_b,y^*)=argmax_ x\ (-η (f( x),y^*)-d(x_b, x)). Since the exponential function is monotonically increasing, this is equivalent to argminx~[ηℓ(f(x~),y∗)+d(xb,x~)]argmin_ x [η (f( x),y^*)+d(x_b, x) ]. ∎ This proposition shows that distance minimization in conventional CE is equivalent to MAP estimation of a Gibbs posterior within the generalized-Bayes framework. This equivalence is not merely a formal correspondence; it strengthens the theoretical justification of CE in the following sense. By reformulating CE as a Gibbs posterior, one can naturally introduce, beyond the single-point MAP estimate, a variety of decision rules that use the full posterior distribution (e.g., Bayes decision, CVaR-CE) as well as distribution-level evaluation metrics (e.g., success probability, stability). In other words, the justification put forward in this study is the structural claim that distance minimization is a special case of probabilistic inference, which in turn provides the theoretical grounds for applying the rich toolkit of Bayesian methods to CE. This expression admits the following interpretation. Remark 1 (Relation between λ and η). There are two possible views on the relationship between λ and η. (1) Under the view that identifies η=λη=λ, x~CE=x^MAP x^CE= x^MAP holds in a set-valued sense. (2) Under the view that treats λ and η as parameters with distinct roles and fixes λ=η=1λ=η=1, the role of the trade-off coefficient λ is absorbed into the scale of the prior (the variance σ2σ^2 of the Gaussian prior); setting d(xb,x~)=12σ2‖x~−xb‖2d(x_b, x)= 12σ^2\| x-x_b\|^2 likewise makes x~CE=x^MAP x^CE= x^MAP hold. For consistency of exposition and comparability with previous studies, we use η=1η=1 as the default. This corresponds to view (2), in which the trade-off is represented through the relative scales of the loss and the prior. The behavior for η≠1η≠ 1 is examined in the η-sensitivity analysis in Section 3 (Figure 3). Under view (1), varying η changes the spread of the posterior distribution. This interpretation also suggests extensions such as methods for selecting η within generalized Bayes. For the Gibbs posterior to be well defined as a probability distribution, however, the unnormalized density exp(−ηℓ−d) (-η -d) must be integrable. If the loss function is nonnegative and measurable, then exp(−ηℓ(f(x~),y∗))≤1 (-η (f( x),y^*))≤ 1, so the condition ∫exp(−d(xb,x~))x~<∞ _X (-d(x_b, x))\,d x<∞ on the prior side (which holds for the norm-based distances used in this study) is a sufficient condition. Note that the loss function need not correspond to a likelihood. The epistemological status of the posterior p(x~|xb,y∗)p( x|x_b,y^*) introduced in this section is discussed in Section 5. 2.4 Counterfactual Explanations Incorporating Model Uncertainty As a natural extension of the CE presented in Section 2.3, we now describe a formulation that accounts for model uncertainty. Here, model uncertainty broadly encompasses uncertainty about the configuration θ of the prediction model, including both the choice of model class (e.g., LightGBM, XGBoost, or Random Forest) and the parameter values within a model class. Our formulation accommodates both discrete θ (model-class choice) and continuous θ (model parameters). When multiple distinct models have comparable predictive performance, a phenomenon known as model multiplicity [4, 11, 7], a CE derived from one model may be invalid under another. Suppose that we place a prior on θ and use a loss L(D,θ)L(D,θ) in place of a likelihood to construct the generalized Bayes posterior p(θ|D)∝exp(−η′L(D,θ))p(θ)p(θ|D) (-η L(D,θ))p(θ). Integrating over this posterior gives p(x~|xb,y∗,D)∝∫exp(−ηℓ(fθ(x~),y∗))p(x~|xb)p(θ|D)θp( x|x_b,y^*,D) (-η\, (f_θ( x),y^*))\,p( x|x_b)\,p(θ|D)\,dθ (7) This is a fully generalized Bayesian CE that requires no likelihood; we refer to it as the Fully Generalized Bayes CE (FG-GBCE). Whereas the GBCE of Section 2.3 generalizes the prior p(x~|xb)p( x|x_b), the FG-GBCE extends the framework along a separate dimension by incorporating uncertainty in the prediction model θ. This construction, however, has a nested structure: it requires sampling from the posterior over θ and then over x~ x for each θ, making it computationally expensive. Moreover, tree-based ensemble models such as LightGBM and XGBoost do not have a fixed-dimensional continuous parameter vector, which makes a direct application of this formulation with continuous θ difficult. We therefore adopt, in this study, a practical special case in which θ is treated as a discrete variable over a finite set of model classes 1,…,K\1,…,K\. Taking a finite set of models f(1),…,f(K) \f^(1),…,f^(K) \, selected by some procedure, as candidates, and setting p(θ=k|D)=wkp(θ=k|D)=w_k, the integral in Eq. (7) reduces to a sum, yielding p(x~|xb,y∗,D)∝∑k=1Kwkexp(−ηℓ(f(k)(x~),y∗))p(x~|xb)p( x|x_b,y^*,D) _k=1^Kw_k (-η\, (f^(k)( x),y^*))\,p( x|x_b) (8) We refer to this formulation as ModelUnc (Model-Uncertainty DP-GBCE). It is the special case of FG-GBCE in which θ is discrete over a finite set of model classes. When accounting for model multiplicity, we set wk=1/Kw_k=1/K; when the models differ in performance, we set wk∝exp(−γRk)w_k (-γ R_k), where RkR_k is the cross-validation loss. Note that this differs from approaches such as stacking, in which multiple models are combined into a single model beforehand. Furthermore, if the losses L(D,θ)L(D,θ) and ℓ(fθ(x~),y∗) (f_θ( x),y^*) in the FG-GBCE (Eq. (7)) are each made to correspond to a likelihood, then, when θ represents the coefficients of a continuous parametric model, one obtains a fully Bayesian counterfactual explanation in the conventional sense (a Fully Bayesian CE) that jointly randomizes the prediction model’s parameters and the counterfactual. This special case, however, requires strong generative-model assumptions (the design of a likelihood). This point is also addressed in the discussion in Section 5. The proposed framework can therefore accommodate various forms of uncertainty and constraints, including model uncertainty. 2.5 Applications of the Posterior Distribution The resulting Gibbs posterior p(x~|xb,y∗)p( x|x_b,y^*) provides more than a single optimal solution: the full distribution can support a variety of decision-making and evaluation tasks. We consider three such uses: (1) selecting a CE from the posterior, (2) evaluating the resulting point estimate, and (3) evaluating the posterior distribution itself. 2.5.1 Selecting a CE from the Posterior Distribution Several methods can be used to select a point-valued CE from the posterior. The appropriate decision rule depends on the problem structure, including the shape of the success region and the cost of failure. We consider three main rules, representing distinct decision-making perspectives: correspondence with conventional CE (MAP), point summarization through minimization of posterior expected decision loss (Bayes decision), and risk aversion (CVaR-CE). (a) MAP decision (the mode) x~MAP=argmaxx~p(x~|xb,y∗) x^MAP=argmax_ x\ p( x|x_b,y^*). Appropriate when the success region is unimodal, or, even if multimodal, when the selected mode has a clear meaning. (b) Bayes decision (minimizing the expected decision loss) x~Bayes∈argminz∈x~∼p(x~|xb,y∗)[Δ(z,x~)] x^Bayes∈ z argmin\ E_ x p( x|x_b,y^*)[ (z, x)] (9) Here, Δ(z,x~) (z, x) is a decision loss on the input space that measures the discrepancy between the decision z and a counterfactual x~ x drawn from the posterior; note that it is not the goal-attainment loss ℓ itself. When Δ=∥⋅∥2 =\|·\|^2, the solution is the posterior mean; when Δ=∥⋅∥ =\|·\|, it is the geometric median. Note, however, that when the posterior distribution is multimodal or curved, the posterior mean can fall in a low-density region that does not attain the goal – a known limitation of point summarization. In such cases, an alternative is the medoid, which constrains the decision to lie within the posterior samples. (c) Risk-averse decision (CVaR-CE) Let rD_r denote the perturbation distribution representing the error incurred when executing the CE (in this study, the isotropic normal distribution (0,σδ2I)N(0, _δ^2I); see Section 3 for the specific settings), and let Lδ(z)=ℓ(f(z+δ),y∗)L_δ(z)= (f(z+δ),y^*) (with δ∼rδ _r) denote the loss under perturbation. We define CVaR-CE as the point that minimizes the conditional expectation (CVaR) of the upper (1−τ)(1-τ) tail of this perturbed loss: x~CVaR∈argminz∈CVaRτ(Lδ(z)),CVaRτ(Lδ(z))=δ∼r[Lδ(z)|Lδ(z)≥VaRτ(Lδ(z))] x^CVaR∈ z argmin\ CVaR_τ\! (L_δ(z) ), _τ\! (L_δ(z) )=E_δ _r\! [L_δ(z)\, |\,L_δ(z) _τ\! (L_δ(z) ) ] (10) Here, VaRτVaR_τ is the τ-quantile of Lδ(z)L_δ(z). In implementation, we restrict the candidate points z to samples from the posterior p(x~|xb,y∗)p( x|x_b,y^*) (see Section 3). This risk-averse rule limits exposure to poor outcomes under execution noise when success is uncertain or failure is costly. Beyond these, many other approaches are conceivable, such as using the medoid; extracting a set of CEs using the HPD (highest posterior density) region when diversity is to be ensured; chance-constrained decisions; and sample-based decisions that sample directly from the distribution. 2.5.2 Pointwise Evaluation Metrics for Individual CEs We define metrics for evaluating the quality of a single-point CE x~ x obtained by a decision rule. These are standard metrics based on existing studies [3, 14, 5], and we use four: achieved loss, change distance, robustness, and plausibility. These four metrics respectively correspond to the main requirements demanded of a CE: goal attainment, minimality of change, stability at execution time, and feasibility with respect to the data. Lpt(x~) L_pt( x) =ℓ(f(x~),y∗) = (f( x),y^*) (11) Dpt(x~) D_pt( x) =d(x~,xb) =d( x,x_b) (12) Rb(x~) Rb( x) =Pδ∼r[ℓ(f(x~+δ),y∗)≤εrb] =P_δ _r[ (f( x+δ),y^*)≤ _rb] (13) Plu(x~) Plu( x) =1q∑i∈Nq(x~)‖x~−xi‖ = 1q _i∈ N_q( x) \| x-x_i \| (14) LptL_pt represents the degree of attainment of the desired output y∗y^*; smaller is better. DptD_pt represents the amount of change from the original input; smaller is better. RbRb is the probability that the prediction remains within the success region when a small perturbation δ∼rδ _r (the same perturbation distribution as in Section 2.5.1(c)) is added to x~ x, and it reflects the stability of the CE (higher is better). PluPlu is the average distance to the q nearest neighbors in the training data, measuring whether the CE is a plausible point with respect to the data distribution (plausibility; smaller is better). 2.5.3 Evaluation Metrics for the Posterior Distribution Itself Independently of the decision rule, it is also important to evaluate the posterior distribution p(x~|xb,y∗)p( x|x_b,y^*) itself. The following metrics quantify the posterior probability of success, tail risk, dispersion, and variable-specific changes. SP SP =Px~∼p[ℓ(f(x~),y∗)≤εsp](success probability: higher is better) =P_ x p[ (f( x),y^*)≤ _sp] (success probability: higher is better) (15) Tail Tail =inft:P(ℓ(f(x~),y∗)≤t)≥1−α(tail of the distribution: smaller is better) = \t:P( (f( x),y^*)≤ t)≥ 1-α\ (tail of the distribution: smaller is better) (16) Stability Stability =tr[Cov(x~)](stability: smaller is better) =tr[Cov( x)] (stability: smaller is better) (17) VarImpj VarImp_j =[|x~j−xb,j|](distributional variable importance) =E[| x_j-x_b,j|] (distributional variable importance) (18) SPSP is the proportion of posterior samples satisfying the success condition ℓ≤εsp ≤ _sp, and it indicates how well the posterior distribution covers the success region. TailTail is the (1−α)(1-α)-quantile of the loss, i.e., VaR1−αVaR_1-α, and it captures the heaviness of the distribution’s tail (the worst-case loss level). StabilityStability is the trace of the covariance matrix of the posterior distribution, representing the magnitude of the CE’s variance – that is, the uniqueness and stability of the solution. VarImpjVarImp_j is the expected amount of change in each variable j, indicating at the distribution level which variables are important in the counterfactual. 2.6 Overview of the Computational Procedure Algorithm 1 summarizes the full procedure described above, from the construction of the DP-GBCE (and the ModelUnc of Section 2.4) through the application of the decision rules to the computation of the evaluation metrics. The experiments in Section 3 follow this procedure, performing an approximate computation via importance sampling. Algorithm 1 Computational procedure for DP-GBCE (and ModelUnc) 1:Base point xbx_b, target value y∗y^*, prediction model f (for ModelUnc, a set of models f(k)k=1K\f^(k)\_k=1^K with weights wk\w_k\), loss ℓ , temperature η, proposal-distribution scale σ (corresponding to the distance prior d(xb,x~)=12σ2‖x~−xb‖2d(x_b, x)= 12σ^2\| x-x_b\|^2), number of candidates NcandN_cand, number of samples NsampN_samp 2:Generate candidate points x~ii=1Ncand∼(xb,σ2I)\ x_i\_i=1^N_cand (x_b,σ^2I) from the proposal distribution ⊳ Match the proposal distribution to the distance prior ∝exp(−d(xb,x~)) (-d(x_b, x)) 3:for i=1,…,Ncandi=1,…,N_cand do 4: if single model (DP-GBCE) then 5: ui←exp(−ηℓ(f(x~i),y∗))u_i← \! (-η\, (f( x_i),y^*) ) 6: else (multiple models, ModelUnc) 7: ui←∑k=1Kwkexp(−ηℓ(f(k)(x~i),y∗))u_i← _k=1^Kw_k \! (-η\, (f^(k)( x_i),y^*) ) 8: end if 9:end for 10:Normalize the weights ui\u_i\ and draw NsampN_samp samples with replacement via importance resampling to obtain posterior samples x~ss=1Nsamp\ x_s\_s=1^N_samp 11:Apply a decision rule (MAP, Bayes decision, or CVaR-CE) to determine a single-point CE (Section 2.5.1) 12:Compute pointwise evaluation metrics for the selected CE (Section 2.5.2) 13:Compute distributional evaluation metrics for x~s\ x_s\ (Section 2.5.3) 14:return the CE selected by each decision rule and the corresponding evaluation metrics 3 Empirical Evaluation This section evaluates the proposed method on simulated and real data. Using two-dimensional simulated data, we first examine whether the DP-GBCE and model-uncertainty posteriors concentrate in the success region (Figures 1 and 2) and compare the trade-offs among the MAP, Bayes, and CVaR-CE decision rules (Table 2). We then examine whether similar trade-offs arise in ten dimensions (Table 3) and analyze the sensitivity of the decision metrics to the temperature parameter η (Figure 3). Finally, we evaluate whether these findings persist under practical conditions using real data. Together, these experiments illustrate the additional information provided by a distributional approach relative to point-estimate CE methods. 3.1 Simulated Data We consider a regression task. The two-dimensional data (d=2d=2, n=3,000n=3,000) were generated from the following nonlinear function: y=2.0sin(x1)+0.8x22−1.2x1x2+ε,ε∼(0, 0.32)y=2.0 (x_1)+0.8\,x_2^2-1.2\,x_1x_2+ , (0,\,0.3^2) The ten-dimensional data (d=10d=10, n=4,000n=4,000) were generated from y=1.5sin(x1)+0.8x22−x1x3+0.5x4+0.3∑j=57xj+ε,ε∼(0, 0.52)y=1.5 (x_1)+0.8\,x_2^2-x_1x_3+0.5\,x_4+0.3 _j=5^7x_j+ , (0,\,0.5^2) In both cases, the input variables were generated independently from (0,1)N(0,1). The target value y∗y^* was set to the 90th percentile of the training outcomes, and squared-error loss, ℓ(y^,y∗)=(y^−y∗)2 ( y,y^*)=( y-y^*)^2, was used. To train the prediction models used to compute the CEs, we applied FLAML AutoML with eight candidate estimators: LightGBM, XGBoost, depth-limited XGBoost, HistGradientBoosting, Random Forest, ExtraTrees, ElasticNet, and SGD. The evaluation metric was MSE, the search-time budget for each estimator was capped at 25 seconds, and the random seed was fixed at 42 throughout. Model selection followed FLAML’s default setting (eval_method="auto"), which uses five-fold cross-validation at the data scale considered here. The best-performing model was used for the Gibbs posterior, whereas the top K=3K=3 models by MSE were used for the model-uncertainty CE. The latter used equal weights, wk=1/Kw_k=1/K, in Eq. (8); we did not examine performance-based weights wk∝exp(−γRk)w_k (-γ R_k) in these experiments. We set η=1η=1. Table 1 shows the top three models selected by AutoML. Table 1: Simulated data: top three models selected by AutoML (ordered by MSE) Data Model MSE 2D ExtraTrees 0.146 XGBoost (depth-limited) 0.152 LightGBM 0.157 10D LightGBM 0.426 XGBoost (depth-limited) 0.436 XGBoost 0.437 We approximated the Gibbs posterior by importance sampling. Candidate points were generated from the distance-based Gaussian prior centered at xbx_b, weighted by exp(−ηℓ) (-η ), and resampled. We used Ncand=2×104N_cand=2× 10^4 and Nsamp=2×103N_samp=2× 10^3 (proposal scale σ=1.0σ=1.0) in 2D, and Ncand=3.5×104N_cand=3.5× 10^4 and Nsamp=3×103N_samp=3× 10^3 (σ=1.0σ=1.0) in 10D. The Bayes decision (Section 2.5.1(b)) was computed with decision loss Δ=∥⋅∥2 =\|·\|^2, taking the sample mean (posterior mean) of the posterior samples as the solution; we refer to this as the Mean decision in the tables and figures below. For CVaR-CE, we set τ=0.9τ=0.9, thereby minimizing the average loss in the worst 10% of outcomes. We used nperturb=64n_perturb=64 perturbations drawn from δ∼(0,σδ2I)δ (0, _δ^2I) (2D: σδ=0.2 _δ=0.2; 10D: σδ=0.15 _δ=0.15) and evaluated at most 800 candidate points in 2D and 1,200 in 10D. We approximated the MAP decision by selecting the posterior sample with the highest estimated density. Density was estimated by Gaussian kernel density estimation (KDE) in 2D and by a k-nearest-neighbor score in 10D, where KDE is less reliable because of the curse of dimensionality. Specifically, the score was rk−dr_k^-d, where rkr_k is the distance to the kkth neighbor and k=40k=40. The resulting MAP should therefore be regarded as an approximation, particularly in 10D. To assess the proposition empirically, we also included a direct-optimization baseline (DirectOpt) that minimizes λℓ(f(x~),y∗)+d(xb,x~)λ (f( x),y^*)+d(x_b, x), with d(xb,x~)=‖x~−xb‖22/(2σ2)d(x_b, x)=\| x-x_b\|_2^2/(2σ^2) matching the Gaussian distance prior used for candidate generation. Because tree-based models are not differentiable, we used the gradient-free Nelder–Mead method. To reduce the risk of convergence to a local optimum, we performed multistart optimization from 20 initial points: xbx_b and 19 points generated using the same proposal scale as in the importance-sampling procedure. We retained the solution with the smallest objective value, allowing at most 3,000 iterations per run. The proposition concerns the exact posterior mode, whereas the reported Gibbs MAP is selected from a finite Monte Carlo sample using an estimated density. Consequently, numerical agreement with DirectOpt is not guaranteed: differences may arise from finite candidate generation, resampling, density estimation, and local numerical optimization. DirectOpt is therefore reported as an objective-matched reference rather than as a validation test that must coincide with the approximate Gibbs MAP. For the evaluation metrics, we used nperturb=200n_perturb=200 perturbations for robustness RbRb and q=20q=20 neighbors for plausibility PluPlu. The success thresholds εrb _rb and εsp _sp were set to 0.25 in 2D and 0.35 in 10D, and the tail probability was set to α=0.1α=0.1. 3.1.1 Two-Dimensional Data (Figures 1, 2, and Table 1) Figure 1: 2D Gibbs posterior and CE decisions (MAP, Mean, CVaR-CE, DirectOpt⋆ ) Figure 1 shows samples from the two-dimensional Gibbs posterior (blue points) and the CE selected by each decision rule. The posterior is multimodal and curved, reflecting the geometry of the regions that attain the target value y∗y^* (the 90th percentile of the training outcomes). Mean (green) lies between the base point and the main low-loss cluster and consequently has a large goal-attainment loss. MAP (red), CVaR-CE (purple), and DirectOpt (black star) lie within or near the main cluster, although they select different locations and therefore exhibit different loss–distance–robustness trade-offs (Table 2). Figure 2: 2D model-uncertainty CE posterior and CE decisions (MAP, Mean, CVaR-CE, DirectOpt⋆ ) Figure 2 shows the posterior distribution of the model-uncertainty CE (a weighted mixture of the top three models) on the same data. Compared with the Gibbs posterior (Figure 1), mixing multiple models slightly increases the spread of the posterior distribution (StabilityStability: 2.511 vs. 2.464) and slightly lowers the success probability (SPSP: 0.434 vs. 0.443). The decision points remain in broadly similar regions. For CVaR-CE, the model-uncertainty mixture yields higher robustness in this run (RbRb: 0.555 vs. 0.465), illustrating that incorporating model uncertainty need not move every point-level metric in a uniformly adverse direction. Table 2: Two-dimensional simulated data: comparison of CE evaluation metrics (VIxjVI_x_j: distributional variable importance) Method Decision rule LptL_pt DptD_pt RbRb PluPlu SPSP TailTail StabilityStability VIx1VI_x_1 VIx2VI_x_2 Gibbs Mean 3.6534 0.5024 0.000 0.0941 0.443 1.7388 2.4636 0.808 1.167 MAP 0.4333 1.1313 0.455 0.1351 CVaR 0.1162 1.3255 0.465 0.1456 ModelUnc Mean 3.4477 0.5337 0.000 0.0978 0.434 1.9443 2.5109 0.849 1.145 MAP 0.4752 1.1205 0.420 0.1313 CVaR 0.0531 1.3980 0.555 0.1510 DirectOpt — 0.0116 1.1616 0.485 0.1305 SP, Tail, Stability, and VIxjVI_x_j are metrics of the posterior distribution as a whole and are common across decision rules within the same method; they are therefore shown only in the first row of each method group. Table 2 reveals a clear trade-off among the decision rules. Mean has the shortest change distance but a large goal-attainment loss (Gibbs: Lpt=3.653L_pt=3.653; ModelUnc: 3.448) and zero perturbation robustness. MAP lowers the loss substantially while retaining a moderate distance. CVaR-CE further improves goal attainment and robustness at a modest additional distance. DirectOpt attains the smallest loss (0.012) with a distance similar to MAP and robustness of 0.485. The difference between DirectOpt and the approximate Gibbs MAP should not be read as contradicting the proposition: DirectOpt minimizes the continuous objective directly, whereas the reported MAP depends on finite posterior sampling and density estimation. The VIx1VI_x_1 and VIx2VI_x_2 columns of Table 2 report the expected change VarImpjVarImp_j in each variable under each posterior. Under the Gibbs posterior, VIx1≈0.808VI_x_1≈ 0.808 and VIx2≈1.167VI_x_2≈ 1.167, indicating a larger change in x2x_2. In the data-generating process y=2.0sin(x1)+0.8x22−1.2x1x2y=2.0 (x_1)+0.8x_2^2-1.2x_1x_2, x2x_2 affects y through both a quadratic term and an interaction term, which may explain this larger posterior movement. 3.1.2 Ten-Dimensional Data (Table 2) Table 3: Ten-dimensional simulated data: comparison of CE evaluation metrics Method Decision rule LptL_pt DptD_pt RbRb PluPlu SPSP TailTail StabilityStability Gibbs Mean 1.0051 0.9314 0.190 2.2121 0.601 1.4130 9.6905 MAP 1.8156 0.9230 0.000 2.0186 CVaR 0.0439 3.2529 0.975 3.2542 ModelUnc Mean 0.9790 0.9015 0.295 2.2075 0.494 2.1290 9.6022 MAP 0.4047 1.7650 0.640 2.1765 CVaR 0.0005 2.8439 0.950 3.6927 DirectOpt — 0.0170 1.4687 0.540 2.2095 SP, Tail, and Stability are metrics of the posterior distribution as a whole and are common across decision rules within the same method; they are therefore shown only in the first row of each method group. Table 3 shows that the trade-off also depends on how the point decision is extracted in ten dimensions. Gibbs CVaR-CE attains low loss (Lpt=0.044L_pt=0.044) and high robustness (Rb=0.975Rb=0.975), while ModelUnc CVaR-CE attains near-zero loss (0.0005) and robustness of 0.950, both at comparatively large distances. DirectOpt also attains low loss (0.017) at a shorter distance than either CVaR decision. In contrast, the density-based Gibbs MAP approximation has high loss (1.816) and zero robustness in this run, whereas the ModelUnc MAP performs better (Lpt=0.405L_pt=0.405, Rb=0.640Rb=0.640). This instability is consistent with the difficulty of mode estimation from finite samples in high dimensions and underscores that the theoretical MAP equivalence does not imply equality between DirectOpt and a numerically estimated sample mode. The Gibbs success probability is higher in 10D than in 2D (0.601 vs. 0.443), but the candidate count and success threshold εsp _sp differ, so this difference should not be attributed solely to dimensionality. Table 4: Ten-dimensional simulated data: distributional variable importance VIxj=[|x~j−xb,j|]VI_x_j=E[| x_j-x_b,j|] Method x1x_1 x2x_2 x3x_3 x4x_4 x5x_5 x6x_6 x7x_7 x8x_8 x9x_9 x10x_10 Gibbs 0.843 1.044 0.763 0.836 0.807 0.779 0.818 0.770 0.789 0.805 ModelUnc 0.849 0.984 0.789 0.807 0.802 0.793 0.796 0.797 0.792 0.799 Table 4 shows the distributional variable importance VarImpj=[|x~j−xb,j|]VarImp_j=E[| x_j-x_b,j|] for the ten-dimensional data. Under the Gibbs posterior, the noise variables x8x_8–x10x_10 provide a prior-driven baseline of approximately 0.77–0.81. Most other variables are close to this range, while x2x_2 has the largest expected change (1.044). Because the data-generating process contains the quadratic term 0.8x220.8x_2^2, this upward deviation indicates that movement in x2x_2 is especially prominent in constructing counterfactuals for this base point. The remaining coordinates should be interpreted cautiously: values near the baseline can reflect the spread of candidate generation as well as weak or diffuse loss constraints, so VarImpjVarImp_j is most informative when read relative to the noise-variable baseline rather than as an absolute importance score. 3.1.3 Sensitivity Analysis for η (Figure 3) Figure 3: η-sensitivity analysis: mean ± SD of stability (left) and success probability (right) (Gibbs posterior; 5 randomly chosen base points xbx_b) Figure 3 shows the trend in the stability (left) and success probability (right) of the posterior distribution as η varies. The mean success probability increases monotonically from about 0.09 at η=0η=0 to 0.68 at η=2η=2. Mean stability drops sharply at first (from about 2.02 to 1.71 by η=1/3η=1/3) and then fluctuates around 1.6–1.7 rather than decreasing monotonically. Thus, larger η more consistently improves concentration on the success region, while the diversity metric exhibits a plateau with Monte Carlo and base-point variation. Figure 4: η-sensitivity analysis: mean CE loss (MAP, Mean, CVaR-CE) Figure 4 shows how the mean CE loss LptL_pt for each decision rule changes with η. MAP loss falls sharply from 8.24 at η=0η=0 to 0.045 at the default η=1η=1 and remains near zero thereafter. CVaR-CE maintains a low but non-monotone loss throughout the grid. Mean loss decreases from 7.62 to 1.05 as η rises to 2, but remains appreciably larger than the MAP and CVaR losses because the posterior mean can lie outside a curved or multimodal success region. These results support η=1η=1 as a reasonable working value in this experiment while also showing that the decision rules respond differently to temperature. 3.2 Real Data (Google Trends, One Piece–Related) Unlike the simulation, in which the true function is known, the real-data experiments focus on verifying whether the proposed method is practically usable in real-world settings where the true function is unknown. We therefore leave the visualization of the posterior distribution and the η-sensitivity analysis to the results on simulated data, and focus here on comparing the evaluation metrics across decision rules. The real-data analysis uses monthly Japanese search-trend data obtained from Google Trends (January 2016 to December 2025; 119 observations). The target variable is search interest in “One Piece” (Wanpiisu). We retrieved the manga/anime topic through Google Trends’ topic-selection feature, thereby excluding searches for the unrelated clothing item. The explanatory variables are one-month-lagged (lag1) search trends for ten related terms: Golden Week, summer vacation, winter vacation, New Year, Jump (the magazine), LINE Manga, free manga, recommended manga, anime, and Netflix (Table 6). In the 10D analysis, these variables correspond, in order, to x1,…,x10x_1,…,x_10. The 2D analysis uses only x1′x_1 , summer vacation_lag1, and x2′x_2 , Jump_lag1. As in the simulation, we trained the prediction models using FLAML AutoML, evaluated them by MSE, and retained the top K=3K=3 models. Table 5 lists the selected models. Table 5: Real data: top three models selected by AutoML (ordered by MSE) Data Model MSE 2D ExtraTrees 58.08 LightGBM 60.34 Random Forest 62.80 10D Random Forest 51.93 HistGradientBoosting 53.35 XGBoost 54.50 Note: MSE is the squared error of the Google Trends index (0–100 scale). The target y∗y^* was set to the 90th percentile of the observed outcomes, and squared-error loss was used. We generated Ncand=2×104N_cand=2× 10^4 candidates in 2D and Ncand=3.5×104N_cand=3.5× 10^4 in 10D, with Nsamp=2×103N_samp=2× 10^3 posterior samples in both cases. Although the 10D simulation used Nsamp=3×103N_samp=3× 10^3, the real-data analysis repeats the procedure with every eligible observation as xbx_b; we therefore used a common sample size to facilitate comparison between 2D and 10D. To reflect the 0–100 Google Trends scale, we set the proposal scale to σ=10.0σ=10.0 in 2D and σ=8.0σ=8.0 in 10D. The evaluation settings were adjusted accordingly: εsp=εrb=15.0 _sp= _rb=15.0 (equivalent under squared-error loss to an absolute prediction error of approximately 3.9 points), α=0.1α=0.1, and r=(0,5.02I)D_r=N(0,5.0^2I) in both dimensions. For CVaR-CE, we used τ=0.9τ=0.9, nperturb=32n_perturb=32, and at most 400 candidate points. For robustness RbRb, we used nperturb=100n_perturb=100; for plausibility PluPlu, we used q=10q=10 neighbors. Candidate points were not clipped to the feasible range [0,100][0,100]. The base points xbx_b comprised all observations whose predicted values were below y∗y^* (2D: nxb=110n_xb=110; 10D: nxb=106n_xb=106), and we report the mean ± SD of each metric. MAP density estimation followed the simulation procedure: KDE in 2D and k-nearest-neighbor density estimation with k=40k=40 in 10D. Table 6: Descriptive statistics of the real data (Google Trends) Variable Mean SD Min Median Max One Piece 30.46 10.27 18.00 28.00 100.00 Golden Week_lag1 6.94 14.28 0.00 1.00 100.00 Summer vacation_lag1 24.50 29.85 4.00 9.00 100.00 Winter vacation_lag1 14.39 20.13 2.00 6.00 100.00 New Year_lag1 8.75 19.40 0.00 1.00 98.00 Jump_lag1 67.29 14.98 38.00 66.00 100.00 LINE Manga_lag1 66.41 15.30 36.00 66.00 100.00 Free manga_lag1 66.05 16.97 40.00 66.00 100.00 Recommended manga_lag1 59.03 16.81 30.00 57.00 100.00 Anime_lag1 64.37 12.73 46.00 61.00 100.00 Netflix_lag1 51.76 22.85 14.00 59.00 89.00 Table 7: Real data (2D): mean ± SD across all samples (VIxjVI_x_j: distributional variable importance) Method Decision rule LptL_pt DptD_pt RbRb PluPlu SPSP TailTail StabilityStability VIx1′VI_x_1 VIx2′VI_x_2 Gibbs Mean 29.395± 115.667 19.134± 6.809 0.510± 0.254 10.080± 4.490 1.000± 0.000 2.054± 0.637 129.296± 91.664 11.222± 4.630 16.498± 5.649 MAP 1.040± 0.684 19.320± 8.345 0.494± 0.200 9.286± 5.346 CVaR 1.532± 0.953 29.292± 9.398 0.759± 0.257 13.514± 6.580 ModelUnc Mean 47.255± 151.110 19.776± 7.071 0.502± 0.276 10.814± 3.539 0.826± 0.167 52.008± 52.741 139.366± 50.700 13.446± 4.149 15.656± 5.338 MAP 7.115± 24.843 19.789± 9.048 0.505± 0.227 10.102± 5.008 CVaR 10.786± 28.001 29.401± 9.099 0.731± 0.277 14.239± 7.146 SP, Tail, Stability, and VIxjVI_x_j are metrics of the posterior distribution as a whole and are common across decision rules within the same method; they are therefore shown only in the first row of each method group. Table 7 reports the 2D Google Trends results as the mean ± SD across 110 base points. For the Gibbs posterior, SP=1.000±0.000SP=1.000± 0.000, meaning that all posterior samples satisfied the success criterion for every base point. The change distance for Gibbs CVaR-CE, Dpt=29.3±9.4D_pt=29.3± 9.4, is much larger than in the simulation (≈3.2≈ 3.2), reflecting the 0–100 scale of the Google Trends variables. Gibbs MAP attains a low loss of Lpt=1.040±0.684L_pt=1.040± 0.684 while limiting the amount of change. The Mean decision, by contrast, has Lpt=29.4±115.7L_pt=29.4± 115.7, indicating unstable goal attainment. Relative to Gibbs, ModelUnc has a lower SPSP (0.826 vs. 1.000) and a substantially higher TailTail (52.0 vs. 2.1), suggesting that mixing multiple models produces a heavier loss tail. The VIxjVI_x_j columns report the expected changes in summer vacation_lag1 (x1′x_1 ) and Jump_lag1 (x2′x_2 ). Under the Gibbs posterior, VIx1′≈11.2VI_x_1 ≈ 11.2 and VIx2′≈16.5VI_x_2 ≈ 16.5, averaged over nxb=110n_xb=110, indicating that changes in Jump-related search demand contribute more strongly to the CE for the One Piece search trend. Table 8: Real data (10D): mean ± SD across all samples Method Decision rule LptL_pt DptD_pt RbRb PluPlu SPSP TailTail StabilityStability Gibbs Mean 24.881± 34.005 13.547± 7.643 0.451± 0.194 33.400± 9.638 1.000± 0.000 1.505± 0.400 607.938± 72.132 MAP 0.345± 0.544 20.456± 10.013 0.464± 0.196 36.569± 9.783 CVaR 0.914± 1.638 33.701± 6.918 0.761± 0.196 44.213± 10.041 ModelUnc Mean 20.446± 24.708 13.826± 7.426 0.449± 0.194 33.362± 9.635 0.841± 0.118 21.433± 11.446 617.107± 85.673 MAP 5.518± 13.223 20.449± 9.899 0.507± 0.178 36.558± 9.203 CVaR 3.909± 9.061 34.689± 6.427 0.767± 0.197 44.868± 9.266 SP, Tail, and Stability are metrics of the posterior distribution as a whole and are common across decision rules within the same method; they are therefore shown only in the first row of each method group. Table 8 reports the 10D results. All Gibbs posterior samples satisfied the success criterion for every base point (SP=1.000±0.000SP=1.000± 0.000), and Gibbs MAP achieved a low loss of Lpt=0.345±0.544L_pt=0.345± 0.544. The Mean decision had the smallest distance (Dpt=13.5D_pt=13.5) but a comparatively large and variable loss (Lpt=24.9±34.0L_pt=24.9± 34.0), showing that the limitations of the posterior centroid in high dimensions also arise in the real-data analysis. Relative to Gibbs, ModelUnc had a lower SPSP (0.841 vs. 1.000) and a substantially higher TailTail (21.4 vs. 1.5). As in the 2D analysis, mixing multiple models therefore produced a heavier loss tail. Table 9: Real data (10D): distributional variable importance VIxjVI_x_j Method x1x_1 x2x_2 x3x_3 x4x_4 x5x_5 x6x_6 x7x_7 x8x_8 x9x_9 x10x_10 Gibbs 6.475 7.448 6.560 6.441 6.652 6.922 11.673 7.899 6.292 6.581 ModelUnc 6.443 7.162 6.482 6.313 6.814 7.109 11.723 8.126 6.562 6.597 For the ten-dimensional data, Table 9 shows VarImpjVarImp_j for each variable (the maximum is x7x_7 = free manga_lag1, at approximately 11.7, followed by x8x_8 = recommended manga_lag1, at approximately 7.9). A larger value means that the variable plays a more important role in constructing the counterfactual. These real-data experiments have several practical implications. First, the Mean decision exhibits substantial variation in LptL_pt (e.g., SD=115.7SD=115.7 in 2D), making it an unstable point estimate that should be used with caution. As noted in Section 2.5.1(b), this behavior reflects a known limitation of point summaries: under a multimodal or curved posterior, the mean can lie in a low-density region. It does not indicate a defect in the posterior distribution itself. MAP, by contrast, attains a low loss while limiting the amount of change, whereas CVaR-CE improves robustness at the cost of a larger change. The choice between MAP and CVaR-CE should therefore depend on the intervention cost and the consequences of failing to attain the target. Second, ModelUnc has a substantially higher TailTail than the single-model Gibbs posterior, suggesting that accounting for model-selection uncertainty can reveal greater tail risk than an analysis based on a single model. Third, the distributional variable importance VarImpjVarImp_j identifies Jump_lag1, free manga_lag1, and recommended manga_lag1 as strong contributors to the counterfactuals. Monitoring these search-demand indicators – interest in the serialization magazine and free or recommendation-oriented manga services – may therefore help inform the timing of measures intended to increase interest in “One Piece.” 4 Related Work Several studies have examined CEs within a Bayesian framework. As noted in Section 1, [12] formulate counterfactual generation as a probabilistic model of perturbations. Their hierarchical Bayesian framework incorporates validity and proximity as likelihood terms and samples diverse counterfactuals from the posterior. In particular, their population–subgroup–instance hierarchy improves robustness by steering counterfactuals toward high-density regions of the data and supports fairness assessments based on comparisons of recourse costs across protected subgroups. Although internally coherent as a probabilistic model, the method requires a computationally demanding likelihood-based inference procedure. Its use of HMC/NUTS also requires a differentiable classifier. By contrast, the DP-GBCE proposed here is based on generalized Bayes and constructs a pseudo-posterior that combines goal-attainment loss with change cost, without requiring an explicit generative model or a specified likelihood. It can therefore represent counterfactuals as a distribution for differentiable models such as neural networks, nondifferentiable learners such as random forests, and black-box predictors accessed through external APIs, provided that the loss can be evaluated. It also enables multiple decision rules to be compared within a common framework and is thus model agnostic. Another relevant study is [9], who define recourse for black-box classifiers in Bayesian terms and derive a framework that identifies points with a high probability of receiving the desired prediction by minimizing posterior odds. To address the degradation of recourse validity under future changes to the classifier (model shift), they introduce a min–max optimization problem over ambiguity sets of class-conditional distributions defined by the Wasserstein distance, thereby obtaining recourse with a high success probability even in the worst case. Their approach, however, does not construct a distribution over counterfactuals from a prior and likelihood. Instead, it uses local sampling and density estimation to evaluate the odds of success and selects a single recourse point with a high success probability. The DP-GBCE proposed here instead places a probability distribution directly over counterfactual candidates. Through decision rules such as MAP, the posterior mean, and CVaR, it provides a unified way to assess and summarize uncertainty arising from multimodality or dispersion in the success region. Thus, whereas [9] focus on robustness to model changes, our focus is robust decision-making under uncertainty about the success region. Research on distributions of CEs can be divided broadly into two categories: work that considers a distribution of counterfactuals for a single observation and work that evaluates distributions of counterfactuals across a population. The first category seeks to represent counterfactual diversity and uncertainty explicitly, thereby avoiding the instability associated with reliance on a single optimum. It includes methods that generate diverse sets of CEs and methods based on probabilistic generative models [[, e.g.,]]mothilal2020explaining, raman2023bayesian. The second category focuses on the fairness of recourse at the population level, including comparisons of the distribution of recourse costs (corresponding to DptD_pt in this study) across subgroups and analyses of institutional fairness [[, e.g.,]]ustun2019actionable. Our study belongs to the first category: it formulates the counterfactual for a single observation as a probability distribution and incorporates decision rules directly into this distributional CE framework. The robustness of a CE may be challenged by several factors. [4] classify robust CEs into four categories: (1) model changes (MC), or robustness to changes caused by retraining or distribution shift; (2) model multiplicity (M), or uncertainty arising from the coexistence of multiple models with comparable performance; (3) noisy execution (NE), or noise and errors introduced when a CE is implemented; and (4) input changes (IC), including the consistency of explanations for similar inputs. This study primarily addresses robustness under M. For this setting, [11] derive theoretical upper bounds on CE costs under predictive multiplicity based on disagreement among classifiers, and empirically evaluate the robustness of existing CE methods, including methods constrained to data-supported regions. The multi-objective optimization approach of [7], meanwhile, constructs a Pareto set by minimizing the losses for all models simultaneously and thereby makes the trade-off structure among solutions explicit. Our approach differs in that it integrates model uncertainty probabilistically through generalized Bayes and yields a distributional CE with embedded decision rules. The relationship between our framework and both the multi-objective approach of [7] and min–max approaches that provide worst-case guarantees requires clarification. The effective loss of the mixture posterior in Eq. (8), −η−1log∑kwke−ηℓk-η^-1 _kw_ke^-η _k, converges to the weighted expected loss ∑kwkℓk _kw_k _k as η→0η→ 0, corresponding to a linear scalarization of a multi-objective problem. As η→∞η→∞, however, it converges to minkℓk _k _k. The mixture posterior is therefore an OR-type aggregation that assigns mass to regions in which at least one model succeeds. Its direction of aggregation differs from the AND-type worst-case guarantee generally provided by min–max optimization, which seeks validity under every model. Consequently, our framework does not provide a worst-case guarantee. Instead, Bayesian marginalization exposes disagreement among models through greater posterior dispersion, reflected in increases in TailTail and StabilityStability, as observed in Section 3. When worst-case robustness is required, a possible extension would replace the mixture with a loss-side aggregation such as exp(−ηmaxkℓk) (-η _k _k). Establishing the properties of such an extension is left for future work. 5 Discussion We first showed that a counterfactual explanation obtained by distance minimization is equivalent to the MAP estimate of a Gibbs posterior over CEs. This result provides a theoretical justification, from the perspective of Bayesian inference, for the conventional cost-minimization formulation. We then extended the framework by proposing a CE posterior that integrates model uncertainty probabilistically. For the resulting posterior, we introduced two decision rules beyond MAP: a Bayes decision that minimizes expected decision loss and CVaR-CE. This provides a unified framework for multiple decision-making objectives. Finally, we proposed metrics for evaluating both individual CEs and the posterior distribution as a whole, and quantified the trade-offs among the decision rules through experiments on simulated data and real Google Trends data. The significance of these results can be summarized in four points. First, we theoretically established that CE distance minimization is equivalent to MAP estimation under a generalized Bayes Gibbs posterior. Conventional CEs have been justified on practical and philosophical grounds as minimum-cost changes that attain a target; by identifying this formulation as a special case of generalized Bayes, namely MAP estimation under a Gibbs measure, our study provides a broader theoretical foundation. Second, the proposed framework is model agnostic: it requires neither an explicit generative model nor a specified likelihood and can be applied to any learner, differentiable or otherwise, provided that its loss can be evaluated. It can therefore construct counterfactuals for black-box predictors that cannot be handled by methods requiring gradient information. Third, the posterior perspective places MAP, Bayes, and CVaR-CE decisions within a single framework, allowing the decision rule to be selected according to the structure of the problem, such as the shape of the success region or the cost of failure. Fourth, for model multiplicity, where several models have comparable predictive performance, the framework naturally yields a distributional CE that incorporates model uncertainty by mixing the model-specific posteriors using Bayesian weights. The interpretation of the posterior p(x~|xb,y∗)p( x|x_b,y^*) underlying these contributions warrants further discussion. The prior p(x~|xb)p( x|x_b) is a candidate distribution based solely on proximity to xbx_b, whereas the prediction model f supplies additional information by imposing the requirement that the target y∗y^* be attained. The posterior is constructed under this requirement. With a single prediction model (Section 2.3), the spread of the posterior reflects only the multiplicity of the success region, which would remain even if the true f were known perfectly. With multiple prediction models (Section 2.4), an additional source of uncertainty arises because the models impose different requirements; in principle, this uncertainty can be reduced by acquiring more information. As discussed in Section 2.4, ModelUnc is a practical approximation to FG-GBCE and has the classical evidence-based structure in which uncertainty about θ can diminish as the amount of data D increases. In the experiments in Section 3, ModelUnc produced larger values of TailTail and StabilityStability than the single-model Gibbs posterior. This reflects the wider set of candidates that remains when the models impose conflicting constraints, indicating that disagreement among the models is appropriately represented in the posterior. In this sense, the posterior is best understood not as an inference about an unknown fact, but as a framework for choosing among multiple valid solutions under a given constraint (the prediction model) using the tools of Bayesian decision theory, including MAP, expected-loss minimization, and CVaR. This study also has several limitations. First, the efficiency of the importance-sampling approximation decreases as the dimension of x increases. Second, we fixed the temperature parameter at η=1η=1 and did not establish a theoretically grounded method for selecting it. Third, our analysis focused primarily on continuous features and did not accommodate binary features or constraints distinguishing mutable from immutable features. Finally, although we fixed the random seed in both the simulated- and real-data experiments to ensure reproducibility, FLAML’s AutoML search depends on its time budget. The selected models and resulting estimates may therefore vary slightly with the budget. Several directions for future work follow from these limitations. More scalable sampling methods, such as Metropolis–Hastings (MH), could improve approximation efficiency in high dimensions. The temperature parameter η might be selected using theoretically grounded methods based on SafeBayes or PAC-Bayes. Extending the framework to binary features and to constraints on mutable and immutable features is also important. Beyond model multiplicity, the framework may be applicable to other forms of robustness, including noisy execution and input changes. A further direction is to use the posterior distribution over CEs directly in downstream decision-making, such as policy planning or personalized recommendations. Finally, an important extension is to make the FG-GBCE described in Section 2.4, and its special case of a fully Bayesian CE, computationally practical. These formulations can represent epistemic uncertainty about x~ x in a more classical sense, but they require a nested procedure that combines sampling from the posterior over the model parameters θ with sampling over x~ x. This procedure is computationally expensive, and tree-based models present the additional structural difficulty that they lack a fixed-dimensional continuous parameter vector θ. Future work should develop approximations that overcome these practical barriers and investigate extensions such as p(x~|D)p( x|D) that incorporate D to assess validity in terms of proximity to the observed data. References [1] P. G. Bissiri, C. C. Holmes, and S. G. Walker (2016) A general framework for updating belief distributions. Journal of the Royal Statistical Society Series B: Statistical Methodology 78 (5), p. 1103–1130. Cited by: §1, §2.2. [2] B. Guedj (2019) A primer on pac-bayesian learning. arXiv preprint arXiv:1901.05353. Cited by: §2.2. [3] R. Guidotti (2024) Counterfactual explanations and how to find them: literature review and benchmarking. Data Mining and Knowledge Discovery 38 (5), p. 2770–2824. Cited by: §1, §1, §2.1, §2.5.2. [4] J. Jiang, F. Leofante, A. Rago, and F. Toni (2024) Robust counterfactual explanations in machine learning: a survey. In Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence, p. 8086–8094. Cited by: §2.4, §4. [5] A. Karimi, G. Barthe, B. Schölkopf, and I. Valera (2022) A survey of algorithmic recourse: contrastive explanations and consequential recommendations. ACM Computing Surveys 55 (5), p. 1–29. Cited by: §1, §1, §2.1, §2.5.2. [6] A. Karimi, J. von Kügelgen, B. Schölkopf, and I. Valera (2020) Algorithmic recourse under imperfect causal knowledge: a probabilistic approach. In Advances in Neural Information Processing Systems, Vol. 33, p. 265–277. Cited by: §1. [7] K. Kinjo (2025) Robust counterfactual explanations under model multiplicity using multi-objective optimization. arXiv preprint arXiv:2501.05795. Cited by: §2.4, §4, §4. [8] D. Lewis (2013) Counterfactuals. John Wiley & Sons. Note: Originally published in 1973 Cited by: §1. [9] T. H. Nguyen, N. Bui, D. Nguyen, M. Yue, and V. A. Nguyen (2022) Robust bayesian recourse. In Uncertainty in Artificial Intelligence, p. 1498–1508. Cited by: §1, §4, §4. [10] M. Pawelczyk, C. Agarwal, S. Joshi, S. Upadhyay, and H. Lakkaraju (2022) Exploring counterfactual explanations through the lens of adversarial examples: a theoretical and empirical analysis. In International Conference on Artificial Intelligence and Statistics, p. 4574–4594. Cited by: §1. [11] M. Pawelczyk, K. Broelemann, and G. Kasneci (2020) On counterfactual explanations under predictive multiplicity. In Uncertainty in Artificial Intelligence, p. 809–818. Cited by: §2.4, §4. [12] N. Raman, D. Magazzeni, and S. Shah (2023) Bayesian hierarchical models for counterfactual estimation. In International Conference on Artificial Intelligence and Statistics, p. 1115–1128. Cited by: §1, §4. [13] B. Ustun, A. Spangher, and Y. Liu (2019) Actionable recourse in linear classification. In Proceedings of the 2019 Conference on Fairness, Accountability, and Transparency, p. 10–19. Cited by: §1. [14] S. Verma, V. Boonsanong, M. Hoang, K. Hines, J. Dickerson, and C. Shah (2024) Counterfactual explanations and algorithmic recourses for machine learning: a review. ACM Computing Surveys 56 (12), p. 1–42. Cited by: §1, §1, §2.1, §2.5.2. [15] S. Wachter, B. Mittelstadt, and C. Russell (2017) Counterfactual explanations without opening the black box: automated decisions and the GDPR. Harvard Journal of Law & Technology 31 (2), p. 841–887. Cited by: §1, §2.1. Appendix Variational Derivation of the Gibbs Posterior Rather than estimating the counterfactual x~ x as a single point, we treat it as a distribution q(x~)q( x). Let L(x~)=ℓ(f(x~),y∗)L( x)= (f( x),y^*) denote the loss measuring attainment of the target y∗y^*, and introduce a prior π(x~|xb)π( x\,|\,x_b) representing natural changes from the base point xbx_b. Following a variational principle, we seek to reduce the expected loss q[L(x~)]E_q[L( x)] while controlling divergence from the prior π(x~|xb)π( x\,|\,x_b) through KL(q∥π)KL(q\|π). The Donsker–Varadhan variational formula implies that, for any distribution q absolutely continuous with respect to π and any η>0η>0, x~∼q[L(x~)]+1ηKL(q∥π)≥−1ηlogx~∼π[e−ηL(x~)].E_ x q[L( x)]+ 1ηKL(q\|π)\,≥\,- 1η _ x π\! [e^-η L( x) ]. (19) The expectation on the right-hand side is the normalizing constant (partition function), namely the prior expectation of the exponentially weighted loss. The functional on the left-hand side, which expresses the trade-off between loss and proximity to the prior, is uniformly bounded below by the right-hand side. Equality holds if and only if q is the Gibbs posterior qη∝πe−ηLq_η π e^-η L. Thus, the Gibbs posterior arises as the distribution that optimally balances the loss against divergence from the prior. Specifically, qη∈argminqq[L(x~)]+1ηKL(q∥π),q_η∈ _q \E_q[L( x)]+ 1ηKL(q\|π) \, (20) with the closed-form expression qη(x~)∝π(x~|xb)exp−ηL(x~).q_η( x) π( x\,|\,x_b) \-η L( x)\. (21) The parameter η can therefore be interpreted as a temperature (or inverse-temperature) parameter controlling the trade-off between loss minimization and proximity to the prior, which represents the naturalness of a change. As η increases, the distribution concentrates in low-loss regions; as it decreases, the distribution remains more diffuse and closer to the prior.