Computationally efficient variational-like approximations of possibilistic inferential models¹¹1This is an extended version of the paper (Cella and Martin, 2024) presented at the 8th International Conference on Belief Functions, in Belfast, UK.

For a long time, despite Bayesians’ foundational advantages, few statisticians were actually using Bayesian methods—the computational burden was simply too high. Things changed significantly when Monte Carlo methods brought Bayesian solutions within reach. Things changed again more recently with the advances in various approximate Bayesian computational methods, in particular, the variational approximations in Blei et al., (2017), Tran et al., (2021), and the references therein. The once clear lines between what was computationally feasible for Bayesians and for non-Bayesians have now been blurred, reinvigorating Bayesians’ efforts in modern applications. Dennis Lindley predicted that “[statisticians] will all be Bayesians in 2020” (Smith, 1995)—his prediction did not come true, but the Bayesian community is arguably stronger now than ever.

While Bayesian and frequentist are currently the two mainstream schools of thought in statistical inference, these are not the only perspectives. For example, Dempster–Shafer theory originated as an improvement to and generalization of both Bayesian inference and Fisher’s fiducial argument. Of particular interest to us here are the recent advances in inferential models (IMs, Martin and Liu, 2013; Martin and Liu, 2015b ; Martin, 2021a ; Martin, 2022b ), a framework that offers Bayesian-like, data-dependent, possibilistic quantification of uncertainty about unknowns but with built-in, frequentist-like reliability guarantees. IMs and other new/non-traditional frameworks are currently facing the same computational challenges that Bayesians faced years ago. That is, we know what we want to compute and why, but we are currently lacking the tools to do so efficiently. While traditional Monte Carlo methods are still useful (see Section 2), the imprecision that is key to the IM’s reliability guarantees implies that Monte Carlo methods alone are not enough. In contrast to Lindley’s prediction, for Efron’s speculation about fiducial and related methods—“Maybe Fisher’s biggest blunder will become a big hit in the 21st century!” (Efron, 1998)—to come true, imprecision-accommodating advances in Monte Carlo computations are imperative. The present paper’s contribution is in this direction.

We start with a relatively simple idea that leads to a general tool for computationally efficient and statistically reliable possibilistic inference. For reasons summarized in Section 2, our focus is on the possibility-theoretic brand of IMs, which are fully determined by the corresponding contour function or, equivalently, the contour function’s so-called $\alpha$-cuts. We leverage a well-known characterization of a possibility measure’s credal set as the collection of probability measures that assign probability at least $1-\alpha$ to the aforementioned $\alpha$-cuts of the contour function. In the IM context, the $\alpha$-cuts are $100(1-\alpha)$% confidence regions and, therefore, the contents of the IM’s credal set can be justifiably interpreted as “confidence distributions” (e.g., Xie and Singh, 2013; Schweder and Hjort, 2002). Specifically, the “most diffuse” element of the credal set, or “inner probabilistic approximation,” would assign probabilities as close to $1-\alpha$ as possible to each $\alpha$-cut. If we could approximate this distinguished member of the IM’s credal set, via Monte Carlo or otherwise, then we would be well on our way to carrying out most—if not all—of the relevant IM computations. The challenge is that, except for very special classes of problems (Martin, 2023a ), this “inner probabilistic approximation” would be rather complex. We can, however, obtain a relatively simple approximation if we only require an accurate approximation of, say, a single $\alpha$-cut. Towards this, we propose to introduce a simple parametric family of probability distributions, e.g., Gaussian, on the parameter space, with parameters partially depending on data, and then choosing those unspecified features of the parameter in such a way that the distribution (approximately) assigns probability $1-\alpha$ to a specified $\alpha$-cut, say, $\alpha=0.1$. This is akin to variational Bayes in the sense that we propose to approximate a complex probability distribution—in our case it is the “inner probabilistic approximation” of the IM’s possibility measure instead of a Bayesian posterior distribution—by an appropriately chosen member of a posited class of relatively simple probability distributions. The specific, technical details our proposal here are inspired by the recent work in Jiang et al., (2023) and the seemingly unrelated developments in Syring and Martin, (2019).

The remainder of the paper is organized as follows. After some background on possibilsitic IMs and their properties in Section 2, we present our basic but very general variational-like IM approximation in Section 3, which relies on a combination of Monte Carlo and stochastic approximation to tune the index of the posited parametric family of approximations. This is ideally suited for statistical inference problems involving low-dimensional parameters, and we offer a few such illustrations in Section 4. A more sophisticated version of the proposed variational-like approximation is presented in Section 5, which is more appropriate when the problem at hand is of higher dimension, but it is suited primarily for Gaussian approximation families—which is no practical restriction. With a more structured approximation, we can lower the computational cost by reducing the number of Monte Carlo evaluations of the IM contour. This is illustrated in several examples, including a relatively high-dimensional problem with lasso penalization and in problems involving nuisance parameters in parametric, nonparametric, and semiparametric models, respectively. The paper concludes in Section 7 with a concise summary and a discussion of some practically relevant extensions.

The first presentation of the IM framework (e.g., Martin and Liu, 2013; Martin and Liu, 2015b ) relied heavily on (nested) random sets and their corresponding belief functions. The more recent IM formulation in Martin, 2022b , building off of developments in Martin, (2015, 2018), defines the IM’s possibility contour using a probability-to-possibility transform applied to the relative likelihood. This small-but-important shift has principled motivations but we can only barely touch on this here. The present review focuses on the possibilistic IM formulation, its key properties, and the available computational strategies.

Consider a parametric statistical model $\{\mathsf{P}_{\theta}:\theta\in\mathbb{T}\}$ indexed by a parameter space $\mathbb{T}\subseteq\mathbb{R}^{d}$. Examples include $\mathsf{P}_{\theta}={\sf Bin}(n,\theta)$, $\mathsf{P}_{\theta}={\sf N}(\theta,1)$, $\mathsf{P}_{\theta}={\sf Gamma}(\alpha,\beta)$ with $\theta=(\alpha,\beta)$, and many others; see Section 4. Nonparametric problems—see Cella, (2024), Cella and Martin, (2022), and Martin, 2023b (, Sec. 5)—can also be handled, but we postpone discussion of this until Section 6. Suppose that the observable data $X^{n}=(X_{1},\ldots,X_{n})$ consists of iid samples from the distribution $\mathsf{P}_{\Theta}$, where $\Theta\in\mathbb{T}$ is the unknown/uncertain “true value.” The model and observed data $X^{n}=x^{n}$ together determine a likelihood function $\theta\mapsto L_{x^{n}}(\theta)$ and a corresponding relative likelihood function

Throughout we will implicitly assume that the denominator is finite, i.e., that the likelihood function is bounded for almost all $x^{n}$.

The relative likelihood, $\theta\mapsto R(x^{n},\theta)$, itself is a data-dependent possibility contour, i.e., it is a non-negative function such that $\sup_{\theta}R(x^{n},\theta)=1$ for almost all $x^{n}$. This contour, in turn, determines a possibility measure, or consonant plausibility function, that can be used for uncertainty quantification about $\Theta$, given $X^{n}=x^{n}$, which has been extensively studied (e.g., Shafer, 1982; Wasserman, 1990; Denœux, 2006, 2014). Specifically, uncertainty quantification here would be based on the possibility assignments

where $H$ represents a hypothesis about $\Theta$. This purely likelihood-driven possibility measure has a number of desirable properties: for example, inference based on thresholding the possibility scores (with model-independent thresholds) satisfies the likelihood principle (e.g., Birnbaum, 1962; Basu, 1975), and it is asymptotically consistent (under standard regularity conditions) in the sense that it converges in $\mathsf{P}_{\Theta}$-probability to a possibility measure concentrated on $\Theta$ as $n\to\infty$. What the purely likelihood-based possibility measure above lacks, however, is a calibration property (relative to the posited model) that gives meaning or belief-forming inferential weight to the “possibility” it assigns to hypotheses about the unknown $\Theta$. More specifically, if we offer a possibility measure as our quantification of statistical uncertainty, then its corresponding credal set should contain probability distributions that are statistically meaningful. We are assuming vacuous prior information, so there are no meaningful or distinguished Bayesian posterior distributions. The only other natural property to enforce on the credal set’s elements is that they be “confidence distributions.” But, as in (5), this interpretation requires the $\alpha$-cuts of the relative likelihood, i.e., $\{\theta:R(x^{n},\theta)>\alpha\}$, to be $100(1-\alpha)$% confidence sets for $\Theta$, which generally they are not. Therefore, the relative likelihood alone is not enough.

Fortunately, it is at least conceptually straightforward to achieve this calibration by applying what Martin, 2022a calls “validification”—a version of the probability-to-possibility transform (e.g., Dubois et al., 2004; Hose, 2022). In particular, for observed data $X^{n}=x^{n}$, the possibilistic IM’s contour is defined as

and the possibility measure—or upper probability—is

A corresponding necessity measure, or lower probability, is defined via conjugacy as $\underline{\Pi}_{x^{n}}(H)=1-\overline{\Pi}_{x^{n}}(H^{c})$. This measure plays a crucial role in quantifying the support for hypotheses of interest (Cella and Martin, 2023); see Example 5.

An essential feature of this IM construction is its so-called validity property:

This has a number of important consequences. First, (2) is exactly the calibration property required to ensure a “confidence” connection that the relative likelihood on its own misses. In particular, (2) immediately implies that the set

indexed by $\alpha\in[0,1]$, is a $100(1-\alpha)$% frequentist confidence set in the sense that its coverage probability is at least $1-\alpha$:

Moreover, it readily follows that

In words, the IM is valid, or calibrated, if it assigns possibility $\leq\alpha$ to true hypotheses at rate $\leq\alpha$ as a function of data $X^{n}$. This is precisely where the IM’s aforementioned “inferential weight” comes from: (4) implies that we do not expect small $\overline{\Pi}_{x^{n}}(H)$ when $H$ is true, so we are inclined to doubt the truthfulness of those hypotheses $H$ for which $\overline{\Pi}_{x^{n}}(H)$ is small. In addition, the above property ensures that the possibilistic IM is safe from false confidence (Balch et al., 2019; Martin, 2019, 2024), unlike all default-prior Bayes and fiducial solutions. An even stronger, uniform-in-hypotheses version of (4) holds, as shown/discussed in Cella and Martin, (2023):

The “for some $H$ with $H\ni\Theta$” event can be seen as a union of every hypothesis $H$ that contains $\Theta$, which makes it a much broader event than that associated with any single fixed $H$ in (4). Consequently, regardless of the number of hypotheses evaluated or the manner in which they are selected—even if they are data-dependent—the probability that even a single suggestion from the IM is misleading remains controlled at the specified level. For further details about and discussion of the possibilistic IM’s properties, its connection to Bayesian/fiducial inference, etc, see Martin, 2022b ; Martin, 2023b ; Martin, 2023a .

In a Bayesian analysis, inference is based on summaries of the data-dependent posterior distribution, e.g., posterior probabilities of scientifically relevant hypotheses, expectations of loss/utility functions, etc. And all of these summaries boil down to integration involving the probability density function that determines the posterior. Virtually the same statement is true of the possibilistic IM: the lower–upper probability pairs for scientific relevant hypotheses, lower–upper expectations of loss/utility functions, etc. boil down to optimization involving the possibility contour $\pi_{x^{n}}$. For example, if the focus is on a certain feature $\Phi=g(\Theta)$ of $\Theta$, then the Bayesian can find the marginal posterior distribution for $\Phi$ by integrating the posterior density for $\Theta$. Similarly, the corresponding marginal possibilistic IM has a contour that is obtained by optimizing $\pi_{x^{n}}$:

Importantly, unlike with Bayesian integration, the IM’s optimization operation ensures that the validity property inherent in $\pi_{x^{n}}$ is transferred to $\pi_{x^{n}}^{\text{marg}}$, which implies that the IM’s marginal inference about $\Phi$ is safe from false confidence.

The above-described operations applied to the IM’s contour function (1), to obtain upper probabilities or to eliminate nuisance parameters, are all special cases of Choquet integrals (e.g. Troffaes and de Cooman, 2014, App. C). These more general Choquet integrals can be statistically relevant in formal decision-making contexts (Martin, 2021b ), etc. For the sake of space, we will not dig into these details here, but the fact that there are non-trivial operations to be carried out even after the IM’s contour function has been obtained provides genuine, practical motivation for us to find the simplest and most efficient contour approximation as possible.

While the IM construction and properties are relatively simple, computation can be a challenge. The key observation is that we rarely have the sampling distribution of the relative likelihood $R(X^{n},\theta)$, under $\mathsf{P}_{\theta}$, available in closed form to facilitate exact computation of $\pi_{x^{n}}(\theta)$ for any $\theta$. So, instead, the go-to strategy is to approximate that sampling distribution using Monte Carlo at each value of $\theta$ on a sufficiently fine grid (e.g., Martin, 2022b ; Hose et al., 2022). That is, the possibility contour is approximated as

where $1(\cdot)$ is the indicator function and $X_{m,\theta}^{n}$ consists of $n$ iid samples from $\mathsf{P}_{\theta}$ for $m=1,\ldots,M$. The above computation is feasible at one or a few different $\theta$ values, but it is common that this needs to be carried out over a sufficiently fine grid covering the (relevant area) of the often multi-dimensional parameter space $\mathbb{T}$. For example, the confidence set in (3) requires that we can solve the equation $\pi_{x^{n}}(\theta)=\alpha$, or at least find which $\theta$’s satisfy $\pi_{x^{n}}(\theta)\geq\alpha$, and a naive approach is to compute the contour over a huge grid and then keep those that (approximately) solve the aforementioned equation. This amounts to lots of wasted computations. More generally, the relevant summaries of the IM output involve optimization of the contour, and carrying out this optimization numerically requires many contour function evaluations. Simple tweaks to this most-naive approach can be employed in certain cases, e.g., importance sampling, but these adjustments require problem-specific considerations and they cannot be expected to offer substantial improvements in computational efficiency. This is a serious bottleneck, so new and not-so-naive computational strategies are desperately needed.

The Monte Carlo-based strategy reviewed above is not the only way to numerically approximate the possibilistic IM. Another option is an analytical “Gaussian” approximation (see below) based on the available large-sample theory (Martin and Williams, 2024). The goal here is to strike a balance between the (more-or-less) exact but expensive Monte Carlo-based approximation and the rough but cheap large-sample approximation. To strike this balance, we choose to focus on a specific feature of the possibilistic IM’s output, namely, the confidence sets $C_{\alpha}(x^{n})$ in (3), and choose the approximation such that it at least matches the given confidence set exactly. Our specific proposal resembles the variational approximations that are now widely used in Bayesian analysis, where first a relatively simple family of candidate probability distributions is specified and then the approximation is obtained by finding the candidate that minimizes (an upper bound on) the distance/divergence from the exact posterior distribution. Our approach differs in the sense that we are aiming to approximate a possibility measure by (the probability-to-possibility transform applied to) an appropriately chosen probability distribution.

Following Destercke and Dubois, (2014), Couso et al., (2001), and others, the possibilistic IM’s credal set, $\mathscr{C}(\overline{\Pi}_{x^{n}})$, which is the set of precise probability distributions dominated by $\overline{\Pi}_{x^{n}}$, has a remarkably simple and intuitive characterization:

(Where convenient, we will replace the subscript “$x^{n}$” by a subscript “$n$”—e.g., $\mathsf{Q}_{n}$ and $\overline{\Pi}_{n}$ instead of $\mathsf{Q}_{x^{n}}$ and $\overline{\Pi}_{x^{n}}$—to simplify the notation in what follows.) That is, a data-dependent probability $\mathsf{Q}_{n}$ is consistent with $\overline{\Pi}_{n}$ if and only if, for each $\alpha\in[0,1]$, it assigns at least $1-\alpha$ mass to the IM’s confidence set $C_{\alpha}(x^{n})$ in (3). Furthermore, the “best” inner probabilistic approximation of the possibilistic IM, if it exists, corresponds to a $\mathsf{Q}_{n}^{\star}$ such that $\mathsf{Q}_{n}^{\star}\{C_{\alpha}(x^{n})\}=1-\alpha$ for all $\alpha\in[0,1]$. For a certain special class of statistical models, Martin, 2023a showed that this best inner approximation corresponds to Fisher’s fiducial distribution and the default-prior Bayesian posterior distribution. Beyond this special class of models, however, it is not clear if a best inner approximation exists and, if so, how to find it. A less ambitious goal is to find, for a fixed choice of $\alpha$, a probability distribution $\mathsf{Q}_{n,\alpha}^{\star}$ such that

Our goal here is to develop a general strategy for finding, for a given $\alpha\in(0,1)$, a probability distribution $\mathsf{Q}_{n,\alpha}^{\star}$ that (at least approximately) solves the equation in (6). Once identified, we can reconstruct relevant features of the possibilistic IM via (Bayesian-like) Monte Carlo sampling from this $\mathsf{Q}_{n,\alpha}^{\star}$.

We propose to start with a parametric class of data-dependent probability distributions $\mathscr{Q}^{\text{var}}=\{\mathsf{Q}_{n}^{\xi}:\xi\in\Xi\}$ indexed by a generic space $\Xi$. An important example is the case where the $\mathsf{Q}_{n}^{\xi}$’s are Gaussian distributions with mean vector and/or covariance matrix depending on (data and) $\xi$ in some specified way. Specifically, since the possibilistic IM contour’s peak is at the maximum likelihood estimator $\hat{\theta}_{n}=\hat{\theta}_{x^{n}}$, it makes sense to fix the mean vector of the Gaussian $\mathsf{Q}_{n}^{\xi}$ at $\hat{\theta}_{n}$; for the covariance matrix, however, a natural choice would be $\xi^{2}\,J_{n}^{-1}$, where $J_{n}=J_{x^{n}}$ is the observed Fisher information matrix depending on data and on the posited statistical model. This Gaussian family is a very reasonable default choice of $\mathscr{Q}^{\text{var}}$ given that $\mathsf{Q}_{n}^{\xi}$ with $\xi=1$ is asymptotically the best inner probabilistic approximation to $\overline{\Pi}_{n}$; see Martin and Williams, (2024). So, our proposal is to insert some additional flexibility in the Gaussian approximation by allowing the spread to expand or contract depending on whether $\xi>1$ or $\xi<1$. While a Gaussian-based approximation is natural, choosing $\mathscr{Q}^{\text{var}}$ to be a Gaussian family is not the only option for $\mathscr{Q}^{\text{var}}$; see Example 4 below. Indeed, if the parameter space has structure that is not present in the usual Euclidean space where the Gaussian is supported, e.g., if $\mathbb{T}$ is the probability simplex, then choosing $\mathscr{Q}^{\text{var}}$ to respect that structure makes perfect sense.

The one high-level condition imposed on $\mathscr{Q}^{\text{var}}$ is that it be sufficiently flexible, i.e., as $\xi$ varies, the $\mathsf{Q}_{n,\alpha}^{\xi}$-probability of the possibilistic IM’s $\alpha$-cut takes values smaller and larger than the target level $1-\alpha$. This clearly holds for the Gaussian approximation described above, since $\xi$ controls the spread of $\mathsf{Q}_{n,\alpha}^{\xi}$ to the extent that the aforementioned $\mathsf{Q}_{n,\alpha}^{\xi}$-probability can be made arbitrarily small or large by taking $\xi$ sufficiently small or large, respectively. This mild condition can also be readily verified for virtually any other (reasonable) approximation family $\mathscr{Q}^{\text{var}}$.

Given such a suitable choice of approximation family $\mathscr{Q}^{\text{var}}$, indexed by a parameter $\xi\in\Xi$, our proposed procedure is as follows. Define an objective function

so that solving (6) boils down to finding a root of $f_{\alpha}$. If the probability on the right-hand side of (7) could be evaluated in closed-form, then one could apply any of the standard root-finding algorithms to solve this, e.g., Newton–Raphson. However, this $\mathsf{Q}_{n,\alpha}^{\xi}$-probability typically cannot be evaluated in closed-form so, instead, $f_{\alpha}$ can be approximated via Monte Carlo with $\hat{f}_{\alpha}$ defined as

where $\Theta_{1}^{\xi},\ldots,\Theta_{M}^{\xi}\overset{\text{\tiny iid}}{\,\sim\,}\mathsf{Q}_{x^{n}}^{\xi}$. Presumably, the aforementioned samples are cheap for every $\xi$ because the family $\mathscr{Q}$ has been specified by the user. But only having an unbiased estimator of the objective function requires some adjustments to the numerical routine. In particular, rather than a Newton–Raphson routine that assumes the function values are noiseless, we must use a stochastic approximation algorithm (e.g., Syring and Martin, 2021, 2019; Martin and Ghosh, 2008; Kushner and Yin, 2003; Robbins and Monro, 1951) that is adapted to noisy function evaluations of $\hat{f}_{\alpha}$. The basic Robbins–Monro algorithm seeks the root of (7) through the updates

where “$\pm$” depends on whether $\xi\mapsto f_{\alpha}(\xi)$ is decreasing or increasing, $\xi^{(0)}$ is an initial guess, and $(w_{t})$ is a deterministic step size sequence that satisfies

Pseudocode for our proposed approximation is given in Algorithm 1. To summarize, we are proposing a simple data-dependent probability distribution for $\Theta$ whose probability-to-possibility contour closely matches the IM’s contour at least at the specified threshold $\alpha$. More specifically, a reasonable choice of probability distribution is a Gaussian, and the approach presented in Algorithm 1 scales the covariance matrix so that its corresponding contour function accurately approximates the IM’s contour at threshold $\alpha$.

It is well-known that, under certain mild conditions, the sequence $(\xi^{(t)}:t\geq 0)$ as defined above converges in probability to the root of $f_{\alpha}$ in (7). If $\hat{\xi}$ is the value returned when the algorithm reaches practical convergence, e.g., when $|\hat{f}_{\alpha}(\xi^{(t)})|$ or the change $|\xi^{(t+1)}-\xi^{(t)}|$ is smaller than some specified threshold, then we set $\widehat{\mathsf{Q}}_{n,\alpha}=\mathsf{Q}_{n,\alpha}^{\hat{\xi}}$. This distribution should be at least a decent approximation of the IM possibility measure’s inner approximation, i.e., the “most diffuse” member of $\mathscr{C}(\overline{\Pi}_{x^{n}})$. Consequently, the probability-to-possibility transform in (1) applied to $\widehat{\mathsf{Q}}_{x^{n},\alpha}$ should be a decent approximation of the exact possibilistic IM contour $\pi_{x^{n}}$, at least in terms of their upper $\alpha$-cuts. The illustrations presented in Section 4 below show that our proposed approximation is much better than “decent” in cases where we can visually compare with the exact IM contour.

As emphasized above, a suitable choice of $\mathscr{Q}^{\text{var}}$ depends on context. An important consideration in this choice is the ability to perform exact computations of the approximate contour determined by the inner approximation $\widehat{\mathsf{Q}}_{n,\alpha}$. This is the case for the aforementioned normal variational family $\mathscr{Q}^{\text{var}}$, with mean $\hat{\theta}_{n}$ and covariance $\hat{\xi}^{2}J_{n}^{-1}$, as

where $G_{d}$ is the ${\sf ChiSq}(d)$ distribution function. This closed-form expression makes summaries of the contour, which are approximations of certain features (e.g., Choquet integrals) of the possibilistic IM $\overline{\Pi}_{x^{n}}$, relatively easy to obtain numerically.

Our first goal here is to provide a proof-of-concept for the proposed approximation. Towards this, we present a few low-dimensional illustrations where we can visualize both the exact and approximation IM contours and directly assess the quality of the approximation. In all but Example 4 below, we use the normal variational family $\mathscr{Q}$ with mean $\hat{\theta}_{n}$ and covariance $\xi^{2}J_{n}^{-1}$ as described above, with $\xi$ to be determined. All of the examples display the variational-like IM approximation $\widehat{\mathsf{Q}}_{n,\alpha}$ based on $\alpha=0.1$, $M=200$ Monte Carlo samples, step sizes $w_{t}=2(1+t)^{-1}$, and convergence threshold $\varepsilon=0.005$.

Our second goal in this section is to provide a more in-depth illustration of the kind of analysis that can be carried out with the proposed variational-like approximate IM. We do this in the context of regression modeling for count data.

Various tweaks to the basic procedure described in Section 3 might be considered in order applications to speed up computations. If there is a computational bottleneck remaining in the proposal from Section 3, then it would have to be that the IM’s possibility contour must be evaluated at $M$-many points in each iteration of the stochastic approximation. While this is not impractically expensive in some applications, including those presented above, it could be more of a burden in other applications. A related challenge is that a scalar index $\xi$ on the variational family, as we have focused on exclusively so far, may not be sufficiently flexible. Increasing the flexibility by considering higher-dimensional $\xi$ would likewise increase the computational burden, so care is needed. Here we aim to address both of the aforementioned challenges.

The particular modification that we consider here is best suited for cases where the variational family $\mathscr{Q}^{\text{var}}$ satisfies the following property: for each $\xi\in\Xi$, the $100(1-\alpha)$% credible set corresponding to $\mathsf{Q}_{n}^{\xi}$ can be expressed (or at least neatly summarized) in closed-form. The idea to be presented here is more general, but, to keep the details as simple and concrete as possible, we will focus our presentation on the case of a Gaussian variational family. Then the credible sets are ellipsoids in $d$-dimensional space.

As a generalization of the scalar-$\xi$-indexed Gaussian family introduced in Section 3, let $\xi\in\Xi=\mathbb{R}_{\geq 0}^{d}$ be a $d$-vector index and take $\mathscr{Q}^{\text{var}}$ to be the $d$-dimensional Gaussian family with mean vector $\hat{\theta}_{n}$, the maximum likelihood estimator, and covariance matrix $J_{n}(\xi)^{-1}$ defined as follows: take the eigen-decomposition of the observed Fisher information matrix $J_{n}$ as $J_{n}=U\Psi U^{\top}$, and then set

That is, the components of $\xi$ act as multipliers on the singular values of $J_{n}^{-1}$. This indeed generalizes our previously-considered Gaussian approximation family since, if $\xi$ is a scalar, then we recover the previous covariance matrix: $J_{n}(\xi)=J_{n}/\xi^{2}$. Let $K_{n,\alpha}^{\xi}$ denote the $100(1-\alpha)$% credible set corresponding to the distribution $\mathsf{Q}_{n,\alpha}^{\xi}\in\mathscr{Q}^{\text{var}}$. In the present $d$-dimensional Gaussian case, this is an ellipsoid in $\mathbb{R}^{d}$, i.e.,

where $\chi_{d}^{2}(p)$ is the $p^{\text{th}}$ quantile of the ${\sf ChiSq}(d)$ distribution. Then we know that $\mathsf{Q}_{n,\alpha}^{\xi}$ assigns at least probability $1-\alpha$ to the IM’s confidence set $C_{\alpha}(x^{n})$ if

or, equivalently, if

For the Gaussian family of approximations we are considering here, the mode of $\pi_{n}$, the maximum likelihood estimator, is in the interior of each credible set, so the behavior of $\pi_{n}$ outside $K_{n,\alpha}^{\xi}$ is determined by its behavior on the boundary. Moreover, since equality in the above display implies a near-perfect match between the IM’s confidence sets and the posited Gaussian credible sets, a reasonable goal is to find a root to the function

where $\text{bdry}(A)$ denotes the boundary of a set $A$. In designing an iterative algorithm that makes steps towards this root, it is desirable to allow for steps of different magnitudes depending on the direction, and this requires care. For our present Gaussian approximation family, a natural—albeit minimal—representation of the boundary of the credible region $K_{n,\alpha}^{\xi}$ is given by the collection of vectors

where $(\psi_{s},u_{s})$ is the eigenvalue–eigenvector pair corresponding to the $s^{\text{th}}$ largest eigenvalue. Then define the vector-valued function $g_{\alpha}:\Xi\to\mathbb{R}^{d}$ with components

At least intuitively, $\hat{g}_{\alpha,s}(\xi)$ being less than or greater than 0 determines if the credible set $K_{n,\alpha}^{\xi}$ is too large or too small, respectively in the $u_{s}$-direction. From here, we can apply the same stochastic approximation procedure to determine a sequence $(\xi^{(t)}:t\geq 1)$, this time of $d$-vectors, i.e.,

that converges to the root of $\hat{g}_{\alpha}$ which we expect to be close to the root of $g_{\alpha}$. The details of this new but familiar strategy are presented in Algorithm 2. The chief advantage of this new strategy compared to that summarized in Algorithm 1 is computational savings: each iteration of stochastic approximation in Algorithm 2 only requires (Monte Carlo) evaluation of the IM contour at $2d$ many points, which would typically be far fewer than then number of Monte Carlo samples $M$ in Algorithm 1 needed to get a good approximation of $\mathsf{Q}_{n,\alpha}^{\xi}\{\pi_{n}(\Theta)>\alpha\}$.

We will illustrate this new version of the Gaussian variational family and approximation algorithm with two examples. The first is in the relatively simple two-parameter gamma model; the second is a relatively high-dimensional model involving a penalty, intended to serve a gateway into IM solutions to high-dimensional problems.