Papers

This paper studies graph matching under the correlated $\text{Erd\H{o}s-R\'{e}nyi}$ (ER) graph pair model. This model first samples an $\mathrm{ER}(n,\frac{\lambda}{ns})$ base graph, whose edges are then independently subsampled twice with probability $s$ to produce two correlated $\mathrm{ER}(n,\frac{\lambda}{n})$ graphs. We propose a graph matching algorithm that has $n^{2+o(1)}$ time complexity and achieves almost exact recovery with high probability under the assumptions $\lambda=(\log n)^{\alpha+o(1)}$ for some $\alpha\in(0,1)$ and $s\in(\sqrt{C_{\mathrm{Otter}}},1]$, where $C_{\mathrm{Otter}}\approx 0.338$ is Otter's tree-counting constant. This is the first algorithm with almost quadratic time complexity in this regime of $\lambda$, while the best known result in this regime is the chandelier-counting algorithm with time complexity $O(n^{c(s)})$, where $c(s)\rightarrow \infty$ as $s$ approaches $\sqrt{C_\mathrm{Otter}}$ from above. The proposed algorithm is based on local tree correlation tests. It uses a rank-based algorithm to match the vertex pairs instead of threshold-based rules in the literature. This avoids the need of computing an explicit threshold, which is computationally difficult to obtain. To prove the almost exact recovery result, we establish a new analysis of tree correlation tests in the diverging-degree regime, where both the mean degree and the tree depth grow with $n$. Based on this new result, we establish the existence of a threshold for a threshold-based graph matching algorithm via local tree correlation tests. Finally, we couple the performance of the rank-based algorithm with the threshold-based algorithm to show almost exact recovery.

We study the problem of recovering the correspondence between a collection of $n$ points in $\mathbb{R}^d$ and a noisy, permuted version of those points. In the high-dimensional regime $d=\omega(\log n)$, under a Gaussian model with noise variance $\sigma^2=d/(b\log n)$, prior work identifies $b=2$ as the threshold for almost exact recovery. We prove that this threshold is all-or-nothing: for every fixed $b<2$, no estimator recovers a positive fraction of the matching, and even estimating the matched point cloud in Euclidean distance is asymptotically no better than ignoring the correspondence. On the other hand, we consider a multi-view generalization of the problem where $K$ noisy, independently permuted copies of the same latent point cloud are observed. Here we show that a simple polynomial-time procedure recovers all relative matchings up to $o(n)$ errors whenever $b>K/(K-1)$. Thus multiple views can break the impossibility barrier $b=2$ for the original matching problem: in particular, for $3/2 < b < 2$, the two-view model has no nontrivial recovery, but a third view makes all latent correspondences efficiently recoverable.

We explain how important classes of Bayesian semiparametric regression fitting and inference procedures can be sped up, significantly, via the use of orthogonalized design matrices. Typically, design matrices in semiparametric regression contain predictor observations and basis functions of such data. In Bayesian semiparametric regression, loop-type approaches such as Gibbs sampling and coordinate ascent variational inference typically are required. We show that pre-loop reformulation of Bayesian semiparametric regression models involving orthogonalized design matrices lead to two orders of magnitude, with respect to column dimension, computational reduction. Our computer experiments reveal that this simple paradigm results in approximately 5- to 60-fold speed-ups.

In exploratory factor analysis (EFA), one aims to describe latent variables by constructing a factor model based on the relationships among manifest variables. For a model to be useful, it is not enough that it is grounded on data; it must also be meaningful. Hence, in practice, one attempts to interpret different factor models to identify a meaningful, coherent, and theoretically defensible latent structure. Doing so, however, is not straightforward, as it is subjective and requires tracking extensive information. Thus, we introduce FactorFlow, a system designed to help researchers perform EFA more effectively. With an interactive dashboard that supports comprehensively visualizing up to two models simultaneously and large language model integration that enables the generation of automated model interpretations written in natural language, FactorFlow substantially aids the crucial step of model interpretation, all the while supporting the end-to-end workflow. Indeed, our usability survey evidences the effectiveness of FactorFlow.

The state of a dynamic system evolves over time, switching among several latent modes that govern its observable behavior. Filtering methods infer the latent state from observations. Classical filtering approaches, including Kalman filters, typically rely on simple observation models, such as linear-Gaussian models, that are incapable of characterizing the increasingly nonlinear and heterogeneous patterns in high-dimensional sensor signals. To tackle the challenge, we propose Generative Bayesian Filtering (GBF), a filtering framework that replaces restrictive observation models with pretrained conditional generative models parametrized by conditional variational autoencoders (CVAE). For online inference, GBF performs a Bayesian prediction-update recursion in which the measurement update is formulated as a posterior sampling problem that combines the dynamical prior with the CVAE-induced likelihood. The resulting filtering problem is then transformed into a score-based sampling problem, which naturally inherits the flexibility from generative models and the uncertainty quantification capabilities from ensembling. Experiments on synthetic datasets and real-world applications involving manufacturing system monitoring and arrhythmia diagnosis demonstrate that GBF improves state estimation accuracy and robustness relative to baseline approaches.

Spectral embedding methods are widely used for dimensionality reduction and clustering of high-dimensional datasets with intrinsic low-dimensional structures. Although many datasets of practical interest exhibit invariance under symmetries such as rotations, standard spectral embedding methods do not account for this, treating symmetry-related data points as unrelated. Our approach to this problem is to incorporate the symmetries directly into the affinity kernels used for spectral embedding. We analyze the case of a Riemannian data manifold $M$ with symmetries given by a compact Lie group~$G$ and prove that, under suitable conditions, graph Laplacians constructed from three types of invariant kernels converge pointwise to explicit second-order differential operators on the quotient space $M/G$. Our analysis implies improved convergence rates, as the effective dimension drops according to the dimension of the group. We validate our approach on datasets with $\mathrm{SO}(2)$ or $\mathrm{SO}(3)$ symmetry, and show that $G$-invariant spectral embedding recovers the intrinsic geometry of the data, in contrast to standard spectral embedding, which fails to do so even in the limit of infinite data.

We study the active learning problem of fixed-confidence top-$k$ identification from noisy pairwise comparisons. In this problem, an algorithm sequentially chooses pairs of items to compare, observes the outcomes, and stops when it can return the set of top-$k$ items with error probability at most $\delta$. The objective is to design such a $\delta$-correct procedure that minimizes the expected number of comparisons (the sample complexity). This problem falls within the broader literature on fixed-confidence pure exploration in bandit models, where a common target is asymptotic optimality: the algorithm's expected sample complexity matches the information theoretic lower bound as $\delta \to 0$. Asymptotically optimal procedures have been developed for a range of fixed-confidence pure-exploration problems, however to the best of our knowledge, for top-$1$, or more generally top-$k$ identification from pairwise comparisons under latent utility models an asymptotically optimal algorithm has not been established. In this setting, we develop such an algorithm. We characterize the structure of the lower bound and formulate it as a saddle-point problem. This structure enables a computationally efficient primal-dual procedure that learns the asymptotically optimal comparison allocation online. We then construct an adaptive comparison-allocation algorithm that tracks the allocation learned by the primal-dual procedure and prove it is asymptotically optimal.

The stochastic linear bandit, where actions are represented as vectors and rewards are linear, is a central paradigm for sequential decision making. We study a partially observed variant of this problem in which the learning agent only sees a random subset of coordinates for each action. Such partial observability arises naturally in settings like recommendation and healthcare, where full action descriptions can be expensive or even impossible to obtain. In general, this makes sublinear regret information-theoretically impossible. However, we show that this barrier can be overcome when the action vectors have low intrinsic dimension. We propose an algorithm, TOFU-POV, that estimates the latent action subspace using the masked actions, imputes current actions using an epoch-wise frozen representation, and runs OFUL in the resulting low-dimensional coordinates. Our theory shows that TOFU-POV enjoys a $\sqrt{T}$ regret that scales with the intrinsic action subspace dimension as opposed to the ambient dimension and quantifies the interaction between these quantities and the missingness, decision set size, and subspace conditioning. We also devise a rank-adaptive algorithm that does not require the knowledge of the intrinsic dimension. We complement these guarantees with a lower bound based on a novel product construction that separates usual reward-learning uncertainty from a missingness-dependent cost intrinsic to partial observation. Synthetic and real data experiments support our theory and show that TOFU-POV can substantially improve upon natural baselines in this challenging problem.

Large language models increasingly provide labels, evaluations, and feedback for tasks specified in natural language. When a specification admits multiple readings but the supervision channel does not reveal which is operative, additional labels reduce sampling error without resolving the resulting identification problem. We introduce Natural Language PAC (NL-PAC), a framework that uses a fixed model's thresholded decoding law to define admissible labels and candidate targets. The probability that multiple labels are admissible equals the diameter of the pointwise-admissible target class, and under target-blind supervision every learner incurs worst-case risk of at least half this diameter, at every sample size; the exact randomized minimax risk over this class is attained by a data-independent strategy. Finite-sample confidence bounds make these quantities certifiable from held-out unlabeled inputs. In a frozen Qwen~2.5–3B audit, one prespecified prompt yields a positive model-relative certificate, whereas a paraphrase and exact-rule controls yield zero. A held-out bridge audit finds that supplied candidate reading clauses fail the admissibility condition needed to transfer the certificate to coherent readings. The guarantee is specific to the audited model, prompt, threshold, and input distribution; extending it to human interpretations requires external validation.

We study nonasymptotic convergence of primal-dual methods for a class of nonconvex constrained optimization problems under local multiplier regularity. Both the objective and functional inequality constraints are compositions of convex Lipschitz outer functions with smooth nonlinear mappings. Unlike existing analyses requiring regularity globally or along the entire iterate sequence, we impose multiplier regularity only near the feasible set. Using an auxiliary minimax problem with a compact dual set, we derive Lyapunov estimates that identify nearly feasible iterates under feasible initialization and a sufficiently large penalty, without assuming multiplier boundedness. Local regularity then bounds the associated multipliers and makes dual truncation inactive, yielding approximate KKT certificates for the original problem. The smoothed prox-linear augmented Lagrangian method with dual regularization yields an $O(K^{-1/3})$ KKT residual rate. Without regularization, piecewise linear outer functions and additional local structural assumptions yield a local dual error bound and the sharper $O(K^{-1/2})$ rate.

From medicine to marketing to social sciences, the promise of tailoring interventions to individuals is undeniable. However, practical applications force weighing personalization's potential benefits with its possible increased cost and fragility. We introduce a statistical hypothesis test that evaluates, given historical data, evidence that a personalized intervention policy's performance will surpass deploying the best single intervention. The test maintains strict type-I error control while achieving asymptotic normality with the minimal possible variance under specified conditions. Results on diverse datasets from job training, depression treatment, education and recommendation systems demonstrate the test's versatility and its superior performance over alternatives. This test can support decision-makers throughout the intervention sciences by providing a simple and powerful quantification of the potential benefits of personalization.

Predictions with covariates in the tails of the covariate distribution often lack accuracy and are highly uncertain but may be of critical importance in many applications. New covariates may even lie beyond the range of training data. The problem can be particularly critical in environmental contexts, for example with climate-change scenarios, where new covariates represent future, more extreme climate conditions. We here propose novel methods to improve generalized additive models (GAMs) with focus on extreme covariates and covariate extrapolation. Adopting a random-design setting, we continuously integrate GAMs for the bulk of covariate distributions with models motivated by multivariate extreme value theory for high covariate values. We consider continuous responses but develop also a new approach for binary responses by assuming a continuous latent response variable. Our framework imposes a specific structure for large values of the covariates motivated by extreme value theory: by combining a transformation to a specific marginal scale with an appropriate link function, the continuous response variable depends linearly on the covariates. In an application to occurrences and sizes of large wildfires in Europe for the period 2008-2023, we explore how the new method can improve predictions, especially during extreme conditions, using environmental and meteorological covariates.

The functional organization of the brain relies on coordinated activity across spatially distributed regions, making the analysis of inter-regional dependencies fundamental. Existing connectivity measures address this predominantly through phase synchronization, which is vulnerable to volume conduction artifacts and discards amplitude-domain coupling. This study introduces the Spatial Neighboring Scattering Transform, which extends the wavelet scattering transform to the multichannel setting, yielding two descriptors that jointly capture amplitude-envelope coupling between channels and its modulation across frequency scales. SNST was evaluated on the BCI Competition IV-2a motor imagery dataset using a bias-corrected, false-discovery-rate-controlled statistical pipeline, with the validation criterion defined as spatial consistency of significant coupling across subjects. The first-order descriptor identified statistically significant amplitude coupling within a central-parietal electrode neighborhood, reproduced consistently across all subjects and both imagery conditions. The second-order descriptor revealed that this coupling is periodically gated by slow rhythms, indicating a cross-frequency amplitude-modulation structure absent from single-frequency connectivity measures. Phase lag index and weighted phase lag index, computed under an identical correction procedure and verified robust to volume conduction, identified negligible significant coupling with zero overlap with SNST findings, demonstrating that amplitude envelope coupling constitutes a largely distinct connectivity signal. These results establish SNST as a cross-channel scattering-based connectivity descriptor that recovers amplitude-envelope and cross-frequency coupling structure systematically, applicable to any multichannel EEG analysis where amplitude-domain inter-regional dependence is of interest.

Score matching controls average error under the forward marginals, but a discretized reverse-time sampler evaluates the learned score along its own trajectory. We show that small forward-marginal error does not guarantee numerical stability. We construct a single smooth score field with arbitrarily small forward-marginal $L^2$ error. The learned reverse-time process is nonexplosive, has moments of every order, and can be arbitrarily close to the exact reverse-time process in path-space total variation. Yet its Euler–Maruyama discretizations converge in probability while every positive moment diverges. Thus weak convergence can hold even though every Wasserstein distance $W_p$, $p\ge1$, diverges. The same failure can occur within one fixed finite neural architecture. We construct a family of bounded, globally Lipschitz denoisers for which both the forward-marginal error and the path-space total variation distance tend to zero, while their Euler–Maruyama endpoints diverge in every $W_p$. For compactly supported data, we also give a simple positive result. Projecting the learned denoiser onto a known bounded closed convex set containing the support preserves pointwise accuracy, gives grid-uniform moment bounds, and yields Wasserstein convergence under mild local regularity. Experiments with a small fixed DiT-style network show large growth along rare numerical trajectories and its suppression by denoiser projection, while overall trajectory errors remain small.

N-of-1 trials are patient centered randomized controlled trials. Although the primary goal of N-of-1 trials is to obtain the results for each patient separately, pooling the results across patients also has relevance. In this paper, we present two analytical strategies to pool the results across N-of-1 trials, when the main outcome for each patient is a binomial variable. Our first method takes into account the extra-binomial variation, while as the second approach takes into account hierarchical clustering in addition to overdispersion. We illustrate the methods using real data analysis and compare the methods using simulations.

The Fay-Herriot (FH) model is widely used in official statistics to produce reliable estimates for domains with small sample sizes. In the classical FH model, the sampling variances are treated as known, even though they are typically estimated from the data. In practice, these variance estimates can be highly variable. To address this issue, practitioners often use Generalized Variance Functions (GVFs) to borrow strength across areas and stabilize estimation. In this work, we propose a new Bayesian model to improve the posterior estimation of the sampling variances in Small Area Estimation (SAE). Our proposed model incorporates Global-Local (GL) priors to improve the level of shrinkage toward the posterior estimates obtained with the GVF function. We study the theoretical properties of the proposed model and develop adaptive Markov chain Monte Carlo (MCMC) algorithms to address computational challenges arising from conditional distributions involving Gamma functions. The performance of the proposed Bayesian model is investigated through simulation studies and compared with competing approaches. Finally, we implement our proposal in two real applications by estimating the Corn production from the U.S. Department of Agriculture and the Prevalence of the Educational Attainment Index at the municipality levels in Colombia.

The collective risk model represents the aggregate loss of an insurance portfolio as a random sum of individual claim severities. When claim counts and severities are dependent, the claims pooled across policies are no longer a sample from the marginal severity distribution. We show that their empirical distribution converges to the law of an arbitrary observed claim, a size-biased mixture of the conditional severity distributions, so any procedure that fits the severity margin directly to pooled claims is inconsistent in general. The same result identifies the distribution that the pooled claims do sample, and we build a composite likelihood estimation procedure on that distribution. We establish consistency and asymptotic normality, with Godambe information in which the policy, rather than the claim, is the sampling unit. In a Sarmanov collective risk model, the observed-claim density and the aggregate mean are in closed form. A simulation study measures the bias of naive pooled-severity fitting, its correction by the composite likelihood, and the coverage of the policy-level standard errors.

The problem of recovering an (approximately) low-rank Hermitian matrix $\pmb{M}_0 \in \mathbb{C}^{n \times n}$ of rank $r$ from quadratic sampling matrices of the form $\{\pmb{a}_k \pmb{a}_k^*\}_{k=1}^m$ arises in a variety of applications, including phase retrieval. To obtain rigorous recovery guarantees, the sampling vectors $\{\pmb{a}_k\}_{k=1}^m$ are typically modeled probabilistically. However, most existing theoretical results rely on Gaussian or sub-Gaussian assumptions, which may not accurately capture practical data models. In many applications, sampling vectors exhibit heavier tails, while theoretical understanding in such regimes remains scarce. In this paper, we bridge this gap. We show that two widely used convex approaches, nuclear norm minimization and semidefinite-constrained empirical risk minimization, achieve uniform, stable, and robust recovery under the mild assumption that the entries of the sampling vectors have only finite $4+\delta$ moments, with the optimal sample complexity $m = \mathcal{O}(rn)$ up to moment-dependent constants. The two main ingredients of our analysis are moment estimates for quadratic forms established via decoupling, together with recent advances in covariance estimation in heavy-tailed settings. As byproducts, we also establish the optimal sample complexity for low-rank matrix recovery under complex projective $4$-design sampling, thereby improving upon previous results, and obtain stability guarantees for phase retrieval under similarly weak moment assumptions.

Economic policies rarely affect only their direct targets. To study these spillovers, researchers summarize who else was treated with a simple exposure measure, such as the share of treated neighbors within a radius. But for many settings, economic theory provides little guidance on choosing the functional form (e.g., ring) of that measure or its parameters (e.g., radius). We show that the data can inform both choices. Correctly specified exposure measures imply orthogonality conditions that can be used for both estimation and testing. We establish consistency and asymptotic normality of the resulting estimator under spatial and network dependence in a design-based framework, with all randomness arising from treatment assignment. We then characterize the efficient moment conditions. Applied to two large-scale anti-poverty programs, the framework supports some prior radius estimates but rejects others. In the latter case, the revised radius yields substantively different policy-effect estimates.

Online change-point detection in high-dimensional nonlinear time series faces two challenges. The underlying distributions are difficult to model, and state changes are difficult to characterize. We propose persistent Laplacian cumulative sum (PL-CUSUM) to address these challenges. PL-CUSUM maps delay-embedded sliding windows to point clouds. It extracts persistent Betti vectors and the positive spectra of persistent Laplacians from the same Vietoris-Rips filtration. A ridge-whitened projection converts these features into a scalar score. Page's CUSUM recursion then accumulates this score over time. The positive spectra capture within-scale connectivity and geometric information that persistent Betti vectors do not record. Under a finite-support local model, we prove that the oracle upper bound on detection delay and a local minimax lower bound have the same order. This common order is the logarithm of the average run length constraint divided by the squared ridge-whitened separation. We also establish finite-horizon false-alarm and expected-delay bounds for plug-in whitened scores under weak dependence within each state. The method has two phases. Phase I estimates the projection parameters, selects the feature configuration, and calibrates the control limit. Phase II updates the resulting CUSUM statistic online. Experiments on simulated and real monitoring data show stable false-alarm control and competitive detection performance.

We study the problem of matching two sets of noisy feature vectors when underlying true features are related by an unknown scalar affine transformation. Our method comprises two primary steps. First, we standardize the feature vectors to estimate the unknown scalar affine transformation. Subsequently, we estimate the permutation by minimizing the Least Sum of Logarithms (LSL) between two sets of observations using the estimated transformation. Our main result shows that the unknown permutation can be perfectly recovered given that the minimal separation distance of true feature vectors scales as $\sqrt{\rho_\sigma} \vee (d\log n)^{1/4} \vee \sqrt{\log n}$, where $d$ is the ambient dimension, $n$ is the sample size, and $\rho_\sigma$ is the maximal ratio of noise magnitudes. Interestingly, the obtained rate, under mild heteroscedasticity, coincides with that of the non-affine setting. We additionally demonstrate that there exist configurations requiring a larger minimal separation distance for perfect recovery. The latter makes the matching problem more challenging from minimax perspective compared to the non-affine setting. Consequently, we show that in the problem of feature matching, standardizing the data implicitly estimates the scalar affine parameters. As part of our analysis, we prove non-asymptotic concentration bounds for the affine parameter estimators in the presence of heterogeneous noise magnitudes.

Functional dependence measures have become an important tool in the analysis of nonlinear time series and are typically formulated with respect to a given innovation representation of the process. This note points out that the probability space on which such representations yield the expected memory loss properties may not always coincide with the natural dynamical probability space of the model. We exhibit classes of uniformly ergodic autoregressive processes for which the behavior of the natural innovation coupling undergoes a qualitative transition as the model parameter varies. For this family of models, this transition coincides with a change in the sign of an associated Lyapunov exponent. In particular, a positive Lyapunov exponent may prevent the forgetting of initial perturbations along trajectories driven by the same innovations, despite uniform ergodicity of the associated Markov chain. These observations highlight the importance of carefully specifying the underlying probability space when interpreting or applying functional dependence measures.

Suppose we observe two sets of $n$ Gaussian vectors in $\mathbb{R}^d$, with the promise that, after applying a permutation of $[n]$ and a rotation of $\mathbb{R}^d$, the two sets are $\rho$-correlated. The Procrustes matching problem asks us to recover the unknown permutation of $[n]$ that aligns the two sets. The problem is well-studied in the low-dimensional regime $d=O(\log n)$, but the high-dimensional regime $d\gg \log n$ has remained largely uncharted: prior matching guarantees require nearly perfect correlation $\rho=1-o(1)$, even for information-theoretic recovery. Our main result is a polynomial-time algorithm for exact recovery at constant correlation. The algorithm works by computing and comparing weighted counts of a specially chosen family of “wide” trees. So long as $d\ge \mathrm{polylog}(n)$, the algorithm succeeds with high probability for any $\rho^2>\sqrt{\alpha}$, where $\alpha\approx 0.338$ is Otter's tree-counting constant. We complement this algorithmic result with an improved information-theoretic guarantee, showing that exact recovery is possible when $\rho^2 \gtrsim \max\{\log n/d,\sqrt{\log n/n}\}$. We also carry out a low-degree advantage calculation, which suggests that the condition $\rho^2 > \sqrt{\alpha}$ is necessary for any tree-counting algorithm.

Automated market maker (AMM) fee rules are often evaluated by liquidity-provider (LP) welfare, but that objective mixes fee revenue, adverse-selection loss (loss-versus-rebalancing, LVR), routing response, and liquidity supply. Fixed-fee Uniswap v3 history cannot separate these channels or identify counterfactual trader-facing dynamic-fee rules. Real fee-related variation nonetheless exists: the Uniswap protocol-fee switch cut LP take-rates with tier-differentiated intensity while leaving trader-facing fees unchanged. Using a pre-specified matched-overlap event-study difference-in-differences design, we estimate the liquidity-supply response to take-rate cuts, the kernel K_L that simulator-based fee-controller evaluations routinely freeze, while reconstructing treatment, event time, unit roles, and outcomes from public logs into a frozen, hash-checked panel before any estimate. We detect no large short-run average response in active liquidity or local depth; LP participation and composition, more precisely estimated, likewise show none, so the result is a non-detection at the design's resolution rather than a precise zero. Token-1 volume and native fee income fail the parallel-trends gate and are reported descriptively. A channel-admissibility audit delimits the estimand: the LP-side response K_L is design-based, while trader-facing dynamic-fee protection is a model-conditioned boundary, not a second estimand.

The link between descriptive and substantive representation is well established in the literature but is hard to trace historically, where class records are thin. We introduce a replicable enduring-elite surname measure, pairing a contemporary socioeconomic criterion with historical elite registers, and apply it across the Chilean Congress, 1834-2020. Against a dynamic population reference built from 22.65 million birth registrations, the enduring-elite share of Congress falls from about half in the 1860s to about 12% in the 2010s, with a sharp drop of 11 to 13 points around the 1925 constitutional reform. In 1910-1950, composition co-moves with the legislative agenda, net of party: common-surname legislators emphasize labor foremost, elite legislators a statecraft agenda of defense, foreign affairs, and administration. Across this window, who sits in Congress moves together with what Congress attends to.

Producing analysis weights for a complex survey requires a sequence of hierarchical adjustments (resolving unknown eligibility, dropping out-of-scope units, restoring within-household selection, correcting for nonresponse, and calibrating to known population totals), after which design-consistent variances must account for the fact that several adjustments were themselves estimated from the sample. Existing R tools cover parts of this workflow, but none expresses the whole cascade as a single auditable object, nor propagates the variability of every stage into the replicate weights. We present weightflow, a dependency-free (base R) package that builds survey weights through a declarative, pipeable, tidymodels-style API: a recipe is defined lazily as a chain of step_*() adjustments and estimated with prep(). Separating definition from application makes the process reproducible and auditable, and lets a rescaling bootstrap and a delete-a-PSU jackknife re-apply the entire recipe on each replicate, so the replicate weights carry the uncertainty of every estimated stage, not only of the final calibration. The package implements raking, post-stratification and linear/GREG calibration (with bounds, ridge penalisation, and domain-partitioned and integrative variants), model-assisted (Wu-Sitter) calibration, weighting-class and machine-learning response-propensity adjustments with cross-fitting, and representativity (R-)indicators. Weights and replicate weights bridge to the survey and srvyr packages for design-based inference. We validate the calibration and variance results against survey, and illustrate the full cascade on a bundled sample and on real household-survey microdata (the Uruguayan ECH), where it recovers a known poverty rate with design-based uncertainty.

Hansen and Pitts (2006) introduced the problem of nonparametric estimation of the service-time distribution of an M/G/1 queue observed through its workload process at discrete times $t=0,1,\ldots,n$. Despite its seemingly simple formulation, obtaining an estimator with sharp risk guarantees for this observation model has remained an open challenge for nearly two decades. In this paper, we construct such an estimator and prove that its $L^1$-risk is $\mathcal{O}\left(\frac{\log n}{\sqrt{n}}\right)$ as $n\to\infty$. Remarkably, this nearly parametric convergence rate is achieved without assuming stationarity, stability, or knowledge of the arrival rate. Our approach is based on a two-stage screening procedure that uncovers a hidden conditionally independent compound Poisson structure within the dependent workload observations. This probabilistic reduction transforms the original estimation problem into a classical decompounding problem, making it possible to leverage existing nonparametric estimation techniques despite the complex dependence induced by the reflected workload process. More broadly, we hope that the proposed screening methodology will provide a useful framework for statistical inference from dependent stochastic systems beyond the classical assumption of stability.

We develop confidence regions for Reeb graphs from finite samples using the interleaving distance. Given a point cloud equipped with a filter function, we construct a finite proximity graph, extend the filter linearly, and use the Reeb cosheaf of the resulting filtered graph as the primary estimator. Mapper graphs are then treated as controlled cover-based coarsenings of this estimator, separating the statistical approximation problem from the visualization problem. We prove stability bounds for the Reeb estimators obtained both using intrinsic and extrinsic metrics, the latter under positive-reach assumptions, and derive interleaving-distance confidence regions from either \((a,b)\)-standard sampling assumptions or subsampling-based Hausdorff scale estimates. We also compare this object-level metric viewpoint with persistence-based guarantees by showing that the extended-persistence pseudometric is bounded by twice the interleaving distance, with sharp constant \(1\) for the \(H_0\)-related components. Numerical experiments illustrate how statistically significant features can be identified and then projected to Mapper graphs for interpretation.

In this paper, we study quantile-based distributional reinforcement learning from the perspective of statistical efficiency. We focus on distributional policy evaluation, whose goal is to characterize the return distribution, namely the distribution of discounted cumulative rewards under a given policy. To obtain a finite-dimensional representation of the return distribution, we consider the quantile fixed point $\eta_m$ induced by the quantile-projected distributional Bellman equation. Assuming access to a generative model, we construct an estimator $\eta_m^{(n)}$ based on an empirical Markov decision process. For a fixed number of quantiles $m$, we establish a non-asymptotic error bound for $\eta_m^{(n)}$ and $\eta_m$ under the supremum $W_\infty$ metric, showing that the estimation error scales as $\widetilde{O}(\sqrt{m/n})$ with respect to $m$ and $n$. This implies that the quantile-based distributional policy evaluation problem can be solved with sample efficiency, achieving the optimal parametric $\sqrt{n}$ convergence rate. We derive the asymptotic distribution of the quantile parameters $\sqrt{n}(\theta_m^{(n)}-\theta_m)$ and characterize the semiparametric efficiency bound, which is attained by our estimator. Beyond the fixed-dimensional setting, we investigate the asymptotic regime in which the number of quantiles diverges. We characterize the limit covariance structure and show that it matches the semiparametric efficiency bound of the nonparametric model for distributional policy evaluation, showing that quantile-based estimators remain asymptotically efficient in the infinite-dimensional limit. Finally, we establish a Berry–Esseen theorem for smooth functionals $\sqrt{n}(\eta_m^{(n)}(s)-\eta_m(s))f$, thereby providing a foundation for statistically valid inference on functionals of the quantile-projected return distribution.

We propose a unified Bayesian framework namely robust DAG-Cholesky horseshoe (R-DACH) for joint directed acyclic graph (DAG) structure learning and precision matrix estimation in the high-dimensional proportional asymptotic regime $p/n \to c \in (0,\infty)$, under the scale mixture of normal errors. The construction places a global-local horseshoe-type prior directly on the strictly lower-triangular entries of the modified Cholesky factor of the DAG-Markov precision matrix, so that sparsity in the Cholesky parameters induces a coherent parent-set selection consistent with a topological ordering of the variables. A per-observation inverse-gamma scale mixture yields automatic robustness to heavy-tailed and contaminated observations and admits Student-$t$, Laplace, and slash distributions as special cases. We design a partially-collapsed blocked Gibbs sampler that traverses the joint space of orderings, sparsity patterns and continuous parameters. Simulations across $(n,p)$ configurations with $p$ up to several hundreds confirm the theoretical rates and demonstrate substantial gains over graphical-horseshoe, DAG-Wishart, and PC-based competitors under contamination. An application to RNA-seq gene-expression data from The Cancer Genome Atlas reveals biologically interpretable regulatory structure that competing methods fail to recover.