Papers

This paper studies the optimal design of Type-I generalized progressive hybrid censoring schemes for life-testing experiments. The design problem involves simultaneously determining the inspection time, the guaranteed number of failures, and the progressive censoring scheme. First we develop a cost-constrained optimization framework for determining the optimal censoring scheme. Structural properties of the A-optimality criterion and the experimental cost with respect to the inspection time and the guaranteed number of failures are established. It reveals that they are conflicting behaviors which enables to develop an efficient search algorithm that substantially reduces the computational burden. Building on these theoretical results, a multi-objective optimization model is proposed to simultaneously minimize A-optimality criterion and the experimental cost. A Variable Neighborhood Search (VNS) algorithm is proposed to efficiently determine the optimal progressive removal vector by exploring the feasible design space while avoiding exhaustive enumeration. The resulting compromise designs simultaneously improve estimation precision and reduce experimental cost. In addition, the Shannon differential entropy of the observed lifetime distribution is derived and employed as a complementary information-theoretic measure for evaluating the selected censoring schemes. Numerical studies show that entropy-optimal designs generally differ from A-optimal designs, indicating that Shannon entropy characterizes uncertainty in the observed data rather than estimation precision. The proposed methodology provides an efficient computational framework for optimal life-test design and offers a foundation for future multi-objective optimization incorporating statistical efficiency, experimental cost, and information-theoretic uncertainty.

We propose a novel measure of the discrepancy between two probability distributions $f$ and $g$ on a graph - which we call the diffusion distance - that measures the rate of convergence of $f$ to $g$ under a graph-constrained Markov chain with stationary distribution $g$. As a default choice for this Markov chain, we use the Metropolis-Hastings transition matrix targeting $g$ with proposals given by a random walk on the graph. Our primary case of interest is when the second distribution $g$ is uniform, in which case the diffusion distance becomes a measure of spatial clustering in $f$. Used in this way, (Metropolis-Hastings) diffusion distance to uniformity extends Moran's $I$-type measures of spatial autocorrelation by incorporating global graph geometry rather than just local patterns. Indeed, Moran's $I$, the most well-known measure of spatial autocorrelation, can be viewed as a one-step heuristic for diffusion distance, so long as specific spatial weights are used. We establish theoretical bounds and a stability result for our measure, connecting it to graph spectra and optimal transport. We then turn our attention to outlining a statistical test for spatial clustering using diffusion distance. Under permutation null models, we derive high-probability bounds on diffusion distance underpinned by exact spectral formulas for convergence of distributions, enabling an efficient statistical test for spatial clustering on large datasets. We empirically compare diffusion distance to Moran's $I$ both as a numerical measure and as a statistical test. We show that diffusion distance exhibits higher power on synthetic data using a stochastic block model. Empirical analysis of Black population distributions for 100 U.S. cities shows that diffusion distance detects subtle differences in urban segregation patterns that Moran's $I$ does not.

In this paper, we propose a novel tamed stochastic gradient Hamiltonian Monte Carlo (tSGHMC) algorithm for sampling and stochastic optimization problems with superlinearly growing stochastic gradients. Under a certain continuity in average condition and a strong convexity condition, we establish a non-asymptotic error bound in Wasserstein-2 distance for tSGHMC with the rate of convergence equal to $1/4$. Then, we derive an upper estimate for the associated expected excess risk, which provides a theoretical guarantee for the performance of tSGHMC. To illustrate the effectiveness of the proposed algorithm, we apply tSGHMC to practical examples, including a newsvendor problem and a Conditional Value-at-Risk minimization problem, using synthetic and real-world datasets. Numerical results support our theoretical findings. Furthermore, we compare tSGHMC with its first-order counterpart, namely, the tamed unadjusted stochastic Langevin algorithm. Simulation results demonstrate that tSGHMC achieves lower root mean square error and expected excess risk across a range of tasks.

Retail store locations are strategically selected rather than randomly distributed, potentially inducing preferential sampling when the latent spatial factors governing placement also affect store survival. We propose a Bayesian hierarchical model that jointly combines a log-Gaussian Cox process for store locations with a probit regression for binary survival outcomes. The two components share a Gaussian process spatial effect, with a loading parameter measuring the association between the latent drivers of store placement and survival. To enable efficient inference for approximately 1,000 observations, we use a nearest-neighbor Gaussian process approximation and a Metropolis-within-Gibbs algorithm. We apply the model to 999 mini-supermarkets in Tokyo's 23 special wards, including 897 operating and 102 closed stores, using seven spatial covariates and a 3,471-point integration grid. The estimated loading is close to zero, with its credible interval including zero, providing no clear evidence of residual preferential sampling. Regression estimates are also stable across models with and without preferential sampling. Simulations show that the method can distinguish absent from strong preferential sampling. Proximity to full-scale supermarkets is the most robust predictor of closure risk, consistent with competitive substitution.

We consider the post-detection analysis of change-points for multivariate time series, with the goal of identifying which coordinates are responsible for a detected change. After a change-point has been located by an offline detection algorithm, we propose post hoc statistical procedures to determine whether the change occurs in either of two predefined blocks of coordinates or in both. Our methods rely on two-sample testing procedures with a particular focus on nonparametric tests; we provide theoretical guarantees for Type I error control. Simulations and a real-data experiment demonstrate the strong performance of the proposed procedures.

We study the worst-case false discovery rate (FDR) of the Benjamini-Hochberg procedure for both one- and two-sided Gaussian tests with general correlation matrices. For each setting and given target FDR level $q$, we construct a family of finite Gaussian models whose FDR divided by $q$ diverges as $q\downarrow0$. This disproves any universal multiplicative FDR bound. Moreover, we establish new FDR lower bounds that are independent of the number of hypotheses. For two-sided tests, the worst-case FDR over the number of hypotheses, mean vector, and correlation matrix is at least an explicit ${\ell_{=}}(q)>q$: \[ \ell_{=}(q)=\frac{q\sqrt{\log(1/q)}}{2\sqrt{\pi}}+c_\ell q+o(q), \qquad c_\ell=0.6492828\ldots. \] For one-sided tests, we show a lower bound $\ell_{\le}(q)>q$: \[ \ell_{\le}(q)=\frac{q\sqrt{\log(1/q)}}{\sqrt\pi} +\frac q2+o(q). \] Finally, we prove the leading constants are sharp for a class of correlation matrices for both one- and two-sided tests.

To better forecast mixed-frequency time series, it is the key to choose a suitable way for frequency alignment. However, the existing methods are all limited to linear transformations, and this may overlook the possible nonlinearity, leading to a worse prediction. We alternatively consider a deep neural network for each frequency alignment, and hence a depth-separable neural network. Moreover, a parameter-sharing mechanism is adopted across the alignment at each stage, making possible a deeper network for a large set of higher-frequency predictors. This paper establishes an approximation theory for the proposed depth-separable network, and a non-asymptotic prediction error bound is also derived. Simulation studies demonstrate the finite-sample performance of the proposed method, and an empirical application to forecasting U.S. quarterly macroeconomic variables using monthly and daily indicators, highlights its superior predictive accuracy over existing mixed-frequency methods.

Temporal Knowledge Graphs (TKGs) record how facts evolve over time, but forecasting future events on a TKG remains difficult for three reasons: (i) long-range temporal dependencies are hard to encode; (ii) events on different chains mutually excite or inhibit one another in ways that snapshot-level models cannot express; and (iii) inter-arrival times are heavy-tailed and statistically sparse, so deterministic time predictors are unreliable. We address these three issues with a single framework, the Group Attention Neural Hawkes Process (GAttNHP), built around three matched components. First, a self-attention encoder casts each subject–relation chain as a continuous-time point process and captures the lingering excitation of distant history. Second, a semantic soft-grouping module turns globally learnable Hawkes priors into an analytical cross-attention mask, so chains share excitation patterns through their latent group memberships rather than through exhaustive pairwise computation. Third, a Non-Crossing Quantile (NCQ) regression head replaces mean-based time prediction, providing calibrated, monotonically ordered quantile estimates that remain stable under heavy-tailed inter-arrival distributions. On six benchmark TKG datasets, GAttNHP improves over state-of-the-art baselines on both entity prediction and time prediction, and ablations confirm that its largest gains arise on the long-tail event chains where existing models fail most severely.

Local SGD, also known as Federated Averaging, is a widely used distributed optimization algorithm. Although Local SGD often outperforms alternatives such as Mini-batch SGD in practice, theory still only partially explains when and why local updates help under realistic data heterogeneity. Recent work by [Patel et al., 2025] shows that a bounded second-order heterogeneity assumption captures the efficiency of Local SGD for strongly convex objectives, and conjectures that the same principle extends to the general convex setting. In this paper, we prove this conjecture by establishing an improved convergence guarantee for Local SGD on general convex objectives under bounded second-order heterogeneity. We also improve the best-known lower bounds for Local SGD in this setting, showing that our upper bounds are nearly tight. Together, these results provide a sharper, more fine-grained convergence theory for Local SGD. As a further application of our techniques, we provide a lower bound for serial SGD with replacement, showing how second-order heterogeneity captures the impact of rare high-curvature clients.

We introduce a new equivalence test to show sufficiently good agreement of observed data with a binary generalized linear model (GLM). The test statistic is constructed via the minimum distance method. The test is developed for the important special case where all covariates are categorical. The critical values can be calculated using an asymptotic approximation or by means of bootstrapping. The application of the test to logistic regression is illustrated on two real data sets. The finite sample performance of the proposed test is studied by simulations which are based on these two data sets.

Wiener-Hermite cross-correlation identification represents a polynomial response in the Hermite basis. Under Gaussian excitation the basis is orthogonal and a diagonal rule recovers it exactly; under non-Gaussian excitation the same basis is kept, but its Gram matrix gains off-diagonal terms and the diagonal rule is no longer the population projection. We give the exact finite-order excess $L^2(P)$ risk of this mismatch: a moment quadratic form from two Hankel-Cholesky factorizations and one diagonal solve, at $O(s^3)$ cost from moments to order $2s$. Closed cumulant forms at orders three and four expose which non-Gaussian features drive it; symmetry protects the Gaussian basis only through order two. A bootstrap decides, from data, whether a matched basis is worth building; on a Wiener-Hammerstein benchmark it separates a near-Gaussian channel (penalty $\approx 10^{-4}$) from a skewed output (penalty $0.05$). The computation is a weighted-$L^2$ projection whose core normal-system correspondence is machine-checked in Lean 4.

We propose an operator-split Bayesian learning framework for second-order uniformly elliptic Dirichlet problems with unequal numbers of interior and boundary observations. The data consist of noisy measurements of the source in the domain and noisy measurements of the boundary values. Independent Bayesian neural-network (BNN) priors are assigned to these two quantities, and the resulting product posterior is pushed forward through the elliptic solution operator. We prove that the posterior induced by this construction contracts around the true solution. The contraction radius separates a domain contribution, governed by the second-order elliptic operator, from a boundary contribution, governed by the intrinsic dimension of the boundary. Together with the minimax lower bound of \cite{ZhaoLu2026}, this yields a near-minimax upper bound up to logarithmic factors. Our numerical experiments illustrate the propagation of source and boundary uncertainty and the effects of unequal sampling budgets on the posterior reconstruction.

I argue that multiverse analysis is highly suited to two undesirable uses: abdication of researcher's responsibility for their conclusion and manufacturing of doubt. A review of multiverse analyses published in 2025 provides tentative empirical support that abdication of responsibility is present in the literature and I mention anecdotal evidence that multiverse has been used for manufacturing of doubt about Covid-19 precautions. To mitigate negative effects if multiverse analysis becomes widely used I suggest the community adopts two conventions for evaluating multiverse analyzes: evaluating multiverses by the single worst universe they contain and considering large size of a multiverse as a sign of weakness rather than a praiseworthy achievement.

Online controlled experiments are the gold standard for hypothesis testing in online platforms. Notwithstanding their ubiquity, they are notoriously expensive to run, and issues of variance hamper statistical power in assessing treatment effects. While standard variance reduction techniques leverage model-based control variates to reduce outcome noise, they remain agnostic to potential structural relationships between competing policies. In this work, we identify a critical inefficiency in the standard A/B-testing protocol: when a treatment and control policy agree on an action, the resulting outcome contributes noise but no signal regarding the treatment effect – unnecessarily inflating confidence intervals. We propose a novel experimental protocol that exploits this policy overlap to accelerate experimentation. The key insight is to frame the randomised treatment assignment mechanism as a meta-policy, and leverage $\Delta$-Off-Policy Estimation methods to obtain unbiased estimates for average treatment effects. We prove analytically that our approach recovers standard A/B-testing practices in the general case, but that its variance scales with the divergence between policies rather than raw outcome variance. Hence, we dominate the standard Difference-in-Means estimator whenever policies have common support, and the improvement is strict whenever the overlap region contributes non-zero residual variance. Empirical results corroborate these theoretical insights – holding promise for significant impact on the real-world evaluation of recommender systems, information retrieval pipelines, and large language model interfaces.

Whether a deep residual architecture trains stably is usually determined by training it, which is expensive and answers the question only for the architecture, depth, and floating-point format tested. In practice, stability is secured by heuristics for where to place normalization and which type to use, supported by experiments but lacking a common principle. We show that stability can instead be certified before training, directly from the architectural primitives of a residual block, as an explicit function of depth and floating-point format. The certificate has two elements: (a) a power-law growth bound, $|v(x)|\le c|x|^q+b$, whose exponent $q$ is computed from the block's primitives by an arithmetic of exponents (forward); and (b) gradient bounds from a Lipschitz condition on the block over reachable states (backward). The certificate determines when the state can reach the largest finite floating-point value $M$: for $q\le1$, overflow requires at least order $\log M$ layers, whereas for $q>1$ it can occur within order $\log_q\log M$ layers; the threshold $q=1$ is sharp. Every normalization and bounded activation sets $q=0$, and the arithmetic identifies minimal modifications that bring a block from $q>1$ to $q=1$ without normalization. Empirically, in forecasting, only $q>1$ blocks blow up, as predicted by the forward certificate. In GPT-2 on OpenWebText, $q=1$ blocks without normalization also diverge: the forward certificate holds, while the backward one fails in the attention block with the largest query-key coefficient. Relaxing normalization from $q=0$ to $q=1$ improves out-of-distribution generalization in operator learning and accuracy in time-series forecasting. As deep models become more complex and costly to train, our results provide a way to assess the trade-off between stability and representational flexibility directly from block primitives, before training.

Estimating spatially heterogeneous elastic properties from low-resolution displacement measurements is a severely ill-posed inverse elasticity problem because low resolution obscures spatial details needed to distinguish heterogeneous property variations, and small measurement perturbations or fitting errors are amplified through inverse estimation. Existing inverse methods often rely on high-fidelity observations and manually prespecified loss weights, limiting their adaptability and making them sensitive to noise and resolution degradation. We propose a Probabilistic Inverse Elasticity Physics-Informed Neural Network (PIE-PINN) framework for robust estimation of Young's modulus and Poisson's ratio from noisy, low-resolution displacement data. PIE-PINN models displacement observation, strain-discrepancy, and equilibrium residuals using Laplace distributions within a unified probabilistic model. To improve robustness, the framework combines a B-spline-guided displacement network with a hierarchical half-Cauchy model for displacement residual scales. The B-spline provides a smooth global representation of the displacement field, while the neural network correction captures local variations. The hierarchical scale model adaptively downweights severe displacement fitting errors, enabling more robust recovery of the latent mean displacement field. An alternating maximum-likelihood training strategy updates the mean through weighted residual minimization and updates the scales to adjust the loss weights. Systematic case studies across varying noise levels and observation resolutions demonstrate the robustness of PIE-PINN.

Interpreting optimizers as gradient-flow discretizations has motivated applying higher-order Runge-Kutta (RK) integrators to neural networks. We build a representative Adam variant (Bogacki-Shampine 3(2) RK pair, FSAL reuse, local-error step control) and evaluate it under a strict compute-matched protocol giving every method the same gradient-evaluation budget - an accounting this literature rarely enforces. Under it the RK variant loses to plain Adam on training loss in both minibatch and full-batch (RK's best-case) training. Instrumenting it shows the "adaptivity" is illusory: normalized error stays far below tolerance, the step size pins at its growth cap from step one (98-100 percent of steps), and no rtol x hmax x h0 setting makes it act; tolerances spanning 100x give bit-identical trajectories. The method is exactly fixed-step Adam with an averaged gradient at 3-4x cost. Repairing it (true reject branch; error on the applied map) reverses the full-batch result - about 40x lower training loss than tuned Adam - and a fixed-step control isolates adaptivity (an emergent warmup-and-growth schedule) as the mechanism. But the gain is fragile to the initial step size and does not reach test accuracy. A pre-registered follow-up rules out the obvious explanations: deeper minimization does not overfit, and an explicit temperature knob only hurts - leaving a trajectory effect, the controller selecting a minimum generalizing 1.3-3.4 points below first-order descent at equal depth. An n=10 study confirms one secondary effect: gradient averaging is a genuine implicit regularizer, beating lr-matched Adam and AdamW on 10/10 seeds - yet RMSprop and NAdam match or beat it at a third the per-step cost. Higher-order adaptive integration buys deeper deterministic minimization and a small regularization effect, but nothing a cheaper, well-tuned first-order baseline does not already provide.

An $n$-point Gauss quadrature rule approximates the weighted integral of a function by a weighted sum of $n$ evaluations of this function and is exact for polynomials of degree at most $2n-1$. Such rules can be highly accurate with relatively few evaluations of this function. For weight functions associated with classical orthogonal polynomials of a continuous variable (such as Legendre, Hermite and Laguerre), these quadrature rules are readily available. We suppose that this is not the case, so that these rules must be custom-made. We present CustomGaussQuadrature, a Julia package that implements two approaches for computing the Gauss quadrature nodes and weights: the moment determinants method and the Stieltjes procedure. The principal contribution is an implementation of the ill-conditioned moment determinants method with adaptively chosen working precision. Comparisons between computations at increasing precisions provide practical error indicators used to target absolute and relative errors in the Gauss rule nodes and weights, respectively, of about $10^{-16}$. The package also adapts the working precision and number of auxiliary quadrature nodes used for the Stieltjes procedure. Numerical examples using scaled chi and Weibull probability density functions as weight functions show close agreement between the results obtained by the two approaches. The package is intended particularly for statistical applications in which a single custom Gauss rule is computed once and then reused to approximate many weighted integrals having the same weight function but different computationally expensive integrands.

Summary: In the absence of a perfect diagnostic test for a target condition, multiple imperfect tests may be used to arrive at a clinical diagnosis. Latent class analysis can be used to model such data with the objective of estimating test accuracy and target condition prevalence. Such models typically assume two latent classes - target condition positive and target condition negative. However, as we will illustrate in this manuscript, this would be an oversimplification if the different tests do not share the target condition as their measurand. We show how a Directed Acyclic Graph (DAG) can be used to illustrate the relationships between the relevant variables - the observed imperfect test results, their latent measurands, the latent target condition of interest and observed covariates - revealing any conditional dependence relations. The DAG helps determine the number of latent classes, underlying the observed data, and their labels. We show how the likelihood function changes due to incorporating the measurand of each test. We study the impact on identifiability of the model. Using simulation studies we show how ignoring the measurand of an imperfect test, when it is distinct from the target condition, can lead to biased estimates of test accuracy and prevalence. We illustrate the value of the proposed approach by re-analyzing two datasets used in previously published latent class analyses of tests for pediatric tuberculosis and leptospirosis.

We study the sample covariance error of centered Gaussians. A remarkable breakthrough [66] established the correct error scaling order and explicitly revealed the critical role of both the effective rank and the true covariance spectrum. In this work, we move beyond scaling characterizations and determine the precise limiting value of the error's spectral norm. To do so, we develop a generic framework based on Random Duality Theory (RDT). Within this framework, we first determine closed-form, explicit RDT-based upper bounds. We then establish complementary lower bounds by introducing a novel bilinear-quadratic RDT lower-bounding mechanism. By combining this mechanism with a two-replica systems bounding strategy, we show that our lower and upper bounds match in large-dimensional contexts. Our theoretical results are supplemented with numerical evaluations and simulations, demonstrating an excellent agreement already for problem sizes on the order of thousands.

Global warming and the El Niño phenomenon cause various climate changes and anomalies worldwide. Recently, extreme weather events such as heatwaves, droughts, and floods have become increasingly frequent. Researchers have studied atmospheric dynamics to reduce potential damage and increase safety. Statistical methods such as principal component analysis and maximum covariance analysis have been widely used to analyze spatial patterns of atmospheric variables. However, when the signal-to-noise ratio is low, the patterns obtained from these methods are often too noisy to be physically meaningful. This dissertation proposes regularization methods that simultaneously incorporate smoothness and sparseness penalties to estimate spatial patterns more interpretably. The proposed methods, called SpatPCA and SpatMCA, can be applied to both regularly and irregularly spaced data. An efficient algorithm based on the alternating direction method of multipliers (ADMM) is developed for computation. Through analysis of sea surface temperature data in the Indian Ocean and precipitation data in East Africa, we demonstrate the effectiveness of the proposed methods in revealing the spatial structure and studying how temperature variations in the Indian Ocean influence precipitation in East Africa.

Lipschitz continuity is a fundamental property of neural networks that characterizes their sensitivity to input perturbations. It plays a pivotal role in deep learning, governing robustness, generalization and optimization dynamics. Despite its importance, research on Lipschitz continuity is scattered across various domains, lacking a unified perspective. This paper addresses this gap by providing a systematic review of Lipschitz continuity in deep learning. We explore its theoretical foundations, estimation methods, regularization approaches, and certifiable robustness. By reviewing existing research through the lens of Lipschitz continuity, this survey serves as a comprehensive reference for researchers and practitioners seeking a deeper understanding of Lipschitz continuity and its implications in deep learning.

This paper illustrates a methodology for developing aggregated combat models using high-resolution simulation and the Markovian Lanchester process. The details of the mathematical models involved in a high-resolution simulation process of a Meeting Engagement tactical scenario are presented, along with a theoretical discussion on the Markovian Lanchester model. The output from the discrete event simulation model is used to estimate the attrition rates for an aggregated Lanchester model. A comparative study of various statistical estimation methods suitable for such estimation is also presented.

The false discovery rate (FDR) is the most widely used error metric in modern multiple testing. We provide a comprehensive analysis of the admissibility of e-value-based procedures with FDR control. We consider both simultaneous and point procedures and introduce strong and weak notions of dominance. Every simultaneous procedure is strongly, and hence weakly, dominated by an admissible weighted-mean closed e-Benjamini–Hochberg ($\overline{\mathrm{eBH}}$) procedure, and thus weighted-mean $\overline{\mathrm{eBH}}$ procedures form a complete class. Every constant-free weighted-mean $\overline{\mathrm{eBH}}$ procedure is admissible at every level, and we propose two ways for choosing constant-free weights based on side information, such as pre-screening data. Within the symmetric class, the usual mean $\overline{\mathrm{eBH}}$ procedure is the largest element if and only if the FDR level is small enough; otherwise this class has no largest element. We also obtain results on the admissibility of symmetric weighted-mean $\overline{\mathrm{eBH}}$ procedures with non-zero constant terms. Point e-testing procedures have a parallel theory of admissibility. These results highlight the central role of weighted-mean $\overline{\mathrm{eBH}}$ procedures in multiple testing.

Large Language Models (LLMs) have revolutionized AI services, but a critical tension emerges: while personalization improves model performance, it consumes scarce computational resources that users must share. When should a user invest in expensive Supervised Fine-Tuning (SFT) versus lightweight In-Context Learning (ICL)? How does congestion from other users' personalization choices reshape these incentives? And what strategies should platforms adopt when offering multiple personalization algorithms? We develop a tractable framework for LLM serving that captures the statistical-economic trade-offs users face. Our analysis yields several surprising insights. First, we show that ICL and SFT dominate in different regimes, determined by an interplay between pretraining coverage and data signal-to-noise ratios, but congestion can flip these rankings. Second, equilibrium resource consumption exhibits pronounced non-monotonicity: improving pretraining precision reduces the congestion, while broader pretraining coverage and harder tasks sometimes increase it. Third, we prove that offering both personalization methods never hurts the platform's maximal profits, despite potentially increasing computational load. Experiments with GPT-2 on linear regression tasks validate our theoretical predictions about algorithm performance. Complementing these results, our review of documentation from 21 major AI platforms shows that the share offering both SFT and ICL increased from 9.5% in 2021 to 71.4% in 2025, consistent with our platform-design implications.

We address fundamental challenges in representing and computing $\mathbb{R}^{d}$-valued predictable square-integrable processes over $[0,T]$, collected in the space $\mathcal{H}^2_T(\mathbb{R}^{d})$. These processes are central to continuous-time stochastic control, reinforcement learning, and mathematical finance. Although Wiener-chaos expansions offer strong theoretical tools, traditional computational methods are hindered by the need for large chaos dictionaries and high-order iterated integrals. To overcome these obstacles, we introduce NeuralChaos – a neural operator architecture that produces elements of $\mathcal{H}^2_T(\mathbb{R}^{d})$ using only finitely many evaluations of the driving Brownian motion, while preserving predictability and square-integrability. We prove that NeuralChaos is dense in $\mathcal{H}^2_T(\mathbb{R}^{d})$ and achieves the best $N$-term chaoslet approximation rates for compressible and Malliavin–Sobolev regular processes. Moreover, compressibility is shown to be typical for processes from $\mathcal{H}^2_T(\mathbb{R}^{d})$ under non-degenerate sub-Gaussian sampling. In contrast, we show that finite-dimensional Markovian neural SDE models constitute a meagre and Gaussian-null subset in $\mathcal{H}^2_T(\mathbb{R}^{d})$, regardless of discretization, whereas compressible processes are generic. Numerical experiments on a stochastic optimal control problem and dynamic hedging highlight the practical effectiveness of our approach. Our results enable more efficient and expressive modelling in stochastic analysis and mathematical finance.

Policy learning methods are increasingly used to inform treatment allocation under budget constraints. Most proposed methods assume complete treatment data, yet applications frequently suffer from missingness that can bias estimates and lead to suboptimal policies. We address this gap by extending efficient estimators for average treatment effect (ATE) estimation to policy value and conditional average treatment effect (CATE) estimation under missing at random (MAR) and missing completely conditionally at random (MCCAR) treatment data. Through asymptotic efficiency analysis, we prove that the MAR estimator, which leverages partially-observed units, is both valid and more efficient than the MCCAR estimator when MCCAR assumptions hold. This result provides formal justification for preferring MAR-based estimation in policy learning under both missing data settings. Our comprehensive experiments using synthetic and semi-synthetic datasets confirm that correctly specifying the missingness mechanism is crucial: misspecified estimators remain biased regardless of sample size, while our estimators achieve near-oracle performance when assumptions are satisfied. Our work provides practitioners with theoretically grounded, empirically validated tools for robust policy learning in the presence of missing treatment data.

We propose a density-based leave-one-out influence score for unsupervised outlier detection. The motivation is that outliers are naturally associated with regions of very small probability density, but direct leave-one-out density refitting can be computationally prohibitive. We use the Linear-Blend Frequency Polygon (LBFP) estimator and define a score that compares the full-sample fitted density at an observation with the fitted density obtained after removing that observation, while keeping the grid and bandwidth fixed. The resulting statistic measures a relative density perturbation at the observation's own location. For the LBFP estimator, this score has an exact closed-form update, so the density estimator does not need to be refitted for each observation. This preserves a direct density interpretation while making the method computationally efficient for large samples. We study the score under contamination and show that regular positive-density observations and contamination-driven observations have distinct asymptotic orders. Simulations over a broad range of contamination models illustrate these theoretical regimes, show competitive performance relative to standard benchmarks, and document computing time. A credit-card fraud application with 29 variables illustrates that the method works well on a large real data set.

We study sparse threshold random geometric graphs generated by high-dimensional spherical or Gaussian latent vectors. Although each edge has marginal probability $p$, shared latent variables make the adjacency entries dependent. At the connectivity scale $np=\Omega(\log n)$, the spherical adjacency matrix satisfies, with high probability,$\|A-\mathbb E A\|_{\mathrm{op}}=O\left(\sqrt{np\log n}+np\tau\right)$, where $\tau$ is the cap threshold; an analogous estimate holds for Gaussian vectors after controlling radial fluctuations. This sharpens the spectral bound in Liu, Mohanty, Schramm, and Yang (2023) under weaker assumptions and strengthens the global-synchronization guarantee of Abdalla, Bandeira, and Invernizzi (2024) for the homogeneous Kuramoto model. The leading eigenspace also estimates the latent geometry. When $np\gg\log n$, vector and relative Gram-matrix errors vanish for$\log(1/p)\ll d\ll np\log(1/p)/\log n$ in the spherical model and $\log^2(1/p)\log n\ll d\ll np\log(1/p)/\log n$ in the Gaussian model, improving the recovery conditions of Li and Schramm (2023). For the Gaussian mixture block model introduced there, a polynomial-time semidefinite program gives, to our knowledge, the first exact-recovery guarantee at the connectivity scale in a moderate-separation regime. At much larger separation, fixed edge density creates isolated vertices and makes exact recovery impossible. Our reusable decoupling and matrix concentration framework avoids trace-moment methods and applies broadly to random graph models with latent vectors.

Mixture models which cluster skewed random matrices can often suffer from over-parameterization in the absence of performing dimension reduction. Even with the use of bilinear factor analyzers, further parameter reduction can be achieved by constraining parameters over clusters. In this manuscript propose a parsimonious family of 256 models for mixtures of skewed matrix variate bilinear factor analyzers, specifically in the case of the skew t distribution. An AECM algorithm for parameter estimation is discussed in detail. Further, extensive simulations are performed, and the method is considered in the case of the MNIST dataset and the Olivetti faces dataset.