Papers

Building on the James–Stein approach to leading eigenvector estimation (Goldberg and Kercheval, Proc. Natl. Acad. Sci. USA 120, e2207046120, 2023), we develop a data-adaptive augmented James–Stein shrinkage framework for estimating leading eigenvectors and eigenspaces under a generalized spiked population model, in the high-dimensional regime where the dimension $p$ and sample size $n$ grow proportionally. For each spiked eigenvector, we construct an augmented target subspace that combines auxiliary information, either from domain knowledge or prior information, with the remaining sample spiked eigenvectors. This augmentation allows information shared across the sample spiked components to be exploited while retaining a fully data-driven shrinkage rule. We show that the resulting eigenvector estimator strictly improves upon standard PCA whenever the target subspace contains nonvanishing information about the population eigenvector, while asymptotically reverting to PCA when the target is uninformative. The individual estimators further yield a nested sequence of estimators for all leading spiked eigenspaces, with analogous dominance properties. The proposed estimator also strictly improves upon the existing HDLSS-motivated shrinkage estimator in the proportional high-dimensional regime. Simulation studies demonstrate substantial finite-sample gains and robustness to misspecification of the number of spikes.

Non-negative matrix factorization (NMF) is a singular statistical model whose Bayesian asymptotics are governed by the real log canonical threshold (RLCT). We study the local geometry of the factorization map and derive an upper bound for the RLCT of NMF. Let $H$ be the model inner dimension and $H_0$ the non-negative rank of the true $M\times N$ matrix. Assuming that the true matrix admits a strictly positive factorization of inner dimension $H_0$ in the interior of the parameter domain, we prove, for smooth positive priors, that $\lambda\leq \{(H-H_0)\min(M,N)+H_0(M+N-H_0)\}/2$. This bound strictly improves the previous bound when $H_0\geq3$. The proof uses a local analytic normal form that separates independent linear coordinates from a residual matrix product. When $H=H_0$ also equals the ordinary rank of the true matrix, we obtain the exact value $\lambda=H_0(M+N-H_0)/2$. Under the standard assumptions of singular learning theory, these results bound the leading coefficients of the expected Bayesian generalization error and the Bayesian free energy.

Conformal prediction provides marginal coverage guarantees, yet practitioners may wonder if the observed coverage is truly abnormal or consistent with sampling variation. Inference for realized coverage has received comparatively little attention, especially for time series. We study split conformal prediction with adjacent calibration and test sets of temporally dependent data. Using the functional dependence measure, we derive non-asymptotic bounds on marginal coverage error without mixing assumptions, which can be difficult to verify and may fail even for simple short-memory models. We establish a Bahadur representation to derive, to our knowledge, the first central limit theorem for realized coverage of split conformal prediction under temporal dependence. A consistent block-based estimator of the standard error yields an asymptotically justified test. We further study long-memory time series, which remain understudied in conformal prediction. For Gaussian linear processes, we show how very strong temporal dependence can lead to a non-Gaussian limiting law of realized coverage and establish block-sampling inference with an estimated normalization. The resulting theory explains how temporal dependence changes coverage uncertainty.

We study conformalized quantile regression and conformalized median regression that fit a model on one block of a time series and calibrate the conformal interval on the adjacent block. The existing theory of conformal prediction for time series rests largely on mixing conditions, which are hard to verify from a time-series model and fail for many standard processes, including simple ones with short memory. We replace this theoretical toolbox with the functional dependence measure, which in principle accommodates long-memory observations. The accuracy of the conformal interval length for time series has been understudied. To the best of our knowledge, this paper is the first work that establishes non-asymptotic coverage guarantees and accuracy of interval length simultaneously for split conformal regression on time series. Furthermore, for Gaussian linear processes with long memory, where both the estimation of the center and its calibration converge slowly, we establish a sharper rate for the length error. We show that the calibrated length converges faster than the estimated center itself, and provide a matching lower bound for the usual centers when the calibration block is sufficiently large relative to the training block. To our knowledge, this is the first theoretical analysis of conformal interval length dedicated to long memory.

Weighted split conformal prediction reweights calibration scores by the likelihood ratio between the test and training covariate distributions and guarantees marginal coverage under covariate shift. We study its coverage conditional on the calibration data. An elementary argument, based on a single concentration inequality at a fixed population quantile, gives explicit training-conditional bounds without unspecified constants, and shows that the relevant scale is not the supremum of the likelihood ratio but a variance proxy built from the chi-squared divergence of the shift and from the average of the ratio over the part of the test population, of probability equal to the miscoverage level, where it is largest. A two-point lower bound shows that the root-m rate and the chi-squared contribution are intrinsic to the shift. Run at an explicitly inflated level, the weighted quantile becomes a deterministic PAC prediction set. We compare it with randomized rejection sampling and with importance-weighted learn-then-test and, through a certified choice of a clipping level for the likelihood ratio, map the regime in which each gives the narrower valid set. The analysis extends to estimated likelihood ratios and to tail functionals estimated from an unlabeled source sample, which yields a fully finite-sample certificate.

Locally stationary Markov chains have attracted a lot of attention in statistics in the past three decades. In this paper we prove a variety of limit theorems for partial sums generated by such chains. We first provide explicit formulas for the asymptotic mean and variance (and also higher moments), and obtain optimal convergence rates towards them. We then prove a Berry–Esseen theorem, a local central limit theorem, Edgeworth expansions, and large and moderate deviations principles. Our approach involves a parametric Perron–Frobenius theorem, which is proved using the theory of complex (Hilbert) projective metrics developed in \cite{Rugh,Dubois}, together with local approximation arguments and ideas in \cite{dolgopyat2023berry}.

Recovering probability measures from their moments has numerous applications, esp. in connection with the method of moments in statistics and optimization. In the setting where measure need not be finitely atomic, but its support is known to be compact and semialgebraic with codimension at least one, the problem is still open. We combine moment-matrix kernel information with the Christoffel–Darboux kernel to provide a discrete approximation of the support. To validate the proposed approach, we test our algorithm on analytically computed moments and pseudo-moments arising from polynomial optimization problems without unique global minimizers. This complements well-known recent work on recovery of measures with algebraic support, where the kernel of a moment matrix can reveal polynomials vanishing on the support, and on recovery of sufficiently regular full-dimensional supports, where estimators constructed by thresholding the Christoffel–Darboux kernel are known to converge asymptotically to the support.

We propose a new $L_2$-type test for white noise which allows the dimension $p$ of the time series to either (i) be a fixed constant, or (ii) diverge with the sample size $n$. The proposed test statistic exhibits an interesting phase transition, following two different regimes of behavior: $p$ is fixed, and $p\rightarrow\infty$. Because identification of the operable regime is difficult, if not impossible in practice, we devise a novel adaptive bootstrap method to construct unified testing procedure across different phases. Numerical experiments confirm the good finite sample performance of the proposed adaptive $L_2$-type test in comparison to the existing methods in the literature. The proposed testing procedure has been implemented in R package HDTSA.

A special semimartingale admits a unique decomposition $X=X_0+M+A$ into a local martingale $M$ and a predictable finite-variation part $A$. We consider the identification of $A$ when $X$ is observed only through repeated cross-sections. The estimand is then the projection of the sampled predictable compensator onto the observable feature filtration, namely the current state together with whatever randomness is shared across the population, so that at a fixed diffusion coefficient the marginal flow identifies the drift only up to a Markovian projection. If the drift is an affine functional of an observed lag window, the joint problem is a convex quadratic programme whose solution is the pseudo-panel regression of econometrics. Our principal concern is the case, which we believe not to have been treated before, in which the drift is the output of a hidden linear dynamical system whose dynamics are themselves to be identified from the marginals. The joint problem is then a bilinear quadratically constrained programme, which we solve to certified global optimality by spatial branch and bound; with unpenalised state disturbances and a drift basis growing with the grid it is NP-hard already in latent dimension one, by reduction from $\ell^1$ rank-one matrix approximation, whereas the complexity of the deterministic system at fixed latent dimension remains open. A block-coordinate decomposition offers a cheaper alternative. For the estimator itself, we obtain rates at a fixed mesh, separated into Monte-Carlo, estimation and grid contributions.

Many generative systems return multiple candidate responses and are evaluated according to the best one. Recent work shows that, when demonstrations are optimal, pass@$k$ can reduce the sample complexity of learning from demonstrations by a logarithmic factor in $k$. We ask what happens when the demonstrator is not assumed to be optimal. We find that multiple responses provide a qualitatively stronger benefit in this setting. In a finite reward-class model with no reward feedback, moving from pass@$1$ to any pass@$k$ with $k\ge2$ changes the worst-case dependence on target accuracy from $1/\varepsilon^2$ to $1/\varepsilon$, uniformly over demonstrator quality. Under standard evaluation, where an unknown reward is fixed before training, increasing $k$ provides an additional and distinct benefit: the optimal dependence on a reward class of size $N$ improves from $\log N$ to $\log N/\log k$. We further show that these two effects can be separated. Under robust evaluation, where one learned policy must compete with the demonstrator simultaneously for every reward in the class, the fast $1/\varepsilon$ dependence persists, while the $1/\log k$ improvement can disappear. We establish matching upper and lower bounds in the corresponding regimes and give a greedy multiplicative-weights learner achieving the upper bounds without any assumption on demonstrator quality.

We study robust reconstruction for the \(q\)-state Potts broadcast process on a Galton-Watson tree in the growing-alphabet regime \(q\to\infty\), with boundary depth growing sufficiently quickly in relation to \(q\). We show that, in the regime \(d\lambda >1\), where \(d\) is the expected number of offspring and \(\lambda\) is the nontrivial eigenvalue of the Potts transition matrix, the root posterior is asymptotically unchanged by a broad class of noise channels applied independently to the boundary labels. This remains true even when the probability of retaining the true label at an individual boundary vertex tends to zero as \(q\to\infty\). More precisely, under suitable moment and noise assumptions, the noiseless and noisy root posteriors converge to one another in expected total variation, and hence have the same asymptotic Bayes-optimal reconstruction accuracy.

Motivated by blind denoising in diffusion models, we study estimation of an unknown Gaussian noise level from a single high-dimensional observation, assuming the signal law P is known. We characterize the minimax mean-squared error under two structural assumptions on P. For signals with covering complexity k, the minimax rate is $\widetilde{\Theta}_\Lambda(\min\{\Delta_\Lambda^2,d^{-1}+k^2d^{-2}\})$, and the constrained MLE attains it up to logarithmic factors. For $\alpha$-strongly log-concave signals, the rate is $\Theta_\Lambda(\min\{\Delta_\Lambda^2,(1+\alpha^{-1})^2d^{-1}\})$, attained up to constants by a $P$-centered norm estimator. These results show that the structure of the signal law determines both the difficulty of blind noise estimation and the appropriate estimator.

We study entity ranking and inference under the Bradley–Terry–Luce (BTL) model using pairwise comparisons collected over sparse comparison graphs. In practice, collecting high-quality human judgments can be expensive and time-consuming, motivating the use of synthetic data augmentation in model training. We analyze an iterative synthetic augmentation workflow in which synthetic comparisons generated from fitted models are successively added to the original dataset. In this process, the proportion of real data may vanish as the number of iterations grows. However, emerging literature has shown that recursive training on synthetic data can lead to model collapse, raising concerns about the statistical reliability of such augmentation procedures. To this end, we systematically analyze the resulting MLE in finite-sample, high-dimensional regimes. For the resulting iterative maximum likelihood estimator (MLE), we derive its optimal finite-sample $\ell_2$ and $\ell_{\infty}$ statistical rates and establish its asymptotic normality under natural identifiability conditions. We further characterize regimes in which model collapse is avoided despite the diminishing fraction of real data. We validate our theoretical findings through large-scale numerical experiments and an application to the Arena Human Preference 140k dataset.

We propose sliced Orlicz-Wasserstein (SOW) distance which is a generalization of sliced Wasserstein (SW) distance. SOW replaces the $L^p$ norm in SW with a Luxemburg norm cost induced by an Orlicz function $\phi$. First, we prove that SOW distance is a metric on the space of measures with finite Orlicz norm, and show that it recovers the SW distance when the Orlicz function is $\phi(x)=x^p$. Next, we derive the topological properties of the SOW distance. In particular, we show that convergence under SOW implies weak convergence, and the converse is true under the compact support condition. We then present the theoretical results for estimating the SOW distance. We derive sample complexity for both the distance itself and the powered functional of the distance, and prove their minimax optimality. In addition, we discuss the computational algorithm for approximating the SOW distance by Monte-Carlo estimation and bisection search, as well as the associated approximation error and computational complexity analysis. Our experimental results reveal the superior computational efficiency of SOW compared with Orlicz-Wasserstein (OW) distance. Also, in the experiments, we demonstrate the favorable flexibility of SOW distance over SW in detecting differences between distributions by comparing their performance in two-sample tests and evaluating generative models on image datasets.

Domain Adaptation (DA) often lacks worst-case guarantees, while Distributionally Robust Optimization (DRO) based only on source data centers its ambiguity set at the source law and ignores available target structure. To bridge this gap, we introduce a doubly-anchored DRO framework whose ambiguity set is the intersection of $\phi$-divergence balls centered at the source law and a source-completed target reference law, the latter pairing the target covariate law with the source conditional law. We derive dual-induced adversarial bridge geometries for symmetric and asymmetric divergence pairings, notably introducing a Kullback–Leibler/squared-Hellinger (KL/HD) bridge. This asymmetric formulation yields a Lambert-$W$ geometry in which source-side exponential risk tilting and target-side Hellinger stabilization enter through distinct terms, attenuating, but not bounding, the effect of large likelihood ratios. Furthermore, without imposing covariate shift, we establish finite-sample generalization bounds for the minimizer of a structural, loss-agnostic density-bridge risk under bounded-overlap conditions; these bounds do not apply directly to the loss-aware DRO min–max estimator. We translate our framework into a bridge-weighted Nadaraya–Watson estimator, proving uniform consistency for the source regression function and pointwise asymptotic normality, with target recovery when the source and target regression functions coincide, as under covariate shift. Finally, an empirical evaluation on a domain-shifted Fashion-MNIST dataset illustrates the finite-sample stability of the asymmetric KL/HD bridge under severe synthetic target-covariate corruption.

For a class of regular estimators, Hájek, in his celebrated “Convolution Theorem,” showed that the asymptotic distribution of a regular estimator is the convolution of the distribution of an efficient estimator and some residual distribution. This result constitutes the foundation of the concept of asymptotic efficiency of regular estimators. In this paper, we provide a second-order version of that classical result. Introducing a class of second-order regular estimators with a valid Edgeworth expansion, we derive their asymptotic distribution under contiguous alternatives and show that it is the convolution of the second-order efficient distribution and some second-order residual distribution. This constitutes a second-order extension of Hájek's convolution theorem. Based on this, we introduce a concept of second-order robustness for second-order regular estimators. For a class of general Bayes estimators and minimum contrast estimators in time series models, this second-order robustness is used (i) in the characterization of second-order robust priors, (ii) in a comparison between the second-order robustness of maximum likelihood and Whittle estimators.

E-processes provide a flexible framework for anytime-valid inference, with the conventional rejection boundary $1/\alpha$ typically justified by Ville's inequality. Such a boundary is universal but can be conservative, as it does not exploit additional information about the null distribution. We propose a reference-null calibration framework that uses independent null samples to construct sharper rejection thresholds while preserving type-I error control. To study the statistical gain from sharper thresholds, we specialize the framework to conformal martingales. Building on histogram-based betting, we incorporate Krichevsky–Trofimov smoothing and develop a restart-mixture construction for distribution shift detection. We establish quantitative results for detection power and detection delay and characterize how the reduction in the rejection boundary translates into power and delay gains relative to the conventional Ville boundary. Numerical experiments corroborate the theoretical findings across a range of distributional changes. Finally, we apply the proposed calibration strategy to an existing e-process for online LLM watermark detection, showing that the method can improve sequential detection performance without modifying the underlying e-process. These results demonstrate that reference-null calibration provides a general and modular way to enhance the efficiency of e-process-based sequential inference.

We introduce curve band depth (CBD), a band-based data depth for samples of unparameterized planar curves. CBD is motivated by band depth and modified band depth for functional data, but targets trajectory data. Unlike the halfspace-based curve depth of \citet{de2021depth} and the curve stabbing depth of \citet{durocher2023csd}, CBD is defined through a geometric band region generated by two curves, and measures the arc-length proportion of a target curve lying inside such bands. We develop a CBD family consisting of an integral version (int-CBD), an infimal version (inf-CBD), and a fast-walk variant (FW-CBD). The fast-walk band is a narrower band construction contained in the global convex-combination band. We establish boundedness, vanishing at infinity, and similarity invariance for these constructions, together with a Borel-measurability result for the induced depth maps under a mild measurability assumption. A length-penalized variant is proposed for samples with heterogeneous curve lengths. We implement the methods via arc-length sampling and polygonal approximations, and evaluate them through classification of overlapping handwriting data and MNIST-derived digit curves, online-signature screening on MOBISIG, and an exploratory clustering task based on decomposed band contributions.

We study smoothed online prediction with bounded adversarial responses. This widely studied framework bridges i.i.d. sampling and adversarial covariate selection through a smoothness parameter $\textsf{C}_{\textsf{cov}}$, which bounds conditional covariate densities relative to a fixed, unknown base measure. We propose Hedge-Cover, an information-theoretic algorithm that achieves sublinear regret $\widetilde{O}(\sqrt{\text{Pdim}(\mathcal{F}) \textsf{C}_{\textsf{cov}} T})$ for function classes with bounded pseudo-dimension. The algorithm aggregates a carefully constructed family of experts using Hedge, with a prior that links regret to the number of disagreements between a consistent selector and a target function. We bound this number by exploiting covariate smoothness. This answers an open problem posed in \cite{blanchard2025agnostic} on the minimax optimal adaptive regret of the smoothed online regression problem. We establish a matching lower bound for the class of linear predictors. The main intricacy of the lower bound lies in explicitly constructing a challenging sequential covariate distribution supported on mutually orthogonal hyperplanes. This construction may be of independent technical interest. Finally, we revisit the well-specified setting and quantify the effect of response noise. For conditionally $\nu^2$-subGaussian responses, we extend the existing lower bound under realizable responses by showing that the minimax expected regret is $\Omega((1\vee \nu)\sqrt{(\textsf{C}_{\textsf{cov}}-1)dT})$ for a function class of VC dimension $d$. A corresponding upper bound for ERM matches this dependence on $\nu$.

Spectral risk measures, including conditional value-at-risk (CVaR), generate a family of convex criteria for portfolio estimation. We study how the criterion should be chosen when different criteria share the same population minimizer, and whether the efficient criterion can itself be learned from data. We develop a general asymptotic theory for empirical spectral-risk minimization under estimated linear constraints. Under normal scale-mixture elliptical returns, all spectral risk measures identify the same population efficient portfolio, but their empirical minimizers have different sampling distributions. Their asymptotic covariance decomposes into a common component and a positive-semidefinite component scaled by a functional of the spectral measure, reducing efficiency to an optimization over probability measures. We characterize the efficiency-optimal spectral measure and show that single-level CVaR is generally inefficient. We then construct a fully data-adaptive estimator that learns the radial distribution and optimal spectral measure from the same observations used for portfolio estimation, yet has the same first-order distribution as the infeasible oracle. Simulations and an empirical application illustrate the method.

Instead of having a single "yes" or "no" result from a test of the global null hypothesis that a function is increasing, we propose a multiple testing procedure of the function's increasingness at several points. If the global null is rejected, then multiple testing provides more information about why. If the global null is not rejected, then multiple testing can provide stronger evidence in favor of increasingness, by rejecting null hypotheses that the function is decreasing. Our approach uses high-level assumptions that apply to a broad class of causal and descriptive statistical models. By inverting the proposed multiple testing procedure that controls the familywise error rate, we also generate "inner" and "outer" confidence sets for the set of points at which the function is increasing. With high asymptotic probability, the inner confidence set is contained within the true set, whereas the outer confidence set contains the true set. We also improve power with stepdown and two-stage procedures. Simulation and empirical examples illustrate the new methodology, and all code is provided.

Generative AI systems are increasingly used, but aligning their outputs with user requirements poses a continuing challenge. Here, we aim to ensure that the distribution of an attribute of an AI-generated output aligns with a user-specified target. This is motivated by examples such as fairness, where we want to ensure that a protected attribute (e.g., gender, race, or age categories) follows a desired distribution, and synthetic data generation, where we want the generated data to be representative of a target distribution. We study the practically important black-box access setting, where a user can repeatedly query a generative AI model. The goal is to return $m\ge 1$ outputs whose joint attribute distribution is as close as possible to this target. For both exact and approximate alignment, we develop algorithms that minimize the expected number of queries to the generator, and we further demonstrate their optimality as the number of requested outputs $m \rightarrow \infty$. Experiments on text-to-image generation and geocoded persona generation tasks show that our post-processing algorithms improve statistical attribute alignment, complementing prompting-based interventions.

Bi-invariant geodesic regression generalizes linear regression to Lie groups. Its main feature is that it respects the symmetries of the group so that the resulting estimator is independent of arbitrary choices such as a reference frame. However, the local existence and uniqueness of the underlying estimator have not been shown until now. Furthermore, the convergence properties of the proposed algorithm for computing the estimator are not known. In this work, we investigate these questions. We prove that, locally, a unique estimator exists and give explicit bounds on the size of this neighborhood. We also show that the proposed iterative algorithm converges linearly to this estimator.

Transformers have become a central architecture for in-context learning (ICL), particularly through their state-of-the-art performance in large language models. This success motivates understanding how transformers exploit task-relevant structure in geometrically heterogeneous data. However, existing nonparametric ICL theory has largely focused on Euclidean domains or single-manifold models. To address this gap, we study the prediction problem under unknown local geometry, modeled by sample size-dependent mixtures of manifolds with heterogeneous dimensions, smoothness, and sampling masses. Under local separation and small-perturbation conditions, we establish a minimax lower bound capturing the aggregate difficulty of the components and construct an oracle tangent local-polynomial estimator with a matching upper bound. This estimator is connected to a structure-informed, two-stage softmax transformer with a geometric preconditioner and chartwise reduced local-polynomial solvers. The transformer achieves negligible approximation error relative to the minimax rate with logarithmic depth and polynomial size. Finally, we derive an in-context generalization bound for near empirical risk minimizers over this class. Together, these results identify conditions under which the resulting predictor exploits local geometry and attains the aggregate minimax rate.

This article is concerned with the asymptotic shape of quantile surfaces, defined as the set of quantiles at a given level $\alpha$ generated by a controlled one-dimensional distribution. Specifically, when the distribution arises as a linear combination of log-normal random variables and the control is a vector of positive coefficients, we prove that quantile surfaces are globally concave in the left tail ($\alpha\to0$) and globally convex in the right tail ($\alpha\to1$). Moreover, these surfaces exhibit asymptotic separation of scale and shape.

To design systems that are protected against events much rarer than the observational record, extreme-value methods are needed to extrapolate distribution tails. Tail-index estimators such as the Hill estimator are central to this extrapolation, but overestimating the tail exponent can lead to substantial underestimation of rare-event probabilities. Motivated by this, we derive large- and moderate-deviation asymptotics for the Hill estimator and use them to construct estimators whose probability of exceeding the true tail index decays at a controlled exponential rate (the decay rate). In the large-deviations regime, we show that a simple rescaled version of the Hill estimator achieves an optimal balance between bias and decay rate among scale-invariant estimators based on the same top $k$ order statistics. Under a second-order condition, we quantify the effect of the Hill bias, analyze a bias-corrected estimator, and identify sufficient conditions for moderate deviations in the boundary case where the second-order parameter $\rho$ equals zero.

Combining $p$-values is a fundamental procedure in global hypothesis testing. In modern high-dimensional settings, however, component $p$-values often exhibit complex dependence and rely on asymptotic approximations rather than exact finite-sample uniform distributions. This paper establishes a unified asymptotic theory for weighted transformation statistics that decouples marginal finite-sample approximation error from the joint dependence structure. We also provide sufficient conditions based on conditional probability bounds to verify the joint-tail conditions. Utilizing this framework, we derive explicit dimension-growth and correlation rates for test statistics operating under asymptotic Gaussian and chi-square calibrations. For non-exact finite-sample statistics, we analyze standardized weighted sums, demonstrating how Cramér moderate deviations control relative tail error. Analytical examples demonstrate why both marginal and joint conditions are mathematically indispensable for valid global inference under dependence, and numerical experiments confirm that our asymptotic framework maintains accurate finite-sample size control at extreme significance levels.

Martingale posteriors and related predictive resampling methods replace the likelihood–prior pair used within Bayesian inference with a predictive model for future observations. These methods are simple to implement and increasingly popular due to their computational efficiency, but little is known about their ability to accurately quantify uncertainty. In this work, we study the concentration and calibration properties of the martingale posterior for general functionals, and show that credible sets can systematically undercover if the predictive algorithms are not carefully tuned. We propose a simple remedy: the bagged martingale posterior. Rather than starting every predictive path from the observed sample, we start the paths from random bootstrap resamples of the data and then amalgamate the resulting draws. Critically, this approach incurs no additional {simulation} cost compared to standard predictive resampling, and delivers conservatively calibrated credible sets. This scalability enables application to a range of challenging examples, including sparse high-dimensional regression and nonparametric conditional quantile models.

Many real-world high-dimensional time series exhibit long-memory, but Gaussian graphical model testing in this regime remains understudied. We develop a direct, data-adaptive test statistic for assessing conditional independence in the graph structure of stationary Gaussian time series. We establish a finite-sample, Berry–Esseen type Gaussian approximation bound for the statistic, which applies to both short-memory and long-memory time series. The testing procedure is fully data-adaptive using block bootstrap method, on which we provide a finite-sample validity result including in the ultra-high-dimensional scenario, and can be extended to comparing graphical structures in two-sample tests. We also develop a consistency-empowered correction to the statistic and show that such tests attain asymptotic consistency in both size and power. Our proposed method is applied to a real-world fMRI data to understand functional connectivities within brain in different periods.

The purpose of this article is to explicitly describe how the integration technique known as the change of variables is used in graduate-level Probability and Statistics courses. Throughout the article, we employ this integration technique seven times across various examples and results. We demonstrate that four applications of this technique are required to derive the probability density function (PDF) of the Student's t-distribution.