Papers

Bayesian optimization uses posterior uncertainty to decide where to explore, but this uncertainty is not spatially neutral. On bounded domains, Gaussian-process posterior variance is typically larger near faces, edges, and corners, giving variance-based acquisition functions a systematic boundary bias. This spatial preference is usually inherited from the surrogate rather than chosen as part of the exploration strategy. We separate these roles by keeping the predictive model fixed and modifying the spatial preference of the exploration term. We control the preference with spatial weight fields and amortize their construction with neural surrogates. Controlled experiments show that the Gaussian process' default boundary bias changes which parts of the domain are explored and which optima are reached. We compare three alternative exploration principles and find, across standard benchmark functions and stochastic robot-pushing tasks, that they can often lower regret, particularly early in the optimization. The spatial structure of posterior uncertainty is therefore not a neutral default: it shapes the exploration strategy and can materially affect optimization performance.

When an outcome is missing not at random, identification is commonly restored by a shadow variable: a covariate that predicts the outcome but is excluded from the response mechanism. The literature treats the first requirement, relevance, as a yes-or-no condition. We make it continuous, parametrizing it by a scalar that scales the shadow variable's entire contribution to the moment conditions, and ask what inference costs as that scalar vanishes. One functional of the inverse weight governs the answer. We characterize it exactly, show that it fixes the effective sample size, and use it to sharpen the rates of the estimator that avoids the weight. The characterization is concrete. Whether the functional saturates or diverges turns not on boundedness of the outcome but on whether the tilting variable has an atom at its essential supremum: a binary outcome saturates, a uniform one does not. It also bounds how far a scanned set reaches, and that reach grows only as the square root of the logarithm of the sample size, so a thousandfold increase buys one unit of tilt while the set to be scanned widens without limit.

Human reference labels used to evaluate medical AI can vary across readers and annotation workflows. We quantified agreement of a fixed pathology AI system with a future qualified reader without treating inferred consensus as error-free biological truth. Linked Evaluation and Unbiased NuCLS multirater data were analyzed using a fixed external HoVer-Net decision vector. The selected frame contained 1,144 anchors, 52 fields of view, five slides/patients, and five primary readers. The target was fixed-AI agreement with a future reader, with descriptive contrast $\Delta=\theta_E-\theta_U$. Analyses included raw agreement and iMRMC benchmarks, reader-only and shared-detectability Bayesian models, structural sensitivities, adversarial simulations, and predictive assessment. Among 5,188 reader-anchor pairs observed in both conditions, 19.62% changed AI-agreement state, while raw agreement was 0.5466 in Unbiased and 0.5473 in Evaluation. The reader-only model gave $\Delta=-0.0008$ (95% CrI, $-0.0033$ to $0.0008$), and the shared-detectability model gave $\Delta\approx0$ ($-0.0140$ to $0.0127$). Across structural perturbations, absolute agreement ranged from 0.420 to 0.646 in Unbiased and 0.409 to 0.640 in Evaluation. Frozen adversarial simulations met the comparative-stability criterion (median $R=0.295$; 14/16 model-scenario cells had $R<1$). Held-out-AI prediction supported shared predictive information, but varied by slide composition. The Evaluation-minus-Unbiased contrast was comparatively stable, whereas absolute agreement was structurally sensitive. These findings do not establish equivalence, causal assistance effects, biological accuracy, or new-patient transportability. Evaluations using variable human references should report reader population, selected frame, observation mechanism, and structural sensitivity.

The flux variability of Active Galactic Nuclei (AGNs) is commonly modeled using the damped random walk (DRW), but increasingly rich multiband time-domain surveys motivate more flexible stochastic models that characterize variability beyond a single timescale and exploit dependence across bands. Multivariate continuous-time autoregressive moving average (MCARMA) processes extend the multivariate DRW to higher-order stochastic dynamics and provide a matrix-valued power spectral density for characterizing marginal variability and cross-band dependence. In a carefully selected Sloan Digital Sky Survey (SDSS) Stripe~82 sample, higher-order models are selected for 109 of 134 AGNs (81.3%), with substantial rather than marginal support in most objects; convergence diagnostics indicate that this fraction is conservative. The fitted models also reveal strong cross-band dependence that generally decreases with wavelength separation. Motivated by this richer structure, we investigate its recoverability under Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) observing conditions. Model-order recovery depends strongly on the underlying stochastic behavior as well as source brightness and observing cadence. Recovery is high across many of the simulated settings, but deteriorates when the characteristic features of higher-order variability become difficult to distinguish from the multivariate DRW. Additional simulations show that this loss of distinguishability can persist over a range of underlying variability behaviors. These results show that stochastic structure beyond the multivariate DRW is prevalent in our selected SDSS sample, while its recoverability with LSST depends on both observational conditions and the form of the underlying stochastic variability. The methodology is implemented in the Python package mcarma.

LASSO is one of the most popular regularization methods for high-dimensional regression models. While cointegration is common in practical applications of predictive regressions with many regressors, the convex penalty of LASSO overly suppresses the coefficients of cointegrated regressors and therefore fails to deliver consistent parameter estimation. We employ folded concave penalties (FCP), including the well-known SCAD and MCP, to address this issue. Using FCP, we establish two asymptotic results: (i) a local oracle property for the penalized regression, and (ii) a global oracle property for regressions with a penalty-level constraint. In particular, the latter guarantees a global solution that consistently estimates the true coefficients and enjoys oracle prediction performance. Simulation studies and an application to predicting the euro area unemployment rate show that FCP improves variable selection and out-of-sample prediction relative to LASSO.

This paper develops strong Gaussian approximations (KMT type) for constant learning rate stochastic gradient descent (SGD). We work with heavy tailed noise and assume that the number of iterates grows and the learning rate shrinks simultaneously. Our work constructs simple Gaussian processes on the same space and bounds, in probability, the maximum difference between the SGD errors and their respective limits. While there exist central limit theorems for constant rate SGD centered at the optimum, these only provide individual iterate inference. By contrast, our quenched KMT approximations allow for finite sample, simultaneous, whole-trajectory-level inference and imply stochastic process convergence (Functional Central Limit Theorems (FCLTs)). For accumulated errors centered at the true optimum, we first obtain a KMT-type approximation bound when the learning rate scales as $\gamma_n=\frac{c}{\sqrt n}+o(n^{-1/2})$ for a constant $c \in \mathbb{R}^+$ and $n$ iterations. The result reveals a trade off between controlling the bias of SGD and the initialization effect, with no learning rate scaling generally making both contributions vanish in the accumulated error limit. To eliminate these contributions, we consider a Burn in scheme and Richardson-Romberg extrapolation and prove respective KMT and FCLT theorems. We also prove a finite-horizon KMT-type approximation and FCLT for the rescaled iterate trajectory to a stationary Ornstein–Uhlenbeck (OU) process. We consider applications in risk control and simultaneous inference and perform simulations to illustrate our results.

The Gromov-Wasserstein (GW) distance provides a principled framework for comparing two probability distributions supported on different metric spaces by matching their individual intrinsic geometric structures. In practice, the GW distance is estimated from finite samples and often incorporates regularization for statistical and computational tractability. This motivates a quantitative understanding of the stability and estimation error of regularized GW distances. In this work, we study a general family of divergence regularized GW distances, including entropy and Tsallis regularizations over the full range of powers, from a nonasymptotic perspective. Specifically, we establish Lipschitz continuity of the value of regularized GW and Hölder continuity of its optimal coupling under Wasserstein perturbations of the input probability measures. Moreover, we derive the parametric rate of convergence for the empirical regularized (2,2)-GW distance on compactly supported Euclidean spaces.

In observational network data, the adoption of a treatment and the formation of social ties often depend on the same characteristics that affect outcomes. We propose a generalized propensity score for each unit's effective treatment, built from a model for treatment assignment and a dyadic model for network formation, which may include latent node heterogeneity. The score is marginalized over the network given covariates rather than conditioned on the realized network, and the estimand is a population effect: the average effect of an exposure level for a unit drawn from the population of units and networks. We relate this estimand to those of methods that condition on the realized network and show what marginalizing costs and buys. It requires network formation to be unconfounded given observed covariates, a stronger assumption than unconfoundedness of treatment given the network. In exchange, it defines effects at exposure levels that poorly connected units cannot reach, and at levels every unit can reach it gives an estimator whose asymptotic variance is no larger than the network-conditional one. We combine the score with an augmented inverse probability weighting estimator, prove asymptotic normality as the number of networks grows with bounded network size, discuss what carries over to a single large sparse network, and apply the method to the effect of a political communication technology on political participation in Uganda.

Determinantal point processes (DPPs) provide tractable distributions for nite random sets, but the usual symmetric-kernel specication imposes nonpositive pairwise dependence. Statistical research on nonsymmetric DPPs, which can accommodate positive dependence, remains limited. We study identication and estimation for this class. First, we establish the equivalence of two sucient conditions for a nonsymmetric kernel to dene a valid DPP. Second, we develop the identication theory for nonsymmetric DPPs. Third, we develop cardinality-truncated maximum-likelihood estimation and estimation by averaged stochastic gradient descent. Under suitable identication and regularity conditions, we establish consistency and asymptotic normality for both estimators. Simulations examine identication, the eects of truncation, and the nite-sample performance of the estimators. We illustrate the methods with two cyber applications at markedly dierent sampling frequencies: weekly sector-by-breach data and second-by-second enterprise authentication data with more than ve million set observations.

The Copula Active Subspaces (CAS) method of Part I estimates a non-Gaussian observation-noise density $\pi_{\boldsymbol{n}}$ from residual samples. Its KL divergence $\mathcal{D}_{\mathrm{KL}}(\pi_{\boldsymbol{n}}\,\|\,\widehat\pi_{\boldsymbol{n}})$ admits an exact decomposition into marginal, projection and reduced-density terms and bounds the posterior KL divergence averaged over the true predictive distribution. We derive a priori bounds for these terms, together with estimates or bounds from residual samples. Weighted nearest-neighbour entropy estimates are used for the marginal and reduced-density terms, with asymptotic standard errors in exact rank coordinates. For the projection error, the score is fitted on one subsample, and its mean squared magnitude outside the retained subspace is estimated on the other. Corrections for score and rank-transform errors then give a theoretical upper bound with probability at least $1-\delta$. Tests on three known noise laws in dimension $d=20$ use the estimated rank transform. In $11$ of $12$ combinations of example and sample size, the sum of the two estimated terms and the trace bound exceeds the reference divergence after averaging over the estimated models. In these tests, score-estimation error is the largest contribution to looseness. The overestimate from the Gaussian log-Sobolev inequality alone is at most a factor of $1.98$.

We study Wasserstein-$p$ goodness-of-fit and two-sample testing under differential privacy. For general compact metric spaces, we establish upper and lower bounds on the minimax separation radii in terms of covering numbers, matching up to logarithmic factors. On the hypercube $[0,1]^d$, these bounds match up to constant factors for two-sample testing in all dimensions and for goodness-of-fit testing whenever $d\neq 4p$, thereby closing logarithmic gaps in existing non-private results. At the critical case $(d,p)=(4,1)$, we further show that a logarithmic factor is unavoidable by establishing matching minimax bounds for non-private Wasserstein-$1$ goodness-of-fit testing under a bounded-density assumption. Finally, for unbounded domains, we establish minimax-optimal separation rates up to logarithmic factors under an exponential tail condition.

As watermarking of LLM-generated text enters deployment, users increasingly need to check text whose source is unknown. Existing detectors read the watermark through a specific tokenizer or embedding pipeline, and even if all providers adopted the same watermarking method, differences in these representations would leave their watermarks incompatible, and the text alone does not reveal which representation produced it. We propose ByteMark, a unified watermarking scheme that defines its signal on the UTF-8 bytes of published text. Models with different tokenizers embed this signal under one public specification and key, and a single detector verifies the text from the text and key alone. When the tokenizer or embedding model of a baseline detector is replaced, its detection degrades sharply, whereas ByteMark, which uses neither, retains its detection rate; ByteMark also remains detectable when participating models continue or rewrite the text.

Modern learning problems increasingly leverage data from heterogeneous sources to improve robustness and generalizability under distribution shift, often with high-dimensional parameters admitting low-rank structure. Existing group distributionally robust optimization (group-DRO) approaches optimize against adverse source mixtures, but extending this principle to low-rank matrix estimation with a broad class of losses leads to a challenging nonconvex-concave problem for which methods and theory remain limited. We propose Prior-Informed generalizable Low-rank representation Learning via Adversarial Reweighting (PILLAR), which smooths the group-DRO objective by regularizing adversarial source weights toward a user-specified reference distribution while learning a shared low-rank representation. By tempering extreme adversarial reweighting, PILLAR improves estimation stability and reduces the potential conservativeness of group-DRO. Using Burer-Monteiro factorization, we optimize the nonconvex-concave objective by bifactored gradient descent with mirror ascent. We establish linear convergence to a statistical neighborhood of the population saddle point and derive an error bound that separates low-rank estimation from adversarial reweighting. We further show that PILLAR improves on vanilla group DRO when the latter's excess conservativeness dominates source heterogeneity and prior misspecification. Simulations and an electronic health record study of suicide-attempt risk demonstrate improved worst-source and transported predictive accuracy.

What does it mean to measure the replicability of scientific results? Is it merely determining whether results are consistent across studies? Much of the existing literature on replicability draws conclusions exclusively from the comparison of study results, in terms of observed estimates, uses a variety of different metrics, and ultimately collapses the comparison into a binary verdict of success or failure. This fails to differentiate between repeated errors and corroboration of the shared generative effect. We develop a unified Bayesian framework in which distinct but related notions of replicability are formally defined and jointly embedded within a single inferential structure. A single model answers four complementary questions: agreement among independent studies; consistency between the individual studies and the shared generative effect; agreement between an existing body of evidence and a new study from the same super-population; inference on the shared generative effect. In practice, we suggest quantifying replicability through posterior probabilities relative to regions of practical relevance. We illustrate these ideas in an eQTL (expression quantitative trait locus) analysis assessing the replicability of a regulatory effect across human ancestral groups.

Neural learners for classification and regression are often criticized as black boxes. We propose a probabilistic decomposition of neural learners that accounts for heterogeneous predictability across the feature space. A training procedure combining the EM algorithm, stochastic gradient descent, and maximum likelihood estimation is developed to jointly learn the predictable component and an input-dependent reliability score. We show that the optimal decomposed learner recovers the response distribution within the predictable component, while the reliability score characterizes the predictability of individual observations. We further extend the framework to incorporate localized feature transformations. Experiments on synthetic and real-world datasets demonstrate the effectiveness of the proposed framework.

Quasi-periodic oscillations are characterized by fitting Lorentzian functions to a nonparametric estimate of the power spectral density, a two-step procedure that separates estimation from modeling and does not accommodate measurement error. We develop a parametric alternative based on autoregressive moving-average (ARMA) processes observed with measurement error. We derive an exact decomposition of the ARMA spectral density into real-valued components, one for each class of autoregressive roots under complex conjugation, plus a polynomial term when the moving-average order is at least the autoregressive order. Each component is a Wrapped Cauchy kernel, and for complex roots it is accompanied by a conjugate-Poisson term that makes the peak asymmetric and displaces its maximum from the frequency given by the argument of the root. We obtain closed-form expressions for the location, height, width, and asymmetry of each peak, and a matching variance decomposition. Parameters are estimated by maximum likelihood through the Kalman filter, and all descriptors follow by plug-in. Simulations show that the closed-form peak location is recovered without appreciable bias. An X-ray binary light curve illustrates the methodology.

Quantile regression has received considerable attention in high dimensional methodology due to its robustness and ability to accommodate heteroskedasticity. Despite this attention, current inferential methods for high dimensional quantile regression have significant limitations—only testing significance, testing at a single quantile at a time, or assuming homoskedasticity. To generalize the inferences that can be drawn in high dimensional quantile regression, we introduce partial penalized Wald and rank score tests for testing linear hypotheses uniformly over a set of quantiles $\mathcal T$ under a heteroskedastic model. We propose globally adaptive local linear approximation (LLA) algorithms for folded concave penalized estimation and show that our global LLA algorithms recover “testing oracles” for the full and reduced models with high probability. We then show that our global LLA estimators achieve the $\sqrt{s/n}$ oracle estimation rate uniformly over $\mathcal T$, establish Bahadur representations, and derive asymptotic distributions for the test statistics under the null and local alternatives. We develop a dual formulation incorporating adaptive penalty weights for the rank scores and use the strong oracle property to establish their oracle equivalence. Using empirical process techniques, we derive reference distributions based on suprema of squared Bessel processes, extending partial penalized inference beyond pointwise chi-square calibration.

Treatment decisions are often sequential and depend on information absent from the recorded history. Latent sequential ignorability models such assignment by assuming that baseline covariates and principal strata, defined by joint potential intermediate responses, account for confounding of later treatment. For two treatment periods with a binary intermediate response, we derive sharp population and principal effect bounds under a specified association between the potential responses and show that the two targets can require different optimizations to obtain sharp bounds. For population effects, optimization over latent treatment probabilities reduces to their feasible endpoints. Stronger selection between response types can be offset by making their outcome distributions more similar, preserving both the observed distribution and population marginal outcome distributions. This replacement can change effects within principal strata, whose sharp bounds may instead require interior treatment probabilities. Restrictions on differences in conditional effects can tighten the bounds without equating outcome levels. For binary outcomes, we give an exact finite algorithm that includes these interior assignments. We also develop inference under smoothness and estimation conditions. We illustrate the bounds in a study of gym attendance after an incentive program.

Sharing time-to-event data can leak sensitive information about individual participants in clinical studies. To protect the privacy of data, differential privacy (DP) framework offers rigorous protection for data sharing. The local variant of differential privacy is particularly appealing because each record is perturbed before it leaves the data owner or contributor. However, the injected perturbation distorts the true distribution of the source data, yielding low utility and invalid statistical inference with the released data. In this work, we propose a nonparametric estimator for the true survivor function via the deconvolution kernel function with locally differential private (LDP) right-censored survival data in which both the event times and the event indicators are privatized. Our estimator simultaneously de-biases the perturbation applied to the times and the errors induced by flipping the event indicators due to LDP. Theoretically, we establish the uniform consistency and point-wise asymptotic normality of the de-biased distribution estimator, and derive its uniform convergence rate together with the optimal bandwidth attaining that rate. Simulation studies and the analysis of a semi-synthetic dataset generated from the METABRIC breast cancer cohort demonstrate the performance of the proposed estimator across privacy budgets, censoring levels, and sample sizes. By the post-processing property of DP, the proposed de-biased estimator inherits the privacy guarantee while restoring valid inference on the survivor function.

We develop a measure-theoretic framework for conditional multivariate total positivity of order two ($\mathrm{cmMTP}_2$), formulated through lattice inequalities without requiring the existence of probability densities. This formulation extends classical density-based notions of positive dependence to general probability measures, including singular distributions. We establish fundamental preservation properties under marginalization, monotone transformations, and weak convergence, thereby providing a unified treatment of dependence structures under distributional limits. We further investigate the behavior of conditional total positivity under Markov-kernel composition and derive characterizations for important copula families. In the Archimedean setting, the defining condition is expressed through log-convexity of generator derivatives, whereas for Gaussian distributions it corresponds to nonpositive off-diagonal entries of the precision matrix, or equivalently, nonnegative partial correlations. To connect these structural results with statistical inference, we introduce empirical lattice contrasts and establish their joint asymptotic normality. The resulting simultaneous testing procedures achieve asymptotic size control and admit a local power analysis. Our results provide a general foundation for density-free characterization, specification testing, and inference concerning positive dependence in multivariate probability models, particularly when conventional density-based methods are inapplicable or insufficient.

The cvxoptdes R-package provides a unified framework to perform convex optimal experiment design and online response optimization for linear, generalized linear and additive, nonlinear, and random intercept mixed models, using exclusively open-source solvers. It supports optimization and augmentation of weighted, penalized and convex constrained continuous designs across a discrete grid of candidate points and rounding to integer-valued designs using point-exchange algorithms or mixed-integer programming for orthogonally constrained designs, including the class of orthogonal minimally aliased response surface (OMARS) designs. Models and design constraints are defined and updated in an object-oriented context with the aim to facilitate the process of sequential (penalized) design augmentation analogous to online design procedures known from active learning and Bayesian optimization.

Augmentation-based self-supervised learning can recover the content shared across views without separating its underlying factors into individual coordinates. We ask whether augmentations can be designed to make those factors identifiable. Our approach allows content factors to fluctuate across augmented views while preserving their mean for each source instance. Each instance receives a fixed, recorded set of augmentation strengths that controls how much its factors fluctuate. When factors respond differently to changes in these strengths, their patterns of variation provide information that distinguishes them. Under suitable decoder and support assumptions, we prove that these differences identify individual factors in a conditional Gaussian model, up to permutation, scaling, and translation. We then develop a learner that uses the recorded strengths to estimate these patterns jointly across instances, enabling learning from few views per instance. Synthetic experiments demonstrate near-perfect factor recovery through nonlinear mixing. Controlled dSprites experiments show improved coordinate alignment using augmentations that target known factors, even without telling the learner which augmentation corresponds to which representation coordinate. These results suggest a route to component identification in self-supervised learning: augmentation design can supply identifying information beyond shared content alone.

Variance-based methods recover a Gaussian linear structural equation model (SEM) with equal or suitably constrained error variances, either by recursively identifying the source nodes with smallest conditional variance (top-down method) or sink nodes with highest conditional variance (bottom-up method). We show that the use of conditional variance is incidental: replacing it with the conditional risk, i.e., the minimal expected loss of linear prediction under the negative log-likelihood of the noise distribution, extends both methods to non-Gaussian noise. Existing methods are the least-squares special case. This yields two main results. First, exploiting the noise distribution improves structure recovery: the matched-loss procedures recover the SEM more accurately than their least-squares counterparts and the LiNGAM methods, which are designed for non-Gaussianity. Second, the two methods, interchangeable in the Gaussian setting, separate under non-Gaussian noise: the top-down method is correct for every noise density with finite differential entropy, whereas the bottom-up method requires a regularity condition and can fail without it, for instance under mean-zero exponential or Cauchy noise. Simulations confirm both results.

This paper analyses spatial trends in armed conflicts in North and West Africa. It shows that all types of violence saw a significant increase in 2025. However, these violent incidents are distributed very unevenly across the region: while conflict is at its lowest in North Africa, the vast majority of violent incidents and victims are concentrated in five West African countries: Nigeria, Cameroon, Mali, Burkina Faso, and Niger. Coastal countries are experiencing a very mixed picture: jihadism is not advancing as a front, from the Sahel toward the Gulf of Guinea, and most of the violence stems from local disputes, ethnic conflicts, and banditry. Second, the spatial indicator developed by the OECD/SWAC to track the evolution of this violence confirms that the conflict landscape is expanding at its margins and intensifying at its center. Even more concerning is that peripheral areas are increasingly being drawn directly into well-established conflict zones, rather than going through a gradual, low-intensity transition phase. The geography of expansion is shifting from a diffuse, low-intensity spread toward peripheral areas to a rapid outward expansion from dense, well-established conflict hotspots. The local indicator used in the report suggests that violence continues to spread to new locations at the margins of established conflict zones. At the heart of these zones, tensions are intensifying relative to their historical baseline levels, and a growing number of regions have experienced clashes for a decade or more. Projections indicate that the intensity of violence is expected to increase in 2026, reaching more than 15,000 incidents and more than 32,500 deaths. Much of the increase is explained by the growing number of conflicts in Nigeria. In the Sahel, fatalities should increase in Mali and decrease in Burkina Faso.

Rating procedures are crucial in psychological and medical diagnosis, educational assessment, peer review, and grant evaluation. A major methodological challenge is disentangling multiple sources of variability in ratings, such as subject ability, rater behavior, and item characteristics. Item Response Theory-based models, most notably many-facet Rasch models and hierarchical rater models, have been widely adopted for jointly estimating subject, rater, and item effects, while retaining high interpretability. We introduce the Dynamic Extended Generalized Many-Facet Rasch model, a general class of Bayesian semiparametric models that integrate and extend existing approaches. Its measurement model allows location and discrimination parameters to depend jointly on item and rater characteristics and incorporates subject-item and rater-subject variation. The dynamic extension introduces a rater-specific cross-lagged effect capturing dependence on the previously rated subject, and a Bayesian P-spline capturing systematic changes in discrimination over the rating sequence, accounting for learning, fatigue, or adaptation effects. Subject heterogeneity is addressed through Bayesian nonparametric mixtures, flexibly learning the latent ability distribution without strict parametric assumptions and enabling model-based subject clustering. We also derive a closed-form, model-based inter-rater consensus coefficient and study its characteristics. We discuss theoretical properties, develop posterior computation strategies, and evaluate the model in simulation studies and two real-data applications.

Generative modelling involves two complementary tasks: learning a probabilistic approximation of an unknown data distribution from observations, and designing a tractable mechanism for sampling from this approximation. We present a stochastic-analysis perspective on score-based generative models, connecting ideas from classical sampling theory, score estimation, and diffusion time reversal. We first explain how score functions link Langevin sampling with score matching, and how Gaussian smoothing leads to denoising score matching and multiscale sampling methods. We then emphasize the distinction between dynamics preserving a fixed distribution and dynamics transporting a prescribed family of marginals. This motivates a path-space treatment of diffusion time reversal. Finally, we connect the ideal reverse diffusion with practical generative algorithms.

We study the convergence of Markov chains in the geometry induced by Hilbertian embeddings associated with positive definite kernels, also called Maximum Mean Discrepancy (MMD). Under a contraction assumption in the kernel-induced geometry, and a standard spectral gap assumption, we establish non-asymptotic, high-probability bounds for the convergence of the empirical occupation measure to the invariant measure. The obtained bounds separate the roles of the two conditions: the spectral gap determines the effective sample size in the leading term, while the kernel contraction enters only through a linear correction. As an application, we study one-step prediction of dynamical systems from a single trajectory, deriving fast learning rates for a nonlinear least-squares estimator. The rates highlight the interplay between mixing properties in the kernel-induced geometry and the classical source and capacity conditions used in the analysis of least-squares.

We study the last iterate of standard asynchronous Q-learning from a single trajectory of a finite Markov decision process. The trajectory follows a fixed behavior policy whose state–action chain is irreducible and may be periodic. For discount factor $\gamma$, write $H=(1-\gamma)^{-1}$, and let $\mu_{\min}$ be the minimum stationary state–action probability and $\mathfrak h_\star$ the largest expected hitting time of the reverse state–action chain. We prove that zero-initialized Q-learning with a constant step size achieves sup-norm error at most $\varepsilon$ with high probability using $\widetilde{O}(H^4/(\mu_{\min}\varepsilon^2)+\mathfrak h_\star)$ transitions, for $0<\varepsilon\le1$. The statistical term retains the quartic effective-horizon dependence of standard Q-learning, while initialization cost is additive in $\mathfrak h_\star$. For aperiodic behavior chains, a comparison with mixing time yields a transient term $\widetilde{O}(t_{\mathrm{mix}}/\mu_{\min})$. An empirical state–action count rule selects a constant step size for a subsequent training run and attains the same rate without prior knowledge of model parameters, target accuracy, or confidence level. The output is the last training iterate, with both stages included in the sample count. The proof uses two deterministic Bellman envelopes to separate statistical error from fixed-policy initialization error and adapts reverse-time hitting-cost compensation from TD learning to state–action propagation. An exploration lower bound holds for arbitrary predictable step sizes. Together with the existing statistical lower bound and its extension to pilot-selected constant steps, it yields matching mean-square rates up to logarithmic factors on the stated large-sample parameter classes.

Many scientific models can be simulated but have intractable likelihoods, making frequentist uncertainty quantification difficult. We construct exact anytime-valid confidence sequences for the parameters of such models. At each observation, fresh simulations from each candidate null conformalize a positive predictable score into a conditional e-value; their product is a test martingale, which we invert over the parameter. Validity holds for any such score, whatever its accuracy. We show that the oracle likelihood-ratio score maximizes expected log growth, which motivates approximating it with an amortized neural ratio estimator. On synthetic benchmarks, an audit of model routing in an LLM API, and microlensing data, the resulting confidence sequences stay valid under continuous monitoring, where fixed-sample methods do not, and neural scores remain competitive with exact likelihood ratios and application-specific scores.

We establish finite-sample guarantees for conditional diffusion models learned from data and implemented with finitely many denoising steps. Our approach aligns the score-training objective with the sampling grid, allowing statistical estimation and numerical discretization errors to be analyzed together. Under a smooth conditional density model, we derive high-probability bounds on the generated distribution's error and quantify the required number of score evaluations. For every smoothness index $\beta>0$, a suitable grid yields averaged conditional Kullback–Leibler (KL) divergence $\widetilde O_p(n^{-2\beta/(d_x+d_y+2\beta)})$ relative to the original target with or without early stopping. Both the full and early-stopped samplers attain average target total-variation (TV) error $\widetilde O_p(n^{-\beta/(d_x+d_y+2\beta)})$. Here $n$ is the sample size, $d_x,d_y$ are the data and conditioning dimensions, and errors are averaged over the conditioning variable. The total-variation bound matches the rate established in existing conditional distribution-estimation theory, up to logarithmic factors.