We study Gaussian Process Thompson Sampling (GP-TS) for sequential decision-making over compact, continuous action spaces and provide a frequentist regret analysis based on fractional Gaussian process posteriors, without relying on domain discretization as in prior work. We show that the variance inflation commonly assumed in existing analyses of GP-TS can be interpreted as Thompson Sampling with respect to a fractional posterior with tempering parameter $α\in (0,1)$. We derive a kernel-agnostic regret bound expressed in terms of the information gain parameter $γ_t$ and the posterior contraction rate $ε_t$, and identify conditions on the Gaussian process prior under which $ε_t$ can be controlled. As special cases of our general bound, we recover regret of order $\tilde{\mathcal{O}}(T^{\frac{1}{2}})$ for the squared exponential kernel, $\tilde{\mathcal{O}}(T^{\frac{2ν+3d}{2(2ν+d)}} )$ for the Matérn-$ν$ kernel, and a bound of order $\tilde{\mathcal{O}}(T^{\frac{2ν+3d}{2(2ν+d)}})$ for the rational quadratic kernel. Overall, our analysis provides a unified and discretization-free regret framework for GP-TS that applies broadly across kernel classes.

Flexible Bayesian models are typically constructed using limits of large parametric models with a multitude of parameters that are often uninterpretable. In this article, we offer a novel alternative by constructing an exponentially tilted empirical likelihood carefully designed to concentrate near a parametric family of distributions of choice with respect to a novel variant of the Wasserstein metric, which is then combined with a prior distribution on model parameters to obtain a robustified posterior. The proposed approach finds applications in a wide variety of robust inference problems, where we intend to perform inference on the parameters associated with the centering distribution in presence of outliers. Our proposed transport metric enjoys great computational simplicity, exploiting the Sinkhorn regularization for discrete optimal transport problems, and being inherently parallelizable. We demonstrate superior performance of our methodology when compared against state-of-the-art robust Bayesian inference methods. We also demonstrate equivalence of our approach with a nonparametric Bayesian formulation under a suitable asymptotic framework, testifying to its flexibility. The constrained entropy maximization that sits at the heart of our likelihood formulation finds its utility beyond robust Bayesian inference; an illustration is provided in a trustworthy machine learning application.

When a statistical model $\{P_θ : θ\in Θ\}$ lacks analytically tractable likelihoods, parametric statistical inference based on data generated from an unknown underlying distribution $P$ can still be performed as long as simulations from the model are possible. This approach is called Simulation Based Inference (SBI). Statistical models are rarely exactly correct (that is, $P \notin \{P_θ: θ\in Θ\}$), and Robust SBI focuses on inferring a reasonable parameter even under model mis-specification. We focus on the setting where $P$ possesses potentially both geometric and Total Variation type discrepancies from $P_{θ^*}$. For this problem, we use a Kullback-Liebler informed robust Optimal Transport divergence, motivated by Empirical Likelihood considerations. We introduce a stochastic sub-gradient ascent algorithm with a convergence guarantee for estimating the semi-discrete version of this robust Optimal Transport divergence, and design a parallelized SBI algorithm which employs the regular bootstrap on top of minimum semi-discrete robust Optimal Transport for parameter uncertainty quantification. We demonstrate mathematically why the divergence is robust under a joint geometric plus Total Variation type contamination and then illustrate the robustness of inferences on a complex benchmark SBI task.

The marginal likelihood or evidence in Bayesian statistics contains an intrinsic penalty for larger model sizes and is a fundamental quantity in Bayesian model comparison. Over the past two decades, there has been steadily increasing activity to understand the nature of this penalty in singular statistical models, building on pioneering work by Sumio Watanabe. Unlike regular models where the Bayesian information criterion (BIC) encapsulates a first-order expansion of the logarithm of the marginal likelihood, parameter counting gets trickier in singular models where a quantity called the real log canonical threshold (RLCT) summarizes the effective model dimensionality. In this article, we offer a probabilistic treatment to recover non-asymptotic versions of established evidence bounds as well as prove a new result based on the Gibbs variational inequality. In particular, we show that mean-field variational inference correctly recovers the RLCT for any singular model in its canonical or normal form. We additionally exhibit sharpness of our bound by analyzing the dynamics of a general purpose coordinate ascent algorithm (CAVI) popularly employed in variational inference.

Many statistical models involve parameter-dependent normalizing constants that are computationally intractable, creating substantial obstacles to standard Bayesian inference. Although existing likelihood-based algorithms can often circumvent these constants, their uncertainty quantification may be poorly calibrated under model misspecification. To address these challenges, we propose SME-BETEL, a semiparametric Bayesian framework that combines score matching estimating equations with Bayesian exponentially tilted empirical likelihood. The resulting posterior avoids evaluation of normalizing constants and does not require learning-rate calibration. We establish consistency and asymptotic normality of the score matching estimator, and prove a Bernstein-von Mises theorem for the SME-BETEL posterior. These results show that SME-BETEL credible sets are asymptotically calibrated to the sampling variability of the score matching estimator, yielding valid frequentist coverage under model misspecification. We further develop a new score matching criterion for mixed-domain data, extending SME-BETEL to models whose observations combine components from different sample spaces. This construction enables robust Bayesian inference for mixed-domain doubly-intractable models, including preferential sampling models with an intractable spatial normalizing integral. Simulation studies show that SME-BETEL remains competitive under correct specification and substantially improves uncertainty quantification under misspecification. An ozone-monitoring application further demonstrates the utility of the mixed-domain construction for spatial preferential modeling.

We study a variant of the Thompson Sampling (TS) algorithm, called α-TS, for solving stochastic generalized linear bandit problems. Existing analyses of TS require inflating the posterior variance to derive near-optimal regret guarantees. We formalize the idea of variance inflation by introducing α-TS that uses a fractional or α-posterior instead of the standard posterior. Our main contribution is to identify general regularity conditions on the prior and reward distributions that enable a regret analysis of α-TS without assuming any tractable approximation of the posterior distribution, unlike previous works. For a specific choice of α∝ d^-1, our general regret bound yields the best known regret bound of O(d^3/2√(T)log T) for both the exponential and sub-Gaussian families of reward distributions. We further provide an α-dependent lower bound showing that the regret constant depends on the product αd, and that when α∝ d^-1 the regret scales as Ω(d^3/2√(T)), explaining the origin of the d^3/2 factor in the upper bound. Our proof technique adapts and combines recent advancements in the analysis of linear bandit problems with first- and second-order posterior concentration theory from the Bayesian statistics literature.

We propose a semiparametric Bayesian methodology for estimating the average treatment effect (ATE) within the potential outcomes framework using observational data with high-dimensional nuisance parameters. Our method introduces a Bayesian debiasing procedure that corrects for bias arising from nuisance estimation and employs a targeted modeling strategy based on summary statistics rather than the full data. These summary statistics are identified in a debiased manner, enabling the estimation of nuisance bias via weighted observables and facilitating hierarchical learning of the ATE. By combining debiasing with sample splitting, our approach separates nuisance estimation from inference on the target parameter, reducing sensitivity to nuisance model specification. We establish that, under mild conditions, the marginal posterior for the ATE satisfies a Bernstein-von Mises theorem when both nuisance models are correctly specified and remains consistent and robust when only one is correct, achieving Bayesian double robustness. This ensures asymptotic efficiency and frequentist validity. Extensive simulations confirm the theoretical results, demonstrating accurate point estimation and credible intervals with nominal coverage, even in high-dimensional settings. The proposed framework can also be extended to other causal estimands, and its key principles offer a general foundation for advancing Bayesian semiparametric inference more broadly.

There is growing interest in developing statistical estimators that achieve exponential concentration around a population target even when the data distribution has heavier than exponential tails. More recent activity has focused on extending such ideas beyond Euclidean spaces to Hilbert spaces and Riemannian manifolds. In this work, we show that such exponential concentration in presence of heavy tails can be achieved over a broader class of parameter spaces called CAT(κ) spaces, a very general metric space equipped with the minimal essential geometric structure for our purpose, while being sufficiently broad to encompass most typical examples encountered in statistics and machine learning. The key technique is to develop and exploit a general concentration bound for the Fréchet median in CAT(κ) spaces. We illustrate our theory through a number of examples, and provide empirical support through simulation studies.

We consider the problem of learning the underlying causal structure among a set of variables, which are assumed to follow a Bayesian network or, more specifically, a linear recursive structural equation model (SEM) with the associated errors being independent and allowed to be non-Gaussian. A Bayesian hierarchical model is proposed to identify the true data-generating directed acyclic graph (DAG) structure where the nodes and edges represent the variables and the direct causal effects, respectively. Moreover, incorporating the information of non-Gaussian errors, we characterize the distribution equivalence class of the true DAG, which specifies the best possible extent to which the DAG can be identified based on purely observational data. Furthermore, under the consideration that the errors are distributed as some scale mixture of Gaussian, where the mixing distribution is unspecified, and mild distributional assumptions, we establish that by employing a non-standard DAG prior, the posterior probability of the distribution equivalence class of the true DAG converges to unity as the sample size grows. This shows that the proposed method achieves the posterior DAG selection consistency, which is further illustrated with examples and simulation studies.

There is extensive literature on accelerating first-order optimization methods in an Euclidean setting. Under which conditions such acceleration is feasible in Riemannian optimization problems is an active area of research. Motivated by the recent success of silver stepsize methods in the Euclidean setting, we undertake a study of such algorithms in the Riemannian setting. We provide the new class of algorithms determined by the choice of vector transport that allows the silver stepsize acceleration on Riemannian manifolds for the function classes associated with the corresponding vector transport. As a core application, we show that our algorithm recovers the standard Wasserstein gradient descent on the 2-Wasserstein space and, as a result, provides the first provable accelerated gradient method for potential functional optimization problems in the Wasserstein space. In addition, we validate the numerical strength of the algorithm for standard benchmark tasks on the space of symmetric positive definite matrices.

Mean-field approximations are widely used for efficiently approximating high-dimensional integrals. While the efficacy of such approximations is well understood for well-behaved likelihoods, it is not clear how accurately it can approximate the marginal likelihood associated with a highly non log-concave singular model. In this article, we provide a case study of the convergence behavior of coordinate ascent variational inference (CAVI) for a general d-dimensional singular model in standard form. We prove that such a model with real log canonical threshold (RLCT). and multiplicity m, the CAVI system converges to one of m locally attracting fixed points. Furthermore, at each of these fixed points, the evidence lower bound (ELBO) recovers the leading-order behavior of the asymptotic expansion of the log marginal likelihood predicted by Watanabe [1, 2, 3]. Our empirical results demonstrate that for models with multiplicity m = 1 the ELBO provides a tighter approximation to the log-marginal likelihood than Watanabe's [1] asymptotic approximation - lambda log n + o(log log n) for a range of sample sizes.

We consider the problem of recovering the true causal structure among a set of variables, generated by a linear acyclic structural equation model (SEM) with the error terms being independent, not necessarily Gaussian, and having equal variances. It is well-known that the true underlying directed acyclic graph (DAG) encoding the causal structure is uniquely identifiable under this assumption. Interestingly, in this setting, it further holds that the sum of minimum expected squared errors for every variable, while predicted by the best linear combination of its parent variables, is minimised if and only if the causal structure is represented by any supergraph of the true DAG. In this work, we propose a Bayesian DAG selection method, where the working model assumes Gaussian SEM with equal error variances, and employ independent g-priors on each set of SEM coefficients. Furthermore, we utilise the aforementioned key property to establish that the proposed method recovers the true graph consistently without any additional distributional assumption, and illustrate it with a simulation study.

The power spectrum of biomedical time series provides important indirect measurements of physiological processes underlying health and biological functions. However, simultaneously characterizing power spectra for multiple time series remains challenging due to extra spectral variability and varying time series lengths. We propose a method for hierarchical Bayesian estimation of stationary time series (HBEST) that provides an interpretable framework for efficiently modeling multiple power spectra. HBEST models log power spectra using a truncated cosine basis expansion with a novel global-local coefficient decomposition, enabling simultaneous estimation of population-level and individual-level power spectra and accommodating time series of varying lengths. The fully Bayesian framework provides shrinkage priors for regularized estimation and efficient information sharing. Simulations demonstrate HBEST's advantages over competing methods in computational efficiency and estimation accuracy. An application to heart rate variability time series demonstrates HBEST's ability to accurately characterize power spectra and capture associations with traditional cardiovascular risk factors.

Finite-rank approximations are widely used to scale Gaussian process (GP) regression, but their posterior behavior can differ from that of the corresponding parent GP prior. We study a class of finite-rank GP priors built from locally supported basis expansions with dependent Gaussian coefficients. Our framework covers finite-element approximations based on the stochastic partial differential equation (SPDE) representation of Matérn GPs and regular-grid GP interpolation schemes. We show that, with a suitable prior on the resolution parameter N, these finite-rank expansions inherit the same posterior contraction rate as the corresponding parent GP prior under the same bandwidth specification used for that parent prior. Consequently, the interpolation construction under a squared-exponential parent GP attains the minimax-optimal rate up to logarithmic factors under a hierarchical prior on the bandwidth parameter and on N, while the SPDE construction attains the same rate under a bandwidth scaling depending on the sample size and the smoothness of the true function, together with a prior on N. We also develop a posterior sampler for the hierarchical interpolation model that jointly updates the resolution and bandwidth parameters, and we provide numerical studies that support the theory.

The orbital superposition method originally developed by Schwarzschild (1979) is used to study the dynamics of growth of a black hole and its host galaxy, and has uncovered new relationships between the galaxy's global characteristics. Scientists are specifically interested in finding optimal parameter choices for this model that best match physical measurements along with quantifying the uncertainty of such procedures. This renders a statistical calibration problem with multivariate outcomes. In this article, we develop a Bayesian method for calibration with multivariate outcomes using orthogonal bias functions thus ensuring parameter identifiability. Our approach is based on projecting the posterior to an appropriate space which allows the user to choose any nonparametric prior on the bias function(s) instead of having to model it (them) with Gaussian processes. We develop a functional projection approach using the theory of Hilbert spaces. A finite-dimensional analogue of the projection problem is also considered. We illustrate the proposed approach using a BART prior and apply it to calibrate the Schwarzschild model illustrating how a multivariate approach may resolve discrepancies resulting from a univariate calibration.

In this work, we offer a thorough analytical investigation into the role of shared hyperparameters in a hierarchical Bayesian model, examining their impact on information borrowing and posterior inference. Our approach is rooted in a non-asymptotic framework, where observations are drawn from a mixed-effects model, and a Gaussian distribution is assumed for the true effect generator. We consider a nested hierarchical prior distribution model to capture these effects and use the posterior means for Bayesian estimation. To quantify the effect of information borrowing, we propose an integrated risk measure relative to the true data-generating distribution. Our analysis reveals that the Bayes estimator for the model with a deeper hierarchy performs better, provided that the unknown random effects are correlated through a compound symmetric structure. Our work also identifies necessary and sufficient conditions for this model to outperform the one nested within it. We further obtain sufficient conditions when the correlation is perturbed. Our study suggests that the model with a deeper hierarchy tends to outperform the nested model unless the true data-generating distribution favors sufficiently independent groups. These findings have significant implications for Bayesian modeling, and we believe they will be of interest to researchers across a wide range of fields.

In recent years, it has become increasingly easy to obtain approximate posterior samples via efficient computation algorithms, such as those in variational Bayes. On the other hand, concerns exist on the accuracy of uncertainty estimates, which make it tempting to consider exploiting the approximate samples in canonical Markov chain Monte Carlo algorithms. A major technical barrier is that the approximate sample, when used as a proposal in Metropolis-Hastings steps, tends to have a low acceptance rate as the dimension increases. In this article, we propose a simple yet general solution named ”graph-accelerated Markov Chain Monte Carlo”. We first build a graph with each node location assigned to an approximate sample, then we run Markov chain Monte Carlo with random walks over the graph. In the first stage, we optimize the choice of graph edges to enforce small differences in posterior density/probability between neighboring nodes, while encouraging edges to correspond to large distances in the parameter space. This optimized graph allows us to accelerate a canonical Markov transition kernel through mixing with a large-jump Metropolis-Hastings step, when collecting Markov chain samples at the second stage. Due to its simplicity, this acceleration can be applied to most of the existing Markov chain Monte Carlo algorithms. We theoretically quantify the rate of acceptance as dimension increases, and show the effects on improved mixing time. We demonstrate our approach through improved mixing performances for challenging sampling problems, such as those involving multiple modes, non-convex density contour, or large-dimension latent variables.

We study the estimation of Zipfian distributions under L1 loss, and provide near minimax optimal bounds in several regimes. Specifically, we assume observations arrive from a known alphabet, and with a known decay rate parametrizing the Zipfian, but we do not know a priori which alphabet elements have larger probability than others. We present a novel Sort and Snap estimator, which uses the empirical proportions to sort the alphabet, and then snaps them to the associated term from the Zipfian distribution. For arbitrary decay rates and smaller alphabet sizes, as well as for large decay rates and large alphabet sizes, we show an exact or minor variant of this estimator is near minimax optimal and has exponential improvement over the standard empirical proportion estimator. However, for small decay rates and larger alphabet sizes a simulation study indicates the standard empirical proportion estimator is competitive with Sort and Snap procedures. In addition to providing nearly tight bounds for important high-dimensional estimation problems, we believe the Sort and Snap estimator, and its analysis, will have independent interest.

Inference in semi-supervised (SS) settings has received substantial attention in recent years due to increased relevance in modern big-data problems. In a typical SS setting, there is a much larger sized unlabeled data, containing observations only for a set of predictors, in addition to a moderately sized labeled data containing observations for both an outcome and the set of predictors. Such data arises naturally from settings where the outcome, unlike the predictors, is costly or difficult to obtain. One of the primary statistical objectives in SS settings is to explore whether parameter estimation can be improved by exploiting the unlabeled data. A novel Bayesian approach to SS inference for the population mean estimation problem is proposed. The proposed approach provides improved and optimal estimators both in terms of estimation efficiency as well as inference. The method itself has several interesting artifacts. The central idea behind the method is to model certain summary statistics of the data in a targeted manner, rather than the entire raw data itself, along with a novel Bayesian notion of debiasing. Specifying appropriate summary statistics crucially relies on a debiased representation of the population mean that incorporates unlabeled data through a flexible nuisance function while also learning its estimation bias. Combined with careful usage of sample splitting, this debiasing approach mitigates the effect of bias due to slow rates or misspecification of the nuisance parameter from the posterior of the final parameter of interest, ensuring its robustness and efficiency. Concrete theoretical results, via Bernstein–von Mises theorems, are established, validating all claims, and are further supported through extensive numerical studies. To our knowledge, this is possibly the first work on Bayesian inference in SS settings, and its central ideas also apply more broadly to other Bayesian semi-parametric inference problems.

As a computational alternative to Markov chain Monte Carlo approaches, variational inference (VI) is becoming more and more popular for approximating intractable posterior distributions in large-scale Bayesian models due to its comparable efficacy and superior efficiency. Several recent works provide theoretical justifications of VI by proving its statistical optimality for parameter estimation under various settings; meanwhile, formal analysis on the algorithmic convergence aspects of VI is still largely lacking. In this paper, we consider the common coordinate ascent variational inference (CAVI) algorithm for implementing the mean-field (MF) VI towards optimizing a Kullback--Leibler divergence objective functional over the space of all factorized distributions. Focusing on the two-block case, we analyze the convergence of CAVI by leveraging the extensive toolbox from functional analysis and optimization. We provide general conditions for certifying global or local exponential convergence of CAVI. Specifically, a new notion of generalized correlation for characterizing the interaction between the constituting blocks in influencing the VI objective functional is introduced, which according to the theory, quantifies the algorithmic contraction rate of two-block CAVI. As illustrations, we apply the developed theory to a number of examples, and derive explicit problem-dependent upper bounds on the algorithmic contraction rate.