A Generalized-Bayes Perspective on Counterfactual Explanations: Posterior-Based Decision-Making and Evaluation

# A Generalized-Bayes Perspective on Counterfactual Explanations: Posterior-Based Decision-Making and Evaluation
Source: [https://arxiv.org/html/2607.29077](https://arxiv.org/html/2607.29077)
###### 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\.

## 1Introduction

Counterfactual explanations have attracted considerable attention as a means of addressing the interpretability challenge in machine learning\[[3](https://arxiv.org/html/2607.29077#bib.bib3),[14](https://arxiv.org/html/2607.29077#bib.bib4),[5](https://arxiv.org/html/2607.29077#bib.bib15)\]\. 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](https://arxiv.org/html/2607.29077#bib.bib13),[6](https://arxiv.org/html/2607.29077#bib.bib14)\], and its connection to adversarial examples has likewise been noted\[[5](https://arxiv.org/html/2607.29077#bib.bib15),[14](https://arxiv.org/html/2607.29077#bib.bib4),[10](https://arxiv.org/html/2607.29077#bib.bib10)\]\.

CE methods typically minimize the cost of changing the original instance into the counterfactual, thereby encouraging similarity between the two\[[15](https://arxiv.org/html/2607.29077#bib.bib16),[3](https://arxiv.org/html/2607.29077#bib.bib3),[14](https://arxiv.org/html/2607.29077#bib.bib4),[5](https://arxiv.org/html/2607.29077#bib.bib15)\]\.

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](https://arxiv.org/html/2607.29077#bib.bib11)\]\. 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](https://arxiv.org/html/2607.29077#bib.bib1)\]\. 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](https://arxiv.org/html/2607.29077#bib.bib9)\], 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](https://arxiv.org/html/2607.29077#bib.bib7)\], 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\.

## 2Counterfactual Explanations via Generalized Bayes

### 2\.1Counterfactual Explanation Formulation

LetD=\{\(xi,yi\)\}i=1nD=\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{n\}, wherexi∈𝒳⊆ℝmx\_\{i\}\\in\\mathcal\{X\}\\subseteq\\mathbb\{R\}^\{m\}is an input andyi∈𝒴y\_\{i\}\\in\\mathcal\{Y\}is its target value\. We take𝒴=ℝ\\mathcal\{Y\}=\\mathbb\{R\}for regression and𝒴=\{0,1\}\\mathcal\{Y\}=\\\{0,1\\\}for binary classification\. Letf:𝒳→𝒴f:\\mathcal\{X\}\\rightarrow\\mathcal\{Y\}denote a prediction model trained onDD\. We writexbx\_\{b\}for the input to be explained,y=f​\(xb\)y=f\(x\_\{b\}\)for its predicted outcome, andy∗∈𝒴y^\{\*\}\\in\\mathcal\{Y\}for the desired output\. A counterfactual explanation forxbx\_\{b\}is denoted byx~∈𝒳\\tilde\{x\}\\in\\mathcal\{X\}\. The functiond​\(xb,x~\)d\(x\_\{b\},\\tilde\{x\}\)represents the cost or distance associated with the change, andℓ​\(f​\(x~\),y∗\)\\ell\(f\(\\tilde\{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](https://arxiv.org/html/2607.29077#bib.bib16)\]:

x~C​E∈argminx~∈𝒳​\[λ​ℓ​\(f​\(x~\),y∗\)\+d​\(xb,x~\)\]\\tilde\{x\}^\{CE\}\\in\\underset\{\\tilde\{x\}\\in\\mathcal\{X\}\}\{\\operatorname\{argmin\}\}\\left\[\\lambda\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\+d\(x\_\{b\},\\ \\tilde\{x\}\)\\right\]\(1\)The problem is therefore to find ax~\\tilde\{x\}that both brings the prediction close to the desired outputy∗y^\{\*\}\(small loss\) and remains close to the original inputxbx\_\{b\}\(small distance\)\. Here,λ≥0\\lambda\\geq 0controls the weight assigned to the loss\. A fixedλ\\lambdais widely used for three reasons\. First, it is often difficult to solve either the constrained problem of minimizingd​\(xb,x~\)d\(x\_\{b\},\\tilde\{x\}\)subject toℓ​\(f​\(x~\),y∗\)≤ε\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\\leq\\varepsilonor the bi\-objective problem of jointly minimizingℓ\\ellanddd; 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λ\\lambdaprovides a soft constraint\. Third,λ\\lambdahas a natural interpretation as the trade\-off between goal attainment and feasibility, and thus expresses the degree of tolerance permitted in the explanation\[[3](https://arxiv.org/html/2607.29077#bib.bib3),[14](https://arxiv.org/html/2607.29077#bib.bib4),[5](https://arxiv.org/html/2607.29077#bib.bib15)\]\.

### 2\.2Generalized 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](https://arxiv.org/html/2607.29077#bib.bib1)\]\.

To review conventional Bayesian inference briefly, letDDdenote the observed data andθ∈Θ\\theta\\in\\Thetathe unknown quantity of interest\. Given a priorp​\(θ\)p\(\\theta\)and a probability modelp​\(y\|x,θ\)p\(y\|x,\\theta\), the posterior is

p​\(θ\|D\)∝p​\(D\|θ\)​p​\(θ\)p\(\\theta\|D\)\\propto p\(D\|\\theta\)p\(\\theta\)\(2\)Here,p​\(D\|θ\)p\(D\|\\theta\)is the likelihood and, under conditional independence, can be written as∏i=1np​\(yi\|xi,θ\)\\prod\_\{i=1\}^\{n\}p\(y\_\{i\}\|x\_\{i\},\\theta\)\.

In generalized Bayes, a loss functionL​\(D,θ\)∈ℝ\+L\(D,\\theta\)\\in\\mathbb\{R\}^\{\+\}that evaluatesθ\\thetaagainstDDreplaces the likelihood, yielding the following generalized Bayes posterior \(or Gibbs posterior\):

p​\(θ\|D\)∝exp⁡\(−η​L​\(D,θ\)\)​p​\(θ\)p\(\\theta\|D\)\\propto\\exp\(\-\\eta L\(D,\\theta\)\)p\(\\theta\)\(3\)Here,η\\etais called the learning rate and controls the concentration of the posterior distribution\.

SettingL​\(D,θ\)=−log⁡p​\(D\|θ\)L\(D,\\theta\)=\-\\log p\(D\|\\theta\)andη=1\\eta=1recovers 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](https://arxiv.org/html/2607.29077#bib.bib2)\]\. These guarantees concern estimation of a parameterθ\\theta; we do not claim that an analogous guarantee applies to the CE posterior overx~\\tilde\{x\}introduced in Section 2\.3\.

### 2\.3Counterfactual Explanation as a Gibbs Posterior

We now apply this framework to CEs and derive a Gibbs posterior over counterfactuals\.

Letx~\\tilde\{x\}be the object of inference and letp​\(x~\|xb\)p\(\\tilde\{x\}\|x\_\{b\}\)be its prior distribution conditional onxbx\_\{b\}\. Usingℓ​\(f​\(x~\),y∗\)\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)as the loss yields

p​\(x~\|xb,y∗\)∝exp⁡\(−η​ℓ​\(f​\(x~\),y∗\)\)​p​\(x~\|xb\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)\\propto\\exp\(\-\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\)p\(\\tilde\{x\}\|x\_\{b\}\)\(4\)We refer to this formulation as the Generalized Bayes CE \(GBCE\)\.

Choosing the distance\-based priorp​\(x~\|xb\)∝exp⁡\(−d​\(xb,x~\)\)p\(\\tilde\{x\}\|x\_\{b\}\)\\propto\\exp\(\-d\(x\_\{b\},\\tilde\{x\}\)\)gives

p​\(x~\|xb,y∗\)∝exp⁡\(−η​ℓ​\(f​\(x~\),y∗\)−d​\(xb,x~\)\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)\\propto\\exp\\\!\\left\(\-\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\-d\(x\_\{b\},\\tilde\{x\}\)\\right\)\(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](https://arxiv.org/html/2607.29077#alg1)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 distributionp​\(x~\|xb\)p\(\\tilde\{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\) priorexp⁡\(−d​\(xb,x~\)\)\\exp\(\-d\(x\_\{b\},\\tilde\{x\}\)\)includes the normal distribution \(for squared Euclidean distance\) and the Laplace distribution \(forL1L\_\{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^M​A​P=argmaxx~​p​\(x~\|xb,y∗\)=argminx~​\[η​ℓ​\(f​\(x~\),y∗\)\+d​\(xb,x~\)\]\\hat\{x\}^\{MAP\}=\\underset\{\\tilde\{x\}\}\{\\operatorname\{argmax\}\}\\ p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)=\\underset\{\\tilde\{x\}\}\{\\operatorname\{argmin\}\}\\left\[\\eta\\,\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\+d\(x\_\{b\},\\tilde\{x\}\)\\right\]\(6\)

###### Proof\.

p​\(x~\|xb,y∗\)∝exp⁡\(−η​ℓ​\(f​\(x~\),y∗\)−d​\(xb,x~\)\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)\\propto\\exp\(\-\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\-d\(x\_\{b\},\\tilde\{x\}\)\), andargmax\\operatorname\{argmax\}is invariant under multiplication by a positive constant, sox^M​A​P=argmaxx~⁡p​\(x~\|xb,y∗\)=argmaxx~⁡exp⁡\(−η​ℓ​\(f​\(x~\),y∗\)−d​\(xb,x~\)\)\\hat\{x\}^\{MAP\}=\\operatorname\{argmax\}\_\{\\tilde\{x\}\}\\ p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)=\\operatorname\{argmax\}\_\{\\tilde\{x\}\}\\ \\exp\(\-\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\-d\(x\_\{b\},\\tilde\{x\}\)\)\. Since the exponential function is monotonically increasing, this is equivalent toargminx~⁡\[η​ℓ​\(f​\(x~\),y∗\)\+d​\(xb,x~\)\]\\operatorname\{argmin\}\_\{\\tilde\{x\}\}\\left\[\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\+d\(x\_\{b\},\\tilde\{x\}\)\\right\]\. ∎

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\.

For consistency of exposition and comparability with previous studies, we useη=1\\eta=1as 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\\eta\\neq 1is examined in theη\\eta\-sensitivity analysis in Section 3 \(Figure 3\)\. Under view \(1\), varyingη\\etachanges the spread of the posterior distribution\. This interpretation also suggests extensions such as methods for selectingη\\etawithin generalized Bayes\.

For the Gibbs posterior to be well defined as a probability distribution, however, the unnormalized densityexp⁡\(−η​ℓ−d\)\\exp\(\-\\eta\\ell\-d\)must be integrable\. If the loss function is nonnegative and measurable, thenexp⁡\(−η​ℓ​\(f​\(x~\),y∗\)\)≤1\\exp\(\-\\eta\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\)\\leq 1, so the condition∫𝒳exp⁡\(−d​\(xb,x~\)\)​𝑑x~<∞\\int\_\{\\mathcal\{X\}\}\\exp\(\-d\(x\_\{b\},\\tilde\{x\}\)\)\\,d\\tilde\{x\}<\\inftyon 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 posteriorp​\(x~\|xb,y∗\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)introduced in this section is discussed in Section 5\.

### 2\.4Counterfactual 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θ\\thetaof 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θ\\theta\(model\-class choice\) and continuousθ\\theta\(model parameters\)\. When multiple distinct models have comparable predictive performance, a phenomenon known as model multiplicity\[[4](https://arxiv.org/html/2607.29077#bib.bib5),[11](https://arxiv.org/html/2607.29077#bib.bib8),[7](https://arxiv.org/html/2607.29077#bib.bib6)\], a CE derived from one model may be invalid under another\.

Suppose that we place a prior onθ\\thetaand use a lossL​\(D,θ\)L\(D,\\theta\)in place of a likelihood to construct the generalized Bayes posteriorp​\(θ\|D\)∝exp⁡\(−η′​L​\(D,θ\)\)​p​\(θ\)p\(\\theta\|D\)\\propto\\exp\(\-\\eta^\{\\prime\}L\(D,\\theta\)\)p\(\\theta\)\. Integrating over this posterior gives

p​\(x~\|xb,y∗,D\)∝∫exp⁡\(−η​ℓ​\(fθ​\(x~\),y∗\)\)​p​\(x~\|xb\)​p​\(θ\|D\)​𝑑θp\(\\tilde\{x\}\|x\_\{b\},y^\{\*\},D\)\\propto\\int\\exp\(\-\\eta\\,\\ell\(f\_\{\\theta\}\(\\tilde\{x\}\),y^\{\*\}\)\)\\,p\(\\tilde\{x\}\|x\_\{b\}\)\\,p\(\\theta\|D\)\\,d\\theta\(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 priorp​\(x~\|xb\)p\(\\tilde\{x\}\|x\_\{b\}\), the FG\-GBCE extends the framework along a separate dimension by incorporating uncertainty in the prediction modelθ\\theta\. This construction, however, has a nested structure: it requires sampling from the posterior overθ\\thetaand then overx~\\tilde\{x\}for eachθ\\theta, 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θ\\thetadifficult\.

We therefore adopt, in this study, a practical special case in whichθ\\thetais treated as a discrete variable over a finite set of model classes\{1,…,K\}\\\{1,\\ldots,K\\\}\. Taking a finite set of models\{f\(1\),…,f\(K\)\}\\left\\\{f^\{\(1\)\},\\ldots,f^\{\(K\)\}\\right\\\}, selected by some procedure, as candidates, and settingp​\(θ=k\|D\)=wkp\(\\theta=k\|D\)=w\_\{k\}, the integral in Eq\. \([7](https://arxiv.org/html/2607.29077#S2.E7)\) reduces to a sum, yielding

p​\(x~\|xb,y∗,D\)∝∑k=1Kwk​exp⁡\(−η​ℓ​\(f\(k\)​\(x~\),y∗\)\)​p​\(x~\|xb\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\},D\)\\propto\\sum\_\{k=1\}^\{K\}\{w\_\{k\}\\exp\(\-\\eta\\,\\ell\(f^\{\(k\)\}\(\\tilde\{x\}\),y^\{\*\}\)\)\\,p\(\\tilde\{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θ\\thetais discrete over a finite set of model classes\. When accounting for model multiplicity, we setwk=1/Kw\_\{k\}=1/K; when the models differ in performance, we setwk∝exp⁡\(−γ​Rk\)w\_\{k\}\\propto\\exp\(\-\\gamma R\_\{k\}\), whereRkR\_\{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 lossesL​\(D,θ\)L\(D,\\theta\)andℓ​\(fθ​\(x~\),y∗\)\\ell\(f\_\{\\theta\}\(\\tilde\{x\}\),y^\{\*\}\)in the FG\-GBCE \(Eq\. \([7](https://arxiv.org/html/2607.29077#S2.E7)\)\) are each made to correspond to a likelihood, then, whenθ\\thetarepresents 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\.5Applications of the Posterior Distribution

The resulting Gibbs posteriorp​\(x~\|xb,y∗\)p\(\\tilde\{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\.1Selecting 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~M​A​P=argmaxx~⁡p​\(x~\|xb,y∗\)\\tilde\{x\}^\{MAP\}=\\operatorname\{argmax\}\_\{\\tilde\{x\}\}\\ p\(\\tilde\{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~B​a​y​e​s∈argminz∈𝒳​𝔼x~∼p​\(x~\|xb,y∗\)​\[Δ​\(z,x~\)\]\\tilde\{x\}^\{Bayes\}\\in\\underset\{z\\in\\mathcal\{X\}\}\{\\operatorname\{argmin\}\}\\ \\mathbb\{E\}\_\{\\tilde\{x\}\\sim p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)\}\[\\Delta\(z,\\tilde\{x\}\)\]\(9\)Here,Δ​\(z,x~\)\\Delta\(z,\\tilde\{x\}\)is a decision loss on the input space that measures the discrepancy between the decisionzzand a counterfactualx~\\tilde\{x\}drawn from the posterior; note that it is not the goal\-attainment lossℓ\\ellitself\. WhenΔ=∥⋅∥2\\Delta=\\\|\\cdot\\\|^\{2\}, the solution is the posterior mean; whenΔ=∥⋅∥\\Delta=\\\|\\cdot\\\|, 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𝒟r\\mathcal\{D\}\_\{r\}denote the perturbation distribution representing the error incurred when executing the CE \(in this study, the isotropic normal distribution𝒩​\(0,σδ2​I\)\\mathcal\{N\}\(0,\\sigma\_\{\\delta\}^\{2\}I\); see Section 3 for the specific settings\), and letLδ​\(z\)=ℓ​\(f​\(z\+δ\),y∗\)L\_\{\\delta\}\(z\)=\\ell\(f\(z\+\\delta\),y^\{\*\}\)\(withδ∼𝒟r\\delta\\sim\\mathcal\{D\}\_\{r\}\) denote the loss under perturbation\. We define CVaR\-CE as the point that minimizes the conditional expectation \(CVaR\) of the upper\(1−τ\)\(1\-\\tau\)tail of this perturbed loss:

x~C​V​a​R∈argminz∈𝒳CVaRτ\(Lδ\(z\)\),CVaRτ\(Lδ\(z\)\)=𝔼δ∼𝒟r\[Lδ\(z\)\|Lδ\(z\)≥VaRτ\(Lδ\(z\)\)\]\\tilde\{x\}^\{CVaR\}\\in\\underset\{z\\in\\mathcal\{X\}\}\{\\operatorname\{argmin\}\}\\ \\mathrm\{CVaR\}\_\{\\tau\}\\\!\\left\(L\_\{\\delta\}\(z\)\\right\),\\qquad\\mathrm\{CVaR\}\_\{\\tau\}\\\!\\left\(L\_\{\\delta\}\(z\)\\right\)=\\mathbb\{E\}\_\{\\delta\\sim\\mathcal\{D\}\_\{r\}\}\\\!\\left\[L\_\{\\delta\}\(z\)\\,\\middle\|\\,L\_\{\\delta\}\(z\)\\geq\\mathrm\{VaR\}\_\{\\tau\}\\\!\\left\(L\_\{\\delta\}\(z\)\\right\)\\right\]\(10\)Here,VaRτ\\mathrm\{VaR\}\_\{\\tau\}is theτ\\tau\-quantile ofLδ​\(z\)L\_\{\\delta\}\(z\)\. In implementation, we restrict the candidate pointszzto samples from the posteriorp​\(x~\|xb,y∗\)p\(\\tilde\{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\.2Pointwise Evaluation Metrics for Individual CEs

We define metrics for evaluating the quality of a single\-point CEx~\\tilde\{x\}obtained by a decision rule\. These are standard metrics based on existing studies\[[3](https://arxiv.org/html/2607.29077#bib.bib3),[14](https://arxiv.org/html/2607.29077#bib.bib4),[5](https://arxiv.org/html/2607.29077#bib.bib15)\], 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\.

Lp​t​\(x~\)\\displaystyle L\_\{pt\}\(\\tilde\{x\}\)=ℓ​\(f​\(x~\),y∗\)\\displaystyle=\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\(11\)Dp​t​\(x~\)\\displaystyle D\_\{pt\}\(\\tilde\{x\}\)=d​\(x~,xb\)\\displaystyle=d\(\\tilde\{x\},x\_\{b\}\)\(12\)R​b​\(x~\)\\displaystyle Rb\(\\tilde\{x\}\)=Pδ∼𝒟r​\[ℓ​\(f​\(x~\+δ\),y∗\)≤εr​b\]\\displaystyle=P\_\{\\delta\\sim\\mathcal\{D\}\_\{r\}\}\[\\ell\(f\(\\tilde\{x\}\+\\delta\),y^\{\*\}\)\\leq\\varepsilon\_\{rb\}\]\(13\)P​l​u​\(x~\)\\displaystyle Plu\(\\tilde\{x\}\)=1q​∑i∈Nq​\(x~\)‖x~−xi‖\\displaystyle=\\frac\{1\}\{q\}\\sum\_\{i\\in N\_\{q\}\(\\tilde\{x\}\)\}\\left\\\|\\tilde\{x\}\-x\_\{i\}\\right\\\|\(14\)
Lp​tL\_\{pt\}represents the degree of attainment of the desired outputy∗y^\{\*\}; smaller is better\.Dp​tD\_\{pt\}represents the amount of change from the original input; smaller is better\.R​bRbis the probability that the prediction remains within the success region when a small perturbationδ∼𝒟r\\delta\\sim\\mathcal\{D\}\_\{r\}\(the same perturbation distribution as in Section 2\.5\.1\(c\)\) is added tox~\\tilde\{x\}, and it reflects the stability of the CE \(higher is better\)\.P​l​uPluis the average distance to theqqnearest 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\.3Evaluation Metrics for the Posterior Distribution Itself

Independently of the decision rule, it is also important to evaluate the posterior distributionp​\(x~\|xb,y∗\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)itself\. The following metrics quantify the posterior probability of success, tail risk, dispersion, and variable\-specific changes\.

S​P\\displaystyle SP=Px~∼p​\[ℓ​\(f​\(x~\),y∗\)≤εs​p\]\(success probability: higher is better\)\\displaystyle=P\_\{\\tilde\{x\}\\sim p\}\[\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\\leq\\varepsilon\_\{sp\}\]\\quad\\text\{\(success probability: higher is better\)\}\(15\)T​a​i​l\\displaystyle Tail=inf\{t:P​\(ℓ​\(f​\(x~\),y∗\)≤t\)≥1−α\}\(tail of the distribution: smaller is better\)\\displaystyle=\\inf\\\{t:P\(\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\\leq t\)\\geq 1\-\\alpha\\\}\\quad\\text\{\(tail of the distribution: smaller is better\)\}\(16\)S​t​a​b​i​l​i​t​y\\displaystyle Stability=tr​\[Cov​\(x~\)\]\(stability: smaller is better\)\\displaystyle=\\mathrm\{tr\}\[\\mathrm\{Cov\}\(\\tilde\{x\}\)\]\\quad\\text\{\(stability: smaller is better\)\}\(17\)V​a​r​I​m​pj\\displaystyle VarImp\_\{j\}=𝔼​\[\|x~j−xb,j\|\]\(distributional variable importance\)\\displaystyle=\\mathbb\{E\}\[\|\\tilde\{x\}\_\{j\}\-x\_\{b,j\}\|\]\\quad\\text\{\(distributional variable importance\)\}\(18\)
S​PSPis the proportion of posterior samples satisfying the success conditionℓ≤εs​p\\ell\\leq\\varepsilon\_\{sp\}, and it indicates how well the posterior distribution covers the success region\.T​a​i​lTailis the\(1−α\)\(1\-\\alpha\)\-quantile of the loss, i\.e\.,VaR1−α\\mathrm\{VaR\}\_\{1\-\\alpha\}, and it captures the heaviness of the distribution’s tail \(the worst\-case loss level\)\.S​t​a​b​i​l​i​t​yStabilityis 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\.V​a​r​I​m​pjVarImp\_\{j\}is the expected amount of change in each variablejj, indicating at the distribution level which variables are important in the counterfactual\.

### 2\.6Overview of the Computational Procedure

Algorithm[1](https://arxiv.org/html/2607.29077#alg1)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 1Computational procedure for DP\-GBCE \(and ModelUnc\)1:Base point

xbx\_\{b\}, target value

y∗y^\{\*\}, prediction model

ff\(for ModelUnc, a set of models

\{f\(k\)\}k=1K\\\{f^\{\(k\)\}\\\}\_\{k=1\}^\{K\}with weights

\{wk\}\\\{w\_\{k\}\\\}\), loss

ℓ\\ell, temperature

η\\eta, proposal\-distribution scale

σ\\sigma\(corresponding to the distance prior

d​\(xb,x~\)=12​σ2​‖x~−xb‖2d\(x\_\{b\},\\tilde\{x\}\)=\\frac\{1\}\{2\\sigma^\{2\}\}\\\|\\tilde\{x\}\-x\_\{b\}\\\|^\{2\}\), number of candidates

NcandN\_\{\\mathrm\{cand\}\}, number of samples

NsampN\_\{\\mathrm\{samp\}\}
2:Generate candidate points

\{x~i\}i=1Ncand∼𝒩​\(xb,σ2​I\)\\\{\\tilde\{x\}\_\{i\}\\\}\_\{i=1\}^\{N\_\{\\mathrm\{cand\}\}\}\\sim\\mathcal\{N\}\(x\_\{b\},\\sigma^\{2\}I\)from the proposal distribution⊳\\trianglerightMatch the proposal distribution to the distance prior∝exp⁡\(−d​\(xb,x~\)\)\\propto\\exp\(\-d\(x\_\{b\},\\tilde\{x\}\)\)

3:for

i=1,…,Ncandi=1,\\ldots,N\_\{\\mathrm\{cand\}\}do

4:ifsingle model \(DP\-GBCE\)then

5:

ui←exp⁡\(−η​ℓ​\(f​\(x~i\),y∗\)\)u\_\{i\}\\leftarrow\\exp\\\!\\left\(\-\\eta\\,\\ell\(f\(\\tilde\{x\}\_\{i\}\),y^\{\*\}\)\\right\)
6:else\(multiple models, ModelUnc\)

7:

ui←∑k=1Kwk​exp⁡\(−η​ℓ​\(f\(k\)​\(x~i\),y∗\)\)u\_\{i\}\\leftarrow\\textstyle\\sum\_\{k=1\}^\{K\}w\_\{k\}\\exp\\\!\\left\(\-\\eta\\,\\ell\(f^\{\(k\)\}\(\\tilde\{x\}\_\{i\}\),y^\{\*\}\)\\right\)
8:endif

9:endfor

10:Normalize the weights

\{ui\}\\\{u\_\{i\}\\\}and draw

NsampN\_\{\\mathrm\{samp\}\}samples with replacement via importance resampling to obtain posterior samples

\{x~s\}s=1Nsamp\\\{\\tilde\{x\}\_\{s\}\\\}\_\{s=1\}^\{N\_\{\\mathrm\{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\}\\\{\\tilde\{x\}\_\{s\}\\\}\(Section 2\.5\.3\)

14:returnthe CE selected by each decision rule and the corresponding evaluation metrics

## 3Empirical 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](https://arxiv.org/html/2607.29077#S3.F1)and[2](https://arxiv.org/html/2607.29077#S3.F2)\) and compare the trade\-offs among the MAP, Bayes, and CVaR\-CE decision rules \(Table[2](https://arxiv.org/html/2607.29077#S3.T2)\)\. We then examine whether similar trade\-offs arise in ten dimensions \(Table[3](https://arxiv.org/html/2607.29077#S3.T3)\) and analyze the sensitivity of the decision metrics to the temperature parameterη\\eta\(Figure[3](https://arxiv.org/html/2607.29077#S3.F3)\)\. 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\.1Simulated 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\.0​sin⁡\(x1\)\+0\.8​x22−1\.2​x1​x2\+ε,ε∼𝒩​\(0,0\.32\)y=2\.0\\sin\(x\_\{1\}\)\+0\.8\\,x\_\{2\}^\{2\}\-1\.2\\,x\_\{1\}x\_\{2\}\+\\varepsilon,\\quad\\varepsilon\\sim\\mathcal\{N\}\(0,\\,0\.3^\{2\}\)The ten\-dimensional data \(d=10d=10,n=4,000n=4\{,\}000\) were generated from

y=1\.5​sin⁡\(x1\)\+0\.8​x22−x1​x3\+0\.5​x4\+0\.3​∑j=57xj\+ε,ε∼𝒩​\(0,0\.52\)y=1\.5\\sin\(x\_\{1\}\)\+0\.8\\,x\_\{2\}^\{2\}\-x\_\{1\}x\_\{3\}\+0\.5\\,x\_\{4\}\+0\.3\\sum\_\{j=5\}^\{7\}x\_\{j\}\+\\varepsilon,\\quad\\varepsilon\\sim\\mathcal\{N\}\(0,\\,0\.5^\{2\}\)In both cases, the input variables were generated independently from𝒩​\(0,1\)\\mathcal\{N\}\(0,1\)\. The target valuey∗y^\{\*\}was set to the 90th percentile of the training outcomes, and squared\-error loss,ℓ​\(y^,y∗\)=\(y^−y∗\)2\\ell\(\\hat\{y\},y^\{\*\}\)=\(\\hat\{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 topK=3K=3models by MSE were used for the model\-uncertainty CE\. The latter used equal weights,wk=1/Kw\_\{k\}=1/K, in Eq\. \([8](https://arxiv.org/html/2607.29077#S2.E8)\); we did not examine performance\-based weightswk∝exp⁡\(−γ​Rk\)w\_\{k\}\\propto\\exp\(\-\\gamma R\_\{k\}\)in these experiments\. We setη=1\\eta=1\.

Table[1](https://arxiv.org/html/2607.29077#S3.T1)shows the top three models selected by AutoML\.

Table 1:Simulated data: top three models selected by AutoML \(ordered by MSE\)We approximated the Gibbs posterior by importance sampling\. Candidate points were generated from the distance\-based Gaussian prior centered atxbx\_\{b\}, weighted byexp⁡\(−η​ℓ\)\\exp\(\-\\eta\\ell\), and resampled\. We usedNcand=2×104N\_\{\\mathrm\{cand\}\}=2\\times 10^\{4\}andNsamp=2×103N\_\{\\mathrm\{samp\}\}=2\\times 10^\{3\}\(proposal scaleσ=1\.0\\sigma=1\.0\) in 2D, andNcand=3\.5×104N\_\{\\mathrm\{cand\}\}=3\.5\\times 10^\{4\}andNsamp=3×103N\_\{\\mathrm\{samp\}\}=3\\times 10^\{3\}\(σ=1\.0\\sigma=1\.0\) in 10D\. The Bayes decision \(Section 2\.5\.1\(b\)\) was computed with decision lossΔ=∥⋅∥2\\Delta=\\\|\\cdot\\\|^\{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\\tau=0\.9, thereby minimizing the average loss in the worst 10% of outcomes\. We usednperturb=64n\_\{\\mathrm\{perturb\}\}=64perturbations drawn fromδ∼𝒩​\(0,σδ2​I\)\\delta\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\delta\}^\{2\}I\)\(2D:σδ=0\.2\\sigma\_\{\\delta\}=0\.2; 10D:σδ=0\.15\\sigma\_\{\\delta\}=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 akk\-nearest\-neighbor score in 10D, where KDE is less reliable because of the curse of dimensionality\. Specifically, the score wasrk−dr\_\{k\}^\{\-d\}, whererkr\_\{k\}is the distance to thekkth neighbor andk=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~\)\\lambda\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)\+d\(x\_\{b\},\\tilde\{x\}\), withd​\(xb,x~\)=‖x~−xb‖22/\(2​σ2\)d\(x\_\{b\},\\tilde\{x\}\)=\\\|\\tilde\{x\}\-x\_\{b\}\\\|\_\{2\}^\{2\}/\(2\\sigma^\{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 usednperturb=200n\_\{\\mathrm\{perturb\}\}=200perturbations for robustnessR​bRbandq=20q=20neighbors for plausibilityP​l​uPlu\. The success thresholdsεr​b\\varepsilon\_\{rb\}andεs​p\\varepsilon\_\{sp\}were set to 0\.25 in 2D and 0\.35 in 10D, and the tail probability was set toα=0\.1\\alpha=0\.1\.

#### 3\.1\.1Two\-Dimensional Data \(Figures 1, 2, and Table 1\)

![Refer to caption](https://arxiv.org/html/2607.29077v1/fig/fig1_2d_gibbs_flaml.jpg)Figure 1:2D Gibbs posterior and CE decisions \(MAP, Mean, CVaR\-CE, DirectOpt⋆\\star\)Figure[1](https://arxiv.org/html/2607.29077#S3.F1)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 valuey∗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](https://arxiv.org/html/2607.29077#S3.T2)\)\.

![Refer to caption](https://arxiv.org/html/2607.29077v1/fig/fig2_2d_modelunc_flaml.jpg)Figure 2:2D model\-uncertainty CE posterior and CE decisions \(MAP, Mean, CVaR\-CE, DirectOpt⋆\\star\)Figure[2](https://arxiv.org/html/2607.29077#S3.F2)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](https://arxiv.org/html/2607.29077#S3.F1)\), mixing multiple models slightly increases the spread of the posterior distribution \(S​t​a​b​i​l​i​t​yStability: 2\.511 vs\. 2\.464\) and slightly lowers the success probability \(S​PSP: 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 \(R​bRb: 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 \(V​IxjVI\_\{x\_\{j\}\}: distributional variable importance\)SP, Tail, Stability, andV​IxjVI\_\{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](https://arxiv.org/html/2607.29077#S3.T2)reveals a clear trade\-off among the decision rules\. Mean has the shortest change distance but a large goal\-attainment loss \(Gibbs:Lp​t=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\.

TheV​Ix1VI\_\{x\_\{1\}\}andV​Ix2VI\_\{x\_\{2\}\}columns of Table[2](https://arxiv.org/html/2607.29077#S3.T2)report the expected changeV​a​r​I​m​pjVarImp\_\{j\}in each variable under each posterior\. Under the Gibbs posterior,V​Ix1≈0\.808VI\_\{x\_\{1\}\}\\approx 0\.808andV​Ix2≈1\.167VI\_\{x\_\{2\}\}\\approx 1\.167, indicating a larger change inx2x\_\{2\}\. In the data\-generating processy=2\.0​sin⁡\(x1\)\+0\.8​x22−1\.2​x1​x2y=2\.0\\sin\(x\_\{1\}\)\+0\.8x\_\{2\}^\{2\}\-1\.2x\_\{1\}x\_\{2\},x2x\_\{2\}affectsyythrough both a quadratic term and an interaction term, which may explain this larger posterior movement\.

#### 3\.1\.2Ten\-Dimensional Data \(Table 2\)

Table 3:Ten\-dimensional simulated data: comparison of CE evaluation metricsSP, 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](https://arxiv.org/html/2607.29077#S3.T3)shows that the trade\-off also depends on how the point decision is extracted in ten dimensions\. Gibbs CVaR\-CE attains low loss \(Lp​t=0\.044L\_\{pt\}=0\.044\) and high robustness \(R​b=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 \(Lp​t=0\.405L\_\{pt\}=0\.405,R​b=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εs​p\\varepsilon\_\{sp\}differ, so this difference should not be attributed solely to dimensionality\.

Table 4:Ten\-dimensional simulated data: distributional variable importanceV​Ixj=𝔼​\[\|x~j−xb,j\|\]VI\_\{x\_\{j\}\}=\\mathbb\{E\}\[\|\\tilde\{x\}\_\{j\}\-x\_\{b,j\}\|\]Table[4](https://arxiv.org/html/2607.29077#S3.T4)shows the distributional variable importanceV​a​r​I​m​pj=𝔼​\[\|x~j−xb,j\|\]VarImp\_\{j\}=\\mathbb\{E\}\[\|\\tilde\{x\}\_\{j\}\-x\_\{b,j\}\|\]for the ten\-dimensional data\. Under the Gibbs posterior, the noise variablesx8x\_\{8\}–x10x\_\{10\}provide a prior\-driven baseline of approximately 0\.77–0\.81\. Most other variables are close to this range, whilex2x\_\{2\}has the largest expected change \(1\.044\)\. Because the data\-generating process contains the quadratic term0\.8​x220\.8x\_\{2\}^\{2\}, this upward deviation indicates that movement inx2x\_\{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, soV​a​r​I​m​pjVarImp\_\{j\}is most informative when read relative to the noise\-variable baseline rather than as an absolute importance score\.

#### 3\.1\.3Sensitivity Analysis forη\\eta\(Figure 3\)

![Refer to caption](https://arxiv.org/html/2607.29077v1/fig/fig3_1_eta_stability_mean_sd.jpg)

![Refer to caption](https://arxiv.org/html/2607.29077v1/fig/fig3_2_eta_sp_mean_sd.jpg)

Figure 3:η\\eta\-sensitivity analysis: mean±\\pmSD of stability \(left\) and success probability \(right\) \(Gibbs posterior; 5 randomly chosen base pointsxbx\_\{b\}\)Figure[3](https://arxiv.org/html/2607.29077#S3.F3)shows the trend in the stability \(left\) and success probability \(right\) of the posterior distribution asη\\etavaries\. The mean success probability increases monotonically from about 0\.09 atη=0\\eta=0to 0\.68 atη=2\\eta=2\. Mean stability drops sharply at first \(from about 2\.02 to 1\.71 byη=1/3\\eta=1/3\) and then fluctuates around 1\.6–1\.7 rather than decreasing monotonically\. Thus, largerη\\etamore consistently improves concentration on the success region, while the diversity metric exhibits a plateau with Monte Carlo and base\-point variation\.

![Refer to caption](https://arxiv.org/html/2607.29077v1/fig/fig3_3_eta_ce_loss_mean.jpg)Figure 4:η\\eta\-sensitivity analysis: mean CE loss \(MAP, Mean, CVaR\-CE\)Figure[4](https://arxiv.org/html/2607.29077#S3.F4)shows how the mean CE lossLp​tL\_\{pt\}for each decision rule changes withη\\eta\. MAP loss falls sharply from 8\.24 atη=0\\eta=0to 0\.045 at the defaultη=1\\eta=1and 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η\\etarises 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\\eta=1as a reasonable working value in this experiment while also showing that the decision rules respond differently to temperature\.

### 3\.2Real 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η\\eta\-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](https://arxiv.org/html/2607.29077#S3.T6)\)\. In the 10D analysis, these variables correspond, in order, tox1,…,x10x\_\{1\},\\ldots,x\_\{10\}\. The 2D analysis uses onlyx1′x\_\{1\}^\{\\prime\}, summer vacation\_lag1, andx2′x\_\{2\}^\{\\prime\},Jump\_lag1\.

As in the simulation, we trained the prediction models using FLAML AutoML, evaluated them by MSE, and retained the topK=3K=3models\. Table[5](https://arxiv.org/html/2607.29077#S3.T5)lists the selected models\.

Table 5:Real data: top three models selected by AutoML \(ordered by MSE\)Note: MSE is the squared error of the Google Trends index \(0–100 scale\)\.

The targety∗y^\{\*\}was set to the 90th percentile of the observed outcomes, and squared\-error loss was used\. We generatedNcand=2×104N\_\{\\mathrm\{cand\}\}=2\\times 10^\{4\}candidates in 2D andNcand=3\.5×104N\_\{\\mathrm\{cand\}\}=3\.5\\times 10^\{4\}in 10D, withNsamp=2×103N\_\{\\mathrm\{samp\}\}=2\\times 10^\{3\}posterior samples in both cases\. Although the 10D simulation usedNsamp=3×103N\_\{\\mathrm\{samp\}\}=3\\times 10^\{3\}, the real\-data analysis repeats the procedure with every eligible observation asxbx\_\{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\\sigma=10\.0in 2D andσ=8\.0\\sigma=8\.0in 10D\. The evaluation settings were adjusted accordingly:εs​p=εr​b=15\.0\\varepsilon\_\{sp\}=\\varepsilon\_\{rb\}=15\.0\(equivalent under squared\-error loss to an absolute prediction error of approximately 3\.9 points\),α=0\.1\\alpha=0\.1, and𝒟r=𝒩​\(0,5\.02​I\)\\mathcal\{D\}\_\{r\}=\\mathcal\{N\}\(0,5\.0^\{2\}I\)in both dimensions\. For CVaR\-CE, we usedτ=0\.9\\tau=0\.9,nperturb=32n\_\{\\mathrm\{perturb\}\}=32, and at most 400 candidate points\. For robustnessR​bRb, we usednperturb=100n\_\{\\mathrm\{perturb\}\}=100; for plausibilityP​l​uPlu, we usedq=10q=10neighbors\. Candidate points were not clipped to the feasible range\[0,100\]\[0,100\]\. The base pointsxbx\_\{b\}comprised all observations whose predicted values were belowy∗y^\{\*\}\(2D:nx​b=110n\_\{xb\}=110; 10D:nx​b=106n\_\{xb\}=106\), and we report the mean±\\pmSD of each metric\. MAP density estimation followed the simulation procedure: KDE in 2D andkk\-nearest\-neighbor density estimation withk=40k=40in 10D\.

Table 6:Descriptive statistics of the real data \(Google Trends\)Table 7:Real data \(2D\): mean±\\pmSD across all samples \(V​IxjVI\_\{x\_\{j\}\}: distributional variable importance\)SP, Tail, Stability, andV​IxjVI\_\{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](https://arxiv.org/html/2607.29077#S3.T7)reports the 2D Google Trends results as the mean±\\pmSD across 110 base points\. For the Gibbs posterior,S​P=1\.000±0\.000SP=1\.000\\pm 0\.000, meaning that all posterior samples satisfied the success criterion for every base point\. The change distance for Gibbs CVaR\-CE,Dp​t=29\.3±9\.4D\_\{pt\}=29\.3\\pm 9\.4, is much larger than in the simulation \(≈3\.2\\approx 3\.2\), reflecting the 0–100 scale of the Google Trends variables\. Gibbs MAP attains a low loss ofLp​t=1\.040±0\.684L\_\{pt\}=1\.040\\pm 0\.684while limiting the amount of change\. The Mean decision, by contrast, hasLp​t=29\.4±115\.7L\_\{pt\}=29\.4\\pm 115\.7, indicating unstable goal attainment\. Relative to Gibbs, ModelUnc has a lowerS​PSP\(0\.826 vs\. 1\.000\) and a substantially higherT​a​i​lTail\(52\.0 vs\. 2\.1\), suggesting that mixing multiple models produces a heavier loss tail\. TheV​IxjVI\_\{x\_\{j\}\}columns report the expected changes in summer vacation\_lag1 \(x1′x\_\{1\}^\{\\prime\}\) andJump\_lag1 \(x2′x\_\{2\}^\{\\prime\}\)\. Under the Gibbs posterior,V​Ix1′≈11\.2VI\_\{x\_\{1\}^\{\\prime\}\}\\approx 11\.2andV​Ix2′≈16\.5VI\_\{x\_\{2\}^\{\\prime\}\}\\approx 16\.5, averaged overnx​b=110n\_\{xb\}=110, indicating that changes inJump\-related search demand contribute more strongly to the CE for the One Piece search trend\.

Table 8:Real data \(10D\): mean±\\pmSD across all samplesSP, 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](https://arxiv.org/html/2607.29077#S3.T8)reports the 10D results\. All Gibbs posterior samples satisfied the success criterion for every base point \(S​P=1\.000±0\.000SP=1\.000\\pm 0\.000\), and Gibbs MAP achieved a low loss ofLp​t=0\.345±0\.544L\_\{pt\}=0\.345\\pm 0\.544\. The Mean decision had the smallest distance \(Dp​t=13\.5D\_\{pt\}=13\.5\) but a comparatively large and variable loss \(Lp​t=24\.9±34\.0L\_\{pt\}=24\.9\\pm 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 lowerS​PSP\(0\.841 vs\. 1\.000\) and a substantially higherT​a​i​lTail\(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 importanceV​IxjVI\_\{x\_\{j\}\}For the ten\-dimensional data, Table[9](https://arxiv.org/html/2607.29077#S3.T9)showsV​a​r​I​m​pjVarImp\_\{j\}for each variable \(the maximum isx7x\_\{7\}= free manga\_lag1, at approximately 11\.7, followed byx8x\_\{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 inLp​tL\_\{pt\}\(e\.g\.,S​D=115\.7SD=115\.7in 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 higherT​a​i​lTailthan 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 importanceV​a​r​I​m​pjVarImp\_\{j\}identifiesJump\_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\.”

## 4Related Work

Several studies have examined CEs within a Bayesian framework\. As noted in Section 1,\[[12](https://arxiv.org/html/2607.29077#bib.bib9)\]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](https://arxiv.org/html/2607.29077#bib.bib7)\], 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](https://arxiv.org/html/2607.29077#bib.bib7)\]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 toDp​tD\_\{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](https://arxiv.org/html/2607.29077#bib.bib5)\]classify robust CEs into four categories: \(1\) model changes \(MC\), or robustness to changes caused by retraining or distribution shift; \(2\) model multiplicity \(MM\), 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 MM\. For this setting,\[[11](https://arxiv.org/html/2607.29077#bib.bib8)\]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](https://arxiv.org/html/2607.29077#bib.bib6)\], 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](https://arxiv.org/html/2607.29077#bib.bib6)\]and min–max approaches that provide worst\-case guarantees requires clarification\. The effective loss of the mixture posterior in Eq\. \([8](https://arxiv.org/html/2607.29077#S2.E8)\),−η−1​log​∑kwk​e−η​ℓk\-\\eta^\{\-1\}\\log\\sum\_\{k\}w\_\{k\}e^\{\-\\eta\\ell\_\{k\}\}, converges to the weighted expected loss∑kwk​ℓk\\sum\_\{k\}w\_\{k\}\\ell\_\{k\}asη→0\\eta\\to 0, corresponding to a linear scalarization of a multi\-objective problem\. Asη→∞\\eta\\to\\infty, however, it converges tomink⁡ℓk\\min\_\{k\}\\ell\_\{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 inT​a​i​lTailandS​t​a​b​i​l​i​t​yStability, as observed in Section 3\. When worst\-case robustness is required, a possible extension would replace the mixture with a loss\-side aggregation such asexp⁡\(−η​maxk⁡ℓk\)\\exp\(\-\\eta\\max\_\{k\}\\ell\_\{k\}\)\. Establishing the properties of such an extension is left for future work\.

## 5Discussion

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 posteriorp​\(x~\|xb,y∗\)p\(\\tilde\{x\}\|x\_\{b\},y^\{\*\}\)underlying these contributions warrants further discussion\. The priorp​\(x~\|xb\)p\(\\tilde\{x\}\|x\_\{b\}\)is a candidate distribution based solely on proximity toxbx\_\{b\}, whereas the prediction modelffsupplies additional information by imposing the requirement that the targety∗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 trueffwere 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θ\\thetacan diminish as the amount of dataDDincreases\. In the experiments in Section 3, ModelUnc produced larger values ofT​a​i​lTailandS​t​a​b​i​l​i​t​yStabilitythan 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 ofxxincreases\. Second, we fixed the temperature parameter atη=1\\eta=1and 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η\\etamight 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 aboutx~\\tilde\{x\}in a more classical sense, but they require a nested procedure that combines sampling from the posterior over the model parametersθ\\thetawith sampling overx~\\tilde\{x\}\. This procedure is computationally expensive, and tree\-based models present the additional structural difficulty that they lack a fixed\-dimensional continuous parameter vectorθ\\theta\. Future work should develop approximations that overcome these practical barriers and investigate extensions such asp​\(x~\|D\)p\(\\tilde\{x\}\|D\)that incorporateDDto 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 Methodology78\(5\),pp\. 1103–1130\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p4.1),[§2\.2](https://arxiv.org/html/2607.29077#S2.SS2.p2.1)\.
- \[2\]B\. Guedj\(2019\)A primer on pac\-bayesian learning\.arXiv preprint arXiv:1901\.05353\.Cited by:[§2\.2](https://arxiv.org/html/2607.29077#S2.SS2.p5.4)\.
- \[3\]R\. Guidotti\(2024\)Counterfactual explanations and how to find them: literature review and benchmarking\.Data Mining and Knowledge Discovery38\(5\),pp\. 2770–2824\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.1),[§1](https://arxiv.org/html/2607.29077#S1.p2.1),[§2\.1](https://arxiv.org/html/2607.29077#S2.SS1.p2.11),[§2\.5\.2](https://arxiv.org/html/2607.29077#S2.SS5.SSS2.p1.1)\.
- \[4\]J\. Jiang, F\. Leofante, A\. Rago, and F\. Toni\(2024\)Robust counterfactual explanations in machine learning: a survey\.InProceedings of the Thirty\-Third International Joint Conference on Artificial Intelligence,pp\. 8086–8094\.Cited by:[§2\.4](https://arxiv.org/html/2607.29077#S2.SS4.p1.3),[§4](https://arxiv.org/html/2607.29077#S4.p5.1)\.
- \[5\]A\. Karimi, G\. Barthe, B\. Schölkopf, and I\. Valera\(2022\)A survey of algorithmic recourse: contrastive explanations and consequential recommendations\.ACM Computing Surveys55\(5\),pp\. 1–29\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.1),[§1](https://arxiv.org/html/2607.29077#S1.p2.1),[§2\.1](https://arxiv.org/html/2607.29077#S2.SS1.p2.11),[§2\.5\.2](https://arxiv.org/html/2607.29077#S2.SS5.SSS2.p1.1)\.
- \[6\]A\. Karimi, J\. von Kügelgen, B\. Schölkopf, and I\. Valera\(2020\)Algorithmic recourse under imperfect causal knowledge: a probabilistic approach\.InAdvances in Neural Information Processing Systems,Vol\.33,pp\. 265–277\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.1)\.
- \[7\]K\. Kinjo\(2025\)Robust counterfactual explanations under model multiplicity using multi\-objective optimization\.arXiv preprint arXiv:2501\.05795\.Cited by:[§2\.4](https://arxiv.org/html/2607.29077#S2.SS4.p1.3),[§4](https://arxiv.org/html/2607.29077#S4.p5.1),[§4](https://arxiv.org/html/2607.29077#S4.p6.8)\.
- \[8\]D\. Lewis\(2013\)Counterfactuals\.John Wiley & Sons\.Note:Originally published in 1973Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p3.1)\.
- \[9\]T\. H\. Nguyen, N\. Bui, D\. Nguyen, M\. Yue, and V\. A\. Nguyen\(2022\)Robust bayesian recourse\.InUncertainty in Artificial Intelligence,pp\. 1498–1508\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p7.1),[§4](https://arxiv.org/html/2607.29077#S4.p2.1),[§4](https://arxiv.org/html/2607.29077#S4.p3.1)\.
- \[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\.InInternational Conference on Artificial Intelligence and Statistics,pp\. 4574–4594\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.1)\.
- \[11\]M\. Pawelczyk, K\. Broelemann, and G\. Kasneci\(2020\)On counterfactual explanations under predictive multiplicity\.InUncertainty in Artificial Intelligence,pp\. 809–818\.Cited by:[§2\.4](https://arxiv.org/html/2607.29077#S2.SS4.p1.3),[§4](https://arxiv.org/html/2607.29077#S4.p5.1)\.
- \[12\]N\. Raman, D\. Magazzeni, and S\. Shah\(2023\)Bayesian hierarchical models for counterfactual estimation\.InInternational Conference on Artificial Intelligence and Statistics,pp\. 1115–1128\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p7.1),[§4](https://arxiv.org/html/2607.29077#S4.p1.1)\.
- \[13\]B\. Ustun, A\. Spangher, and Y\. Liu\(2019\)Actionable recourse in linear classification\.InProceedings of the 2019 Conference on Fairness, Accountability, and Transparency,pp\. 10–19\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.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 Surveys56\(12\),pp\. 1–42\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p1.1),[§1](https://arxiv.org/html/2607.29077#S1.p2.1),[§2\.1](https://arxiv.org/html/2607.29077#S2.SS1.p2.11),[§2\.5\.2](https://arxiv.org/html/2607.29077#S2.SS5.SSS2.p1.1)\.
- \[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 & Technology31\(2\),pp\. 841–887\.Cited by:[§1](https://arxiv.org/html/2607.29077#S1.p2.1),[§2\.1](https://arxiv.org/html/2607.29077#S2.SS1.p2.12)\.

## Appendix

### Variational Derivation of the Gibbs Posterior

Rather than estimating the counterfactualx~\\tilde\{x\}as a single point, we treat it as a distributionq​\(x~\)q\(\\tilde\{x\}\)\. LetL​\(x~\)=ℓ​\(f​\(x~\),y∗\)L\(\\tilde\{x\}\)=\\ell\(f\(\\tilde\{x\}\),y^\{\*\}\)denote the loss measuring attainment of the targety∗y^\{\*\}, and introduce a priorπ​\(x~\|xb\)\\pi\(\\tilde\{x\}\\,\|\\,x\_\{b\}\)representing natural changes from the base pointxbx\_\{b\}\. Following a variational principle, we seek to reduce the expected loss𝔼q​\[L​\(x~\)\]\\mathbb\{E\}\_\{q\}\[L\(\\tilde\{x\}\)\]while controlling divergence from the priorπ​\(x~\|xb\)\\pi\(\\tilde\{x\}\\,\|\\,x\_\{b\}\)throughKL​\(q∥π\)\\mathrm\{KL\}\(q\\\|\\pi\)\. The Donsker–Varadhan variational formula implies that, for any distributionqqabsolutely continuous with respect toπ\\piand anyη\>0\\eta\>0,

𝔼x~∼q​\[L​\(x~\)\]\+1η​KL​\(q∥π\)≥−1η​log⁡𝔼x~∼π​\[e−η​L​\(x~\)\]\.\\mathbb\{E\}\_\{\\tilde\{x\}\\sim q\}\[L\(\\tilde\{x\}\)\]\+\\frac\{1\}\{\\eta\}\\mathrm\{KL\}\(q\\\|\\pi\)\\,\\geq\\,\-\\frac\{1\}\{\\eta\}\\log\\mathbb\{E\}\_\{\\tilde\{x\}\\sim\\pi\}\\\!\\left\[e^\{\-\\eta L\(\\tilde\{x\}\)\}\\right\]\.\(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 ifqqis the Gibbs posteriorqη∝π​e−η​Lq\_\{\\eta\}\\propto\\pi e^\{\-\\eta L\}\. Thus, the Gibbs posterior arises as the distribution that optimally balances the loss against divergence from the prior\. Specifically,

qη∈arg⁡minq⁡\{𝔼q​\[L​\(x~\)\]\+1η​KL​\(q∥π\)\},q\_\{\\eta\}\\in\\arg\\min\_\{q\}\\left\\\{\\mathbb\{E\}\_\{q\}\[L\(\\tilde\{x\}\)\]\+\\frac\{1\}\{\\eta\}\\mathrm\{KL\}\(q\\\|\\pi\)\\right\\\},\(20\)
with the closed\-form expression

qη​\(x~\)∝π​\(x~\|xb\)​exp⁡\{−η​L​\(x~\)\}\.q\_\{\\eta\}\(\\tilde\{x\}\)\\propto\\pi\(\\tilde\{x\}\\,\|\\,x\_\{b\}\)\\exp\\\{\-\\eta L\(\\tilde\{x\}\)\\\}\.\(21\)
The parameterη\\etacan 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η\\etaincreases, the distribution concentrates in low\-loss regions; as it decreases, the distribution remains more diffuse and closer to the prior\.