Papers

Some of the sharpest challenges in sampling from the energy functions of physical systems arise at phase transitions, where the density of states changes abruptly and many sampling algorithms stall. Nested sampling is a particle method that traverses the density of states under a hard energy constraint and is known to be robust to such transitions, but its application in high dimension is limited by the difficulty of sampling under that constraint. In this work we introduce Quenched Ensemble Sampling, which generalises the hard constraint to a family of repulsive potentials at the energy boundary. This preserves the quenched path of monotonically decreasing energy while making the constrained target amenable to scalable gradient-based kernels. We demonstrate on synthetic models of phase transitions that our method estimates the marginal likelihood and draws posterior samples across a first-order transition where popular alternatives such as tempering fail. We apply the procedure to marginal likelihood estimation in Bayesian neural networks, enabling model comparison between network architectures. Finally, in a high-dimensional continuous lattice field theory, we show that this method traverses a first-order transition and estimates the partition function.

Markov chain Monte Carlo is widely used for model assessment and parameter fitting in astronomy and astrophysics, leading to numerous ready-to-use packages implementing various sampling techniques. For large-scale analyses involving many individual targets, computational efficiency becomes paramount, as it dictates the overall runtime. We introduce Nii-MALA, a high-performance C implementation of the Metropolis-Adjusted Langevin Algorithm, which leverages Message Passing Interface for parallelization and incorporates automatic differentiation backends for gradient evaluation. The code offers flexible control via a configuration file, allowing users to specify the proposal step sizes for all model parameters across parallel chains. We validate its effectiveness on 15- and 30-dimensional Gaussian distributions and a bimodal distribution, and further benchmark it against several other codes using radial velocity data from 51 Pegasi b. Our benchmarks confirm the efficiency of Langevin sampling over random-walk sampling in obtaining the effective sample size, particularly as the dimensionality of the target function increases. This algorithmic advantage can be effectively leveraged when paired with an efficient automatic differentiation library. However, for complex functions commonly encountered in astronomy—such as radial velocity orbital fitting—the derivative-dependent Langevin sampler shows lower efficiency than parallel tempering or alternative sampling approaches, as revealed by effective sample size analysis.

We study lazy single-site Metropolis dynamics $P_{\beta}$ at inverse temperature $\beta \geq 0$ on $\{-1,+1\}^d$ for an Ising star with additional signed interactions among the leaves. If the absolute row sums of the leaf-interaction matrix are at most $\kappa\le1/2$, the worst-case total-variation mixing time of $P_{\beta}$ is at least of order $d\exp\{c\beta(d-1)\}$ for $\beta\ge1$, with universal $c>0$. Averaging over global spin reversal reduces the mixing time to $\mathcal O(d^2(1+\beta))$ for both $GP_\beta G$ and $(P_\beta+G)/2$, where $G$ is the Gibbs kernel induced by the partition of the state space into spin-reversal orbits $\{-x,x\}$. Partition-function interpolation gives the lower bound for $P_\beta$. The upper bounds follow from Wu's Dobrushin inequality and a decomposition into projection and restriction chains. This gives an explicit example in which group averaging turns an exponentially slow chain into a polynomially fast one.

The problem of partial identification concerns bounding causal quantities that remain unidentifiable given the observed distribution and the causal diagram of the underlying structural causal model (SCM). While bounds have been developed for certain causal quantities with continuous outcomes, they are necessarily univariate and are unable to reflect differing levels of confidence. Building on previous work in partial identification, we propose a simple Bayesian method for deriving probabilistic bounds on a large class of causal quantities with potentially multivariate and/or continuous outcome variables.

The basic question in perturbation analysis of Markov chains is how small changes in their transition kernels affect their stationary distributions. Classical perturbation bounds typically require the kernel error to be much smaller than $1/\tau$, where $\tau$ is a mixing or relaxation time. Although this scaling is sharp for Markov chains in general, we investigate a general "square-rooting" phenomenon in which one-step errors of order roughly $1/\sqrt{\tau}$ can be sufficient for local updates. We proved a form of this phenomenon in Lin, Liu and Smith (2025) under strong assumptions. Here we substantially weaken these assumptions, and prove this phenomenon occurs using three distinct approaches. First, block factorization applies under structural assumptions on the stationary measures of both chains. Second, approximate block-update arguments extend the result to statistically relevant Markov chain Monte Carlo (MCMC) settings, where structural guarantees are available for the exact posterior and its associated sampler, but not for the perturbed posterior. Third, we use direct calculations for a class of models with hard constraints where neither general result is directly available. We illustrate these results in three MCMC settings and show how they directly inform the tuning of approximate MCMC algorithms.

We study the mixing time of a systematic-scan analogue of Kac's walk that was proposed as a fast surrogate for Haar-distributed orthogonal matrices in randomized high-dimensional algorithms and was conjectured to approach Haar measure after only logarithmically many sweeps. We show that this conjectured speed-up does not occur for convergence of the full matrix law to Haar measure in Frobenius Wasserstein distance. At fixed normalized accuracy, the mixing time lies between order $n/\log n$ and order $n$ sweeps; at fixed absolute Frobenius accuracy, the corresponding bounds are between order $n$ and order $n\log n$. More strongly, below the scale $n/\log n$, the normalized Wasserstein distance remains asymptotically at its extremal value. We also show that the output law is singular with respect to Haar measure for fewer than $n/2$ sweeps. Thus the sampler may provide effective application-specific randomization without exhibiting the much faster full-Haar mixing.

Computing matrix inverses (or, e.g., left inverses) using known algorithms is in general computationally prohibitive in high dimension, such as in an operational weather forecasting context. In this article, a swift, scalable algorithm to calculate the left inverse (through singular values decompositions techniques) of a square root of a horizontally localized ensemble covariance matrix is introduced. In this algorithm, calculations involving model variables located in different grid vertical columns can be performed separately. This is the aspect that makes this algorithm scalable by allowing for parallel computing. The algorithm is developed in the case when the number of ensemble members equals the number of model variables in each grid vertical column (this number is of order $10^3$ for current high-resolution models). The algorithm is tested in a synthetic experiment on a three-dimensional grid, and it is found to have very good accuracy. Additionally, a version of the algorithm outputting the left-multiplication of a generic vector by that above-mentioned left inverse is also introduced. For operational applications, if executed on a single processor, that version of the algorithm is already $9$ orders of magnitude faster than a traditional algorithm, and, with multiple, parallel processors, it scales until the number of processors reaches the order of $10^6$. In addition, an application of the algorithm is pointed out: performing the transform from model variables to an alternative set of variables, defined herein as standardized variables, being climatologically uncorrelated and each having unit climatological error variance. Furthermore, it is pointed out that the standardized variables allow constructing a new Hybrid covariance model, which is positive definite and climatologically unbiased, and it features a geographically dependent climatological error covariance matrix.

This note presents new findings on the distribution of the volume and surface area of both the Poisson polyhedron and the Crofton polyhedron. It also considers the distribution of the area of the Poisson polygon and the Crofton polygon in the plane. This note motivates the modeling of the distributions of these convex body functionals as mixture distributions and the essential inclusion of the generalized gamma distribution therein as mixture components. Finally, model parameters are estimated and models are selected based on extensive stochastic simulations of these random polytopes. This includes that for the first time, this note gives very accurate estimates of the side number probabilities for both the Poisson polyhedron and the Crofton polyhedron, as these are used as mixing weights in the distribution models.

Medical image classification is frequently complicated by transitional categories whose feature distributions overlap those of adjacent classes, producing ambiguous decision boundaries. Conformal prediction returns uncertainty-aware prediction sets, but these are not directly actionable in clinical screening, where a single decision is required. This work proposes adaptive conformal redistribution (AdaConRed), a label-free post-conformal decision rule that converts ambiguous prediction sets into refined class assignments. A five-stage pipeline is developed. Vision-language generative augmentation addresses minority-class scarcity; a frozen DermFoundation encoder provides embeddings; a lightweight multi-layer perceptron performs classification; an entropy-modulated, margin-aware nonconformity score constructs adaptive prediction sets; samples predicted as transitional with multi-label sets are reassigned to the most probable alternative class within the set, using only model outputs at inference. Evaluation uses the OSCC oral lesion and ISIC skin lesion benchmarks at a miscoverage level of 0.2. On the 3-class OSCC benchmark, overall accuracy improves from 73.54% to 77.38%, with oral cancer accuracy rising from 64.29% to 82.14% and benign accuracy from 56.57% to 70.20%. Reassignment of transitional samples reduces OPMD accuracy from 84.78% to 80.16%, consistent with the asymmetric cost of missed malignancy. On ISIC, overall accuracy improves from 85.83% to 87.19%, melanoma accuracy rising from 66.04% to 68.34%. AdaConRed outperforms LAC, APS and RAPS under an identical backbone and redistribution rule. Conformal prediction can be extended beyond uncertainty quantification toward actionable decision support where transitional disease categories are present, with gains concentrated in the clinically critical malignant categories. Code repository: https://github.com/saibal436ghosh/AdaConRed.

We demonstrate explicit transformations of three independent uniforms with exact Gamma$(a,1)$, $a>0$, Poisson$(\lambda)$, $\lambda>0$, and Binomial$(N,p)$, $N\geq 1$, $0\leq p \leq 1$ laws using only elementary operations.

Accurate estimation of aboveground biomass density is essential for quantifying forest carbon stocks. NASA's GEDI mission provides valuable canopy structure data, but its sparse sampling necessitates the use of simulators to calibrate biomass models at field inventory locations. The widely used simulator of Hancock et al. (2019) emulates GEDI waveforms from airborne LiDAR point clouds, yet it has never been validated over European temperate forests. Here, we compare approximately 9,500 pairs of observed and simulated GEDI relative height (RH) profiles across French forests using the national airborne LiDAR program as input. We separate two sources of error: waveform modeling differences, assessed by referencing both simulated and real RH metrics to a common ALS-derived ground elevation, and ground detection bias, evaluated by comparing each GEDI L2A processing algorithm against the ALS reference. Under the baseline configuration, the mean absolute bias across the full RH profile reaches 0.69 m in leaf-on and 1.28 m in leaf-off conditions. Switching to intensity-based return weighting and selecting the a3 L2A algorithm reduces these biases to 0.44 m and 0.40 m respectively. The a3 algorithm also achieves near-unbiased ground detection (-0.03 m versus -0.88 m for the default), directly reducing a previously overlooked source of error. We also show that leaf-off acquisitions and low-sensitivity shots, both typically excluded from standard biomass products, are simulated as reliably as their counterparts, substantially expanding the potential calibration and inference datasets.

Cluster-based permutation tests are widely used for analyzing MEG data, even though they are limited to cluster-level inference and do not provide spatially resolved effect estimates. We propose a regression-based framework for modeling brain activity directly on the cortical surface. Spatial effects are represented using Wendland radial basis functions, and model-based gradient boosting is employed for data-driven selection and estimation of localized activation patterns. This yields interpretable, spatially resolved effect estimates while mitigating the need for extensive multiple testing correction. In a simulation study with heterogeneous signal structures, the proposed approach recovers localized effects that are difficult to detect using cluster-based methods. An application to experimental MEG data further illustrates its ability to reveal spatially specific activation patterns. Overall, this unified framework provides a flexible and interpretable approach for modeling cortical surface data within a unified statistical setting.

Constructing a confidence interval for the difference between two independent binomial proportions involves a nuisance direction that is not identified by the estimand. The one-sample Wilson score interval inverts a scalar score test, but has no direct bivariate analogue isolating the difference: inverting the joint normal approximation yields an elliptical region in the unit square, whereas the estimand \(p_1-p_2\) is one-dimensional. We define the Elliptically Optimal (EO) confidence interval as the range of \(p_1-p_2\) over this region and solve the resulting optimization problem in closed form, obtaining explicit bounds in six mutually exclusive and exhaustive cases. The solution admits a compact characterization: the EO interval is the score interval obtained by maximizing over the nuisance variance rather than estimating it. It is therefore the shortest interval obtained by projecting the elliptical region, and inherits its coverage guarantee. We derive the exact coverage excess, \(2[\Phi(z\mathcal{R})-\Phi(z)]\), where \(\mathcal{R}\) is the ratio of the least-favourable to the true standard deviation. The excess vanishes on an explicit line through the parameter space, is bounded by \(\alpha\), and is invariant under proportional scaling of the sample sizes. Exact enumeration of the binomial coverage shows that the Wald interval, whose variance estimator is downward biased by a factor \(1-1/n\) under balanced allocation, falls below nominal coverage almost everywhere. The EO interval never under-covers under the normal approximation and always yields admissible, non-degenerate bounds. Its price is over-coverage when both proportions are extreme, which we quantify exactly.

Recovering the linear relationships that govern a system from noisy measurements is a basic task across the physical and engineering sciences. Because every measured variable may carry an unknown amount of noise, classical regression must commit in advance to a set of structural assumptions: ordinary least squares requires a declared input-output partition with input variables being noise-free, total least squares assumes equal noise variance across all variables, and generalized total least squares additionally requires the noisy-variable partition and variances to be known beforehand. Kalman~\cite{Kalman:1982} showed that any procedure returning a unique linear model from inexact data must rest on such unverifiable a priori assumptions – “prejudices” – that cannot be checked against the data itself, and that removing them leaves the identification problem fundamentally indeterminate. Whether these prejudices can instead be resolved directly from the data has remain unresolved. Here we show that an iterative generalized-eigenvalue algorithm, QZ-IPCA, recovers the noisy-variable partition, noise variances, number of linear relations, and regression coefficients of a multivariate linear system simultaneously, using only the raw data. Across all possible exhaustive noise configurations of a five-variable benchmark network, QZ-IPCA correctly identifies model structure and recovers coefficients with error below 6.4%. It outperforms ordinary least squares even when given the best partition, and succeeds in rank identification precisely where standard total least squares falls once noise variances differ across variables. These results show that the assumptions conventionally required for multivariate regression are not necessary, recasting model identification as a problem solvable from data geometry alone.

We consider the problem of estimating a high-dimensional covariance matrix from a very limited number of samples. This problem is ubiquitous in computational fluid dynamics, where a small number of fluid snapshots must be used to construct a Gramian matrix determining a reduced-order model, as well as in computational geoscience, where a small ensemble of Earth system forecasts must be used to estimate the covariance matrix associated with the forecast uncertainty. It is common practice to regularize the small-sample covariance by imposing a "localization" structure that enforces a physically realistic correlation length scale, imposing a sparsity constraint, "shrinking" towards a prescribed target, or attenuating small correlations. We propose an alternate technique that regularizes the small-sample covariance by imposing hierarchical rank structure. Compared to regularization methods that assume sparsity such as spatial localization, hierarchical rank structure accommodates a wider range of covariance matrices, roughly corresponding to situations where long-range correlations vary more smoothly than short-range ones. It also results in a data-sparse matrix format that permits highly efficient matrix-vector products. We present theory and algorithms which show how to efficiently estimate a high-dimensional, hierarchically rank structured covariance matrix from limited samples. Through an error analysis and numerical experiments with a variety of model problems, we demonstrate that these techniques are effective at reducing sampling errors, and that in many cases they achieve smaller estimation error than conventional techniques.

A latent space network model places the nodes in a metric space and lets the probability of a tie decrease with distance. In a space of constant curvature, pairwise distances determine the positions up to an isometry. In the products and in the three remaining three-dimensional model geometries they do not. We study the three geometries that are neither of constant curvature nor products: the Heisenberg group, the solvable group and the universal cover of the unit tangent bundle of the hyperbolic plane. We ask what one network identifies about positions in them and when their geometry is detectable. Two anchors remove the isometry ambiguity. Small configurations are not determined by their distances, and generic local identification holds beyond a finite threshold, certified in the Heisenberg group. We derive the posterior on the quotient by the isometry group. For small configurations, the divergence to the nearest product or constant-curvature competitor is the stress component of the curvature difference and vanishes at high order in the scale; undirected ties therefore detect the geometry only in large networks with large enclosed areas. Directed ties expose it at first order: in the two twisted geometries, asymmetric preferences circulate around triangles in proportion to enclosed area, which no additive ranking produces. On dense competitive-game counter networks, the coupled model beats rankings on every geometry, degree-corrected rankings and free antisymmetric terms. A Euclidean model with the same coupled term matches it, so the gain is the coupling of similarity and circulation through shared coordinates. An additive-and-multiplicative-effects model predicts better still.

As an extension of the Point Estimate Method (PEM) to evaluate probabilistic moments of quantities of interest (QoI) in general $n$-dimensional spaces, the Quadratic Point Estimate Method (QPEM) has been recently developed. This new method is defined to fully represent up to fifth-order input moments in the Gaussian space, providing general analytical expressions for sample locations and weights, without requiring any numerical optimization. The QPEM can significantly improve the estimation accuracy of the output QoI moments, in relation to PEM-based methods whose numbers of sigma points grow linearly with the problem dimension, while at the same time having an affordable and competitive computational cost up to a considerable number of dimensions. The QPEM is further enhanced in this work by enabling copula integration into the framework, which enables effective modeling of the joint input probability density function by estimating marginals and the dependence structure of the involved random variables. The validity and efficient performance of the copula-based QPEM are showcased against numerous other sampling methods in various examples considering two practical scenarios: (i) when the joint dependence structure can be inferred from data, and (ii) when only marginal distributions and correlation matrices are known.

A growing body of work establishes that large language models are not mere statistical memorizers, but are capable of in-context learning: performing inference at test time using only examples provided in the prompt, without any parameter updates. Prior theoretical work has shown that this capability extends to supervised learning tasks such as linear regression. We prove that in-context learning extends further to data generation: frozen transformers can simulate iterative generative samplers from in-context samples. We first show that transformers can realize closed-form and smoothed closed-form diffusion samplers. The construction identifies a concrete generative role for softmax attention: it computes responsibility weights and weighted empirical averages, while feedforward layers implement Euler updates. To empirically relate these constructions to pretrained language models, we study semantic-topic sampling: prompts consisting of words drawn from a common semantic category, such as animals, foods, or cities. Across transformer layers, the normalized hidden states exhibit a two-stage geometry: they move toward a uniform spherical reference in intermediate layers and then return to structured, topic-dependent representations near the output. We further measure an interacting-particle energy on these hidden-state clouds and observe the same U-shape pattern. We then prove that transformers can approximate an energy-based sampler, constructing the same U-shape energy across the layers.

Bayesian palaeoclimate reconstruction from fossil pollen counts relies on a modern pollen-climate calibration data set to infer the pollen-climate relationships used to reconstruct past climates. While geographically large calibration data sets improve coverage of climate space and reduce unreliable extrapolation, they also introduce substantial heterogeneity, structural zeros, and complex pollen-climate relationships. We propose a Bayesian modular framework for palaeoclimate reconstruction from count-compositional pollen data that addresses these challenges, and provides coherent uncertainty quantification. The framework employs the zero-and-$N$-inflated multinomial logistic-normal distribution to describe the compositional pollen counts coupled with Bayesian additive regression tree priors to model the nonlinear effects and interactions among the climate covariates. Inference is formulated through a cut posterior distribution that modularises the analysis into forward and reconstruction modules. The forward module is fitted once using a large modern calibration data set and its posterior uncertainty is subsequently propagated to reconstruct climate variables from fossil pollen counts. For the reconstruction module, we develop and compare three inverse posterior sampling schemes. Simulation studies and empirical validation on the modern data set demonstrate that a combination of sampling importance resampling with a multiple imputation technique and a continuous uniform prior over the domain of the modern climate variables achieves the best predictive performance, with well-calibrated uncertainty quantification for the climate reconstruction. In our motivating case study, we further illustrate the proposed methodology by reconstructing a three-dimensional climate vector from fossil pollen records collected at Lago Grande di Monticchio in southern Italy.

Outlier detection is a fundamental challenge in data processing, with critical implications for robustness across statistical modeling, machine learning and exploratory data analysis. However, existing proposals rarely offer a universal, domain-agnostic definition of an outlier, often relying on heuristic trimming quotas that lack a statistical interpretation. To address this, we propose a nonparametric framework built on a pseudo-isolation outlier score. This score enables a formal, probabilistic definition of an anomaly tied to a false-alarm rate $\alpha$, which extends into a rigorous geometric classification of internal and external outliers. We show that this mechanism seamlessly embeds into any objective-based clustering framework to identify cluster-specific outliers. Here, we integrate it into $K$-means to create ODK-means. The inferential capabilities and topological properties of this framework are explored both theoretically through formal propositions and empirically through extensive simulations and methodological tutorials, highlighting the practical actionability and intuitive appeal of the proposed outlier detection logic.

Group testing pools specimens to cut the cost of screening, but pooling positive specimens with negative ones lowers assay sensitivity, so that the sensitivity of a pooled test depends on how many of its members are positive. Regression models for group testing accommodate this dilution effect through submodels that fix a one-parameter shape for that dependence. We propose BADGER (Bayesian Analysis of Dilution in Group tEsting Regression), a regression model incorporating a shape-free dilution curve. The pooled sensitivity is modeled as a nondecreasing function represented by nonnegative increments with a Dirichlet prior, so that the parametric submodels become special cases. By an appropriate augmentation, every full conditional of BADGER is in closed form and inference is carried out by an exact Gibbs sampler. The evidence for dilution is measured by a Bayes factor computed from the posterior draws. We validate the BADGER model through simulation studies across different dilution shapes, prevalences and pool sizes. We further illustrate the method on hepatitis E serology from a national health survey.

Slice sampling is a Markov chain Monte Carlo algorithm that draws its next state uniformly from a "slice"—a super-level set of the target density function—at each iteration, thereby providing automatic local adaptivity to the scale of the target. In practice the exact slice is not known, so general-purpose implementations use an approximate slice that is grown from a starting interval of length $w>0$, with a computational cost that depends on $w$. This work presents an analysis of the average per-iteration number of target density evaluations, as a function of $w$, of hybrid slice sampling with various slice-finding schemes for targets with contiguous slices. The paper uses the results of the analysis to develop automated, slice-adaptive tuning schemes along with suboptimality bounds and asymptotic convergence guarantees. Simulations demonstrate that the tuning schemes reliably yield near-optimal slice-adaptive tuning with essentially no dependence on the initial setting of $w$.

Traditional Bayesian optimal experimental design (OED) selects measurements that best inform a model's parameters. However, such measurements can be suboptimal for downstream predictions. Goal-oriented OED targets the prediction directly. However, the existing goal-oriented criteria value all reductions in predictive uncertainty equally, with no way to prioritize rare, high-consequence outcomes. In this article, we develop a risk-aware framework that composes risk at three levels, each generalizing an ingredient of classical $I$- and $G$-optimal design: a deviation measure of the posterior predictive uncertainty (generalizing the predictive variance), a risk measure across the prediction domain (interpolating $I$-optimal averaging and $G$-optimal worst-case selection), and a risk measure over datasets (generalizing the expectation). We generate each level from a regret function in the risk quadrangle, so that one triple specifies a practitioner's risk preference. We relax the design to continuous weights on the unit simplex and construct a nested-quadrature estimator that is differentiable in the design variable. This enables solving the optimal design problem with gradient-based methods, avoiding a combinatorial search over candidate designs. For a linear-Gaussian lognormal model and a nonlinear extension, we derive closed-form objectives. These give exact references against which we verify that the estimator converges. We demonstrate this framework for finding optimal sensor placements in an inverse problem governed by an advection-diffusion equation. We find that the risk-aware designs substantially outperform the expected-information-gain baseline, which is statistically indistinguishable from a random allocation.

Markov chain central limit theorems (CLTs) and their associated variances are very important for implementing Markov chain Monte Carlo algorithms among other applications. Häggström and Rosenthal (2007) presented various results regarding the equality of different formulae for this variance and also posed seven open problems. We resolve all seven in this paper. For stationary, ergodic and reversible chains, we prove that whenever the normalized partial sums satisfy a $\sqrt{n}$-CLT, the variance limit is finite if the function is square-integrable, otherwise undefined. Moreover, failure of the $\sqrt{n}$-CLT forces the normalized partial sums to be non-tight. We also show that Roberts' holding-probability condition precludes a CLT even without assuming reversibility or square-integrability. Finally, we develop a general principle that expresses Fourier coefficients as the autocovariances of an ergodic nonreversible Markov chain.

The concept of Markov chain Monte Carlo (MCMC) cycles, an analogy to cyclic processes in heat engines, is presented in order to examine Bayesian inference problems. In this effort, we develop adaptive ensemble schedulers that allow the tuning of external parameters of a Bayesian canonical ensemble during an MCMC run, realising the MCMC cycles in practice. We run these cycles on different statistical models. As a fundamental insight, we find (both theoretically and in practice) that such systems can produce a non-zero net work output if and only if the considered model is non-Gaussian. As such, they may serve as a measure of non-Gaussianity in Bayesian inference, which we test on an example from supernova cosmology.

Model and variable selection are important topics in Bayesian inference. In particular, the selection of different fixed, random effects and hyperparameters in hierarchical models can be difficult because of their complex structure. In this paper, we introduce the use of approximate Bayesian inference for model and variable selection for hierarchical Bayesian models. The approach is based on the application of Markov chain Monte Carlo methods on the model space. In particular, the Metropolis-Hastings algorithm is used to sample from the set of model indices. In this way, the marginal likelihood is used to compute the acceptance probability so that it is not required to estimate the parameter models directly. For this, the integrated nested Laplace approximation (INLA) is used because it provides accurate estimates of the marginal likelihood and it can also provide estimates of the posterior marginal of the model parameters. This method can be applied not only to variable selection but to a wide range of problems subject to model uncertainty. To illustrate the potential of this approach, examples on variable selection, changepoint models and log-Gaussian Cox processes are developed.

We propose a Bayesian generalized network autoregressive (BGNAR) model for multivariate time series whose component series are associated with the nodes of a known network. The proposed framework combines the parsimonious network structure of the generalized network autoregressive (GNAR) model with structured shrinkage and persistence priors adapted from Bayesian vector autoregressive (BVAR) modeling. We adapt Minnesota-type shrinkage to both own-lag and network-lag coefficients, with prior variances decreasing over temporal lags and, for network effects, neighborhood orders. A hierarchical prior on the own-lag coefficients allows information to be shared across nodes while retaining node-specific heterogeneity. We further adapt the sum-of-coefficients and dummy-initial-observation priors to the GNAR parameterization. Posterior inference is performed using a Gibbs sampler. Simulation studies show that BGNAR can use the same deliberately over-specified temporal and neighborhood structure across datasets without dataset-specific BIC order selection, while maintaining forecasting accuracy comparable to BIC-selected GNAR and outperforming the unrestricted BVAR benchmark in the settings considered. The structured prior regularizes weakly supported coefficients toward zero within this fixed model. Posterior distributions for the dynamic coefficients and posterior predictive distributions for future observations provide direct quantification of parameter and predictive uncertainty. An application to a wind-speed network demonstrates that BGNAR achieves point-forecast performance comparable to GNAR while additionally providing posterior inference on own-lag and network-lag effects and posterior predictive uncertainty.

We study trend estimation in state-space models in which the trend has a fractional stochastic difference of order $d>0$ and the observation errors form a short-range-dependent stationary process. Using finite-sequence fractional summation and differencing operators, we analyze the penalized least-squares estimator obtained by shrinking the fractional differences of the trend. We derive asymptotic mean squared error (MSE) approximations for all $d>0$ and identify a sharp phase transition at $d=1/2$. When $d>1/2$, the estimator is consistent and its optimally balanced MSE has order $n^{-(2d-1)/(2d)}$. At the boundary $d=1/2$, we obtain a refined finite-sample approximation and show that the MSE decreases at the slower order $\log\log n/\log n$. When $0<d<1/2$, the MSE converges to an explicit positive limit, so consistent recovery of the trend is impossible under the considered scaling. We also describe a practical criterion for choosing the penalty parameter and differencing order, and numerical experiments illustrate the MSE approximations and the behavior of the selection procedure.

Clinical prediction models are increasingly fitted by penalized logistic regression, because collinearity or many candidate predictors makes maximum likelihood unstable or impossible. Calibration is then almost always assessed by a grouped goodness-of-fit test such as the Hosmer-Lemeshow test. We show that this combination is invalid. Under ridge regression the grouped standardized residuals acquire a non-centrality induced by shrinkage, so the reference distribution used in practice is wrong, and at the penalty that most improves the fitted probabilities the test rejects correctly specified models between 92 and 100 per cent of the time. We derive the corrected law and define the shrinkage-corrected Hosmer-Lemeshow test, which subtracts an estimate of that non-centrality, restoring the maximum likelihood reference exactly to first order, and is made valid by prepivoting at a power cost we measure. We also give the attenuation law governing what any such test can detect once the linear predictor must be estimated. In glaucoma diagnosis by confocal laser tomography, where the maximum likelihood estimate does not exist, the corrected test finds the evidence for misfit weaker by more than three orders of magnitude: the fitted risks are too flat rather than mis-ordered, so the model needs recalibration rather than rebuilding.

Modern statistical models are growing increasingly complex in an effort to realistically capture system dynamics. Using standard simulation-based inference, these models may be computationally prohibitive, necessitating the use of model calibration methods. Bayesian score calibration is a computationally efficient framework for model calibration with strong theoretical guarantees. This framework learns an appropriate correction for an approximate model using a small number of simulations from the data-generating process. Currently, only a location-scale transformation has been explored, which may lack the flexibility to correct the complex error introduced by some approximate models. In this paper, we develop two flexible transformations for use in the Bayesian score calibration framework. The first is a polynomial extension, which can appropriately adjust approximate models with location-varying error. The second is a sequential application of Bayesian score calibration, which can accommodate approximate models with posteriors that have low support for the true parameter values. We also discuss an additional diagnostic for use with this framework. We demonstrate the increased flexibility these two approaches provide over Bayesian score calibration in two illustrative simulation studies.