Stochastic simulation aims to compute output performance for complex models that lack analytical tractability. To ensure accurate prediction, the model needs to be calibrated and validated against real data. Conventional methods approach these tasks by assessing the model-data match via simple hypothesis tests or distance minimization in an ad hoc fashion, but they can encounter challenges arising from non-identifiability and high dimensionality. In this article, we investigate a framework to develop calibration schemes that satisfies rigorous frequentist statistical guarantees, via a basic notion that we call eligibility set designed to bypass non-identifiability via a set-based estimation. We investigate a feature extraction-then-aggregation approach to construct these sets that target at multivariate outputs. We demonstrate our methodology on several numerical examples, including an application to calibration of a limit order book market simulator (ABIDES).

The bootstrap is a versatile method for quantifying statistical uncertainty. Among its variants, a popular approach, the studentized bootstrap, provably achieves higher-order coverage error reduction compared to other benchmarks. However, its implementation typically requires an analytical form of the standard error, or otherwise an additional layer of resampling effort which can be computationally expensive. In this paper, we introduce what we call the studentized cheap bootstrap that achieves the same higher-order coverage accuracy as the conventional studentization, but substantially thinning the computational effort in the additional resampling layer to only very few Monte Carlo replications. Intriguingly, while conventional wisdom views "studentization" as an informal link between the bootstrap and t-distribution, we provide a first recognition that this link is in fact formal, notably with a distinct insight that the degree of freedom in the t-distribution corresponds to the Monte Carlo computation effort in the additional resampling layer, rather than the data size as in traditional thinking. Moreover, our desirable higher-order coverage accuracy builds crucially on this insight, as well as explicit calculations and geometric analyses of higher-order terms in the Edgeworth and Cornish-Fisher expansions tailored to limiting t-distributions.

We consider stochastic optimization where the goal is not only to optimize an average-case objective, but also to mitigate the occurrence of rare catastrophic events. This problem is motivated by safety-aware decision-making and AI training. We first argue that, in the presence of a simulation model, natural attempts to integrate variance reduction into optimization, even executed in a reasonable adaptive fashion, encounter fundamental challenges in guaranteeing realistic runtime when using common stochastic gradient descent algorithms. This challenge arises from the extreme sensitivity of tail-based objectives with respect to the decision variables, which renders a dichotomic failure of convergence regardless of what step size we select. We offer remedies based on a new notion of safe start that allows for efficient finite-time error control, and show how the sampling complexity scales favorably under the combination of safe start and variance reduction. We illustrate our methodologies on examples in portfolio optimization and robust classification with neural networks.

Distributionally robust optimization (DRO) is a worst-case framework for stochastic optimization under uncertainty that has drawn fast-growing studies in recent years. When the underlying probability distribution is unknown and observed from data, DRO suggests computing the worst-case distribution within a so-called uncertainty set that captures the involved statistical uncertainty. In particular, DRO with uncertainty set constructed as a statistical divergence neighborhood ball has been shown to provide a tool for constructing valid confidence intervals for nonparametric functionals and bears a duality with the empirical likelihood (EL). In this paper, we show how adjusting the ball size of such type of DRO can reduce higher-order coverage errors similar to the so-called Bartlett correction. Our correction, which applies to general von Mises differentiable functionals, is more general than the existing EL literature that only focuses on smooth function models or M-estimation. Moreover, we demonstrate a higher-order "self-normalizing" property of DRO regardless of the choice of divergence. Our approach builds on the development of a higher-order expansion of DRO, which is obtained through an asymptotic analysis on a fixed-point equation arising from the Karush-Kuhn-Tucker conditions.

Modeling multivariate distributions with nonlinear dependence, multimodality, and tractable analytical structure for downstream applications is a central challenge in uncertainty quantification. Sliced Normal (SN) distributions were introduced in prior works at the National Aeronautics and Space Administration (NASA) to address this need by representing densities through polynomial feature maps. This construction provides a compact algebraic alternative to more opaque generative models, while retaining the ability to capture nonlinear parameter dependencies and multi-modal behavior. In this paper, we build on the SN framework and develop several improvements that make the approach more reliable and scalable. First, we reformulate SN parameter estimation as a convex optimization problem over a positive semidefinite matrix, replacing the original nonconvex likelihood search with a formulation amenable to standard optimization tools. Second, we clarify the expressive power of the SN class by connecting polynomial log-density modeling to a Stone–Weierstrass-type universal approximation argument on compact domains. Third, we propose a high-dimensional fitting procedure that partitions variables into approximately independent groups, fits SN models within each subgroup, and then assembles the subgroup models through a cross-block completion step to recover residual dependence. We demonstrate the resulting SN modeling pipeline on NASA loss-of-control flight data, where the method captures nonlinear dependence patterns in both low-dimensional slices and a higher-dimensional block-assembled model.

We study stochastic gradient estimation in black-box environments where only noisy simulation observations of function values are available. Finite-difference (FD) methods are among the most widely used zeroth-order gradient estimators in such settings, by measuring the change in function values against a perturbation size. While the optimal order in choosing this perturbation size with respect to the simulation budget is well understood, the optimal constant factor relies on model characteristics that are typically unknown and viewed to be as difficult to estimate as the gradient itself. Consequently, FD estimators are often based on ad hoc tuning of the perturbation size, which may exhibit highly unstable performance across problem instances. In this paper, we challenge this conventional wisdom from both theoretical and practical perspectives. We show that, by pilot-estimating these model quantities using a negligible fraction of the simulation budget, substantial robustness is attained in the resulting FD estimators. Theoretically, we show that using a perturbation size governed by this pilot estimation can already achieve an MSE that is first-order identical to the ``oracle" MSE as if the optimal perturbation size is known in advance. Moreover, we show how such an approach is competitive against any choices of prescribed perturbation size, even if they are designed to be minimax-optimal over reasonable classes of target functions and FD schemes. Our proposed pilot estimation is practically easy to run, and a variety of numerical experiments demonstrate both the robustness and near-oracle optimality of our estimator relative to conventional FD schemes based on ad hoc tuning.

Conformal prediction is a popular method to construct prediction intervals with marginal coverage guarantees from black-box machine learning models. In applications with potentially high-impact events, such as flooding or financial crises, regulators often require very high confidence for such intervals. However, if the desired level of confidence is too large relative to the amount of data used for calibration, then classical conformal methods provide infinitely wide, thus, uninformative prediction intervals. In this paper, we propose a new method to overcome this limitation. We bridge extreme value statistics and conformal prediction to provide reliable and informative prediction intervals with high-confidence coverage, which can be constructed using any black-box extreme quantile regression method. A weighted version of our approach can account for nonstationary data. The advantages of our extreme conformal prediction method are illustrated in a simulation study and in an application to flood risk forecasting.

Recent proliferation of data-optimization integration has led to a range of methods that aim to improve the statistical performance of data-driven optimization decisions. However, while many of these methods are motivated intuitively from a robustness or regularization perspective, their resulting statistical benefits are often unclear and, even if available, are established on a case-by-case basis. We provide a systematic dissection of data-driven optimization formulations using the view of "directionally perturbed" empirical optimization (EO). Specifically, this umbrella of formulations, which we call "EO+", covers many existing data-driven optimization methods, including regularization, distributionally robust optimization, transfer learning, and analogous methods for contextual optimization. On the one hand, we argue that without additional, correctly specified, side information, any EO+ method can result in at most second-order improvements. This provides a negative conclusion, namely “no free lunch is possible", on the statistical power of EO+. On the other hand, we show that when leveraging side information that is geometrically effective, achieving first-order improvements is possible by choosing hyperparameters that are significantly larger than what is typically suggested in the literature. Moreover, we construct a principled methodology based on excess risk estimation, via either system knowledge or bootstrap resampling, to maximize the first-order gain. We demonstrate how this gain connects to the control-variate principle, a variance reduction technique in the Monte Carlo simulation literature, which helps explain why geometrically effective side information is necessary.

Constrained simulation optimization (CSO) is a general framework for optimizing stochastic systems under performance constraints. It arises widely in practice where objective and constraint evaluations are available only through noisy simulation outputs. Compared with the unconstrained setting, the lack of accessible analytical gradients for simulation-based constraints makes it more challenging to develop efficient solution methods and establish non-asymptotic guarantees. To address this gap, we propose a novel single-loop algorithm, called min-max gradient search (MGS), which integrates a primal-dual framework with stochastic gradient estimators. Unlike conventional stochastic approximation methods based on gradient descent for solving simulation optimization problems, such as Zhou and Bhatnagar (2017) and Hu and Fu (2025), MGS performs alternating gradient descent and ascent on the primal and dual variables, which improves the objective while penalizing constraint violations. For the first time, we establish a finite-time convergence guarantee for single-loop CSO algorithms by showing that MGS converges to a stationary solution (a Karush-Kuhn-Tucker point under mild conditions) at a rate of Õ(T^-1/3), where T is the number of iterations. Numerical experiments on a serial queuing system and a 2000-dimensional optimization problem demonstrate the superior performance and scalability of MGS.

Direct Preference Optimization (DPO) has recently emerged as a popular approach to improve reinforcement learning from human feedback (RLHF), leading to better techniques to fine-tune large language models (LLM). A weakness of DPO, however, lies in its lack of capability to characterize the diversity of human preferences. Inspired by Mallows' theory of preference ranking, we develop in this paper a new approach, the *MallowsPO*. A distinct feature of this approach is a *dispersion index*, which reflects the dispersion of human preference to prompts. We show that existing DPO models can be reduced to special cases of this dispersion index, thus unified with MallowsPO. More importantly, we demonstrate empirically how to use this dispersion index to enhance the performance of DPO in a broad array of benchmark tasks, from synthetic bandit selection to controllable generation and dialogues, while maintaining great generalization capabilities. MallowsPO is also compatible with other SOTA offline preference optimization methods, boosting nearly 2\% extra LC win rate when used as a plugin for fine-tuning Llama3-Instruct.

Bagging has emerged as an effective tool for reducing variance and enhancing stability in model training, via repeated data resampling followed by a suitable aggregation. Recently, it has also been used to obtain performance bounds for data-driven solutions in stochastic optimization. However, quantifying statistical uncertainty for bagged estimators can be challenging, as standard bootstrap would require resampling at both the bagging and the bootstrap stages-leading to multiplicative computation costs that can be prohibitively large. In this work, we propose a practical and theoretically justified approach using the cheap bootstrap methodology, which enables valid confidence interval construction for bagged estimators under a controllable number of model evaluation. We establish asymptotic validity of our approach and demonstrate its empirical performance through simulation experiments. Our results show that the proposed method achieves nominal coverages with significantly reduced computational burden than other benchmarks.

Importance Sampling (IS) is a widely used variance reduction technique for enhancing the efficiency of Monte Carlo methods, particularly in rare-event simulation and related applications. Despite its effectiveness, the performance of IS is highly sensitive to the choice of the proposal distribution and often requires stochastic calibration. While the design and analysis of IS have been extensively studied in estimation settings, applying IS within stochastic optimization introduces a fundamental challenge: the decision variable and the importance sampling distribution are mutually dependent, creating a circular optimization structure. This interdependence complicates both convergence analysis and variance control. We consider convex stochastic optimization problems with linear constraints and propose a single-loop stochastic approximation algorithm, based on a joint variant of Nesterov's dual averaging, that jointly updates the decision variable and the importance sampling distribution, without time-scale separation or nested optimization. The method is globally convergent and achieves minimal asymptotic variance among stochastic gradient schemes, matching the performance of an oracle sampler adapted to the optimal solution.

Estimating the false discovery rate (FDR) is one of the key steps in ensuring appropriate error control in the analysis of shotgun proteomics data. Traditional estimation methods typically rely on decoy sequence databases or spectral libraries, which may not always provide satisfactory results due to limitations of decoy construction methods. This study introduces the query mix-max (QMM) method, a decoy-free alternative for FDR estimation in proteomics. The QMM framework builds upon the existing mix-max procedure but replaces decoy matches with entrapment queries to estimate the number of false positive discoveries. Through simulations and real data set analyses, the QMM method was demonstrated to provide reasonably accurate FDR estimation across various scenarios, particularly when smaller sample-to-entrapment spectra ratios were achieved. The QMM method tends to be conservatively biased, particularly at higher FDR values, which can ensure stringent FDR control. While flexible, the protocol's effectiveness may vary depending on the evolutionary distance between the sample and entrapment organisms. It also requires a sufficient number of entrapment queries to provide stable FDR estimates, especially for low FDR values. Despite these limitations, the QMM method is a promising alternative as one of the first query-based FDR estimation approaches in shotgun proteomics.

This paper provides an overview of black-box rare-event simulation methods applicable to the safety testing of artificial intelligence agents. We explore the challenges and efficiency criteria in black-box simulation, especially emphasizing the subtle occurrence and control of underestimation errors. The paper reviews various adaptive methods, such as the cross-entropy method and adaptive multilevel splitting, highlighting both their empirical effectiveness and theoretical limitations. Additionally, it offers a comparative analysis of different confidence interval constructions for crude Monte Carlo methods, aiming to mitigate underestimation errors through effective uncertainty quantification. The paper concludes with a certifiable deep importance sampling approach, using deep neural networks to develop conservative estimators that address underestimation issues.

Data-driven stochastic optimization is ubiquitous in machine learning and operational decision-making problems. Sample average approximation (SAA) and model-based approaches such as estimate-then-optimize (ETO) or integrated estimation-optimization (IEO) are all popular, with model-based approaches being able to circumvent some of the issues with SAA in complex context-dependent problems. Yet the relative performance of these methods is poorly understood, with most results confined to the dichotomous cases of the model-based approach being either well-specified or misspecified. We develop the first results that allow for a more granular analysis of the relative performance of these methods under a local misspecification setting, which models the scenario where the model-based approach is nearly well-specified. By leveraging tools from contiguity theory in statistics, we show that there is a bias-variance tradeoff between SAA, IEO, and ETO under local misspecification, and that the relative importance of the bias and the variance depends on the degree of local misspecification. Moreover, we derive explicit expressions for the decision bias, which allows us to characterize (un)impactful misspecification directions, and provide further geometric understanding of the variance.

The field of simulation optimization (SO) encompasses various methods developed to optimize complex, expensive-to-sample stochastic systems. Established methods include, but are not limited to, ranking-and-selection for finite alternatives and surrogate-based methods for continuous domains, with broad applications in engineering and operations management. The recent advent of large language models (LLMs) offers a new paradigm for exploiting system structure and automating the strategic selection and composition of these established SO methods into a tailored optimization procedure. This work introduces SOCRATES (Simulation Optimization with Correlated Replicas and Adaptive Trajectory Evaluations), a novel two-stage procedure that leverages LLMs to automate the design of tailored SO algorithms. The first stage constructs an ensemble of digital replicas of the real system. An LLM is employed to implement causal discovery from a textual description of the system, generating a structural `skeleton' that guides the sample-efficient learning of the replicas. In the second stage, this replica ensemble is used as an inexpensive testbed to evaluate a set of baseline SO algorithms. An LLM then acts as a meta-optimizer, analyzing the performance trajectories of these algorithms to iteratively revise and compose a final, hybrid optimization schedule. This schedule is designed to be adaptive, with the ability to be updated during the final execution on the real system when the optimization performance deviates from expectations. By integrating LLM-driven reasoning with LLM-assisted trajectory-aware meta-optimization, SOCRATES creates an effective and sample-efficient solution for complex SO optimization problems.

When the underlying probability distribution in a stochastic optimization is observed only through data, various data-driven formulations have been studied to obtain approximate optimal solutions. We show that no such formulations can, in a sense, theoretically improve the statistical quality of the solution obtained from empirical optimization. We argue this by proving that the first-order behavior of the optimality gap against the oracle best solution, which includes both the bias and variance, for any data-driven solution is second-order stochastically dominated by empirical optimization, as long as suitable smoothness holds with respect to the underlying distribution. We demonstrate this impossibility of improvement in a range of examples including regularized optimization, distributionally robust optimization, parametric optimization and Bayesian generalizations. We also discuss the connections of our results to semiparametric statistical inference and other perspectives in the data-driven optimization literature.

Data-driven optimization aims to translate a machine learning model into decision-making by optimizing decisions on estimated costs. Such a pipeline can be conducted by fitting a distributional model which is then plugged into the target optimization problem. While this fitting can utilize traditional methods such as maximum likelihood, a more recent approach uses estimation-optimization integration that minimizes decision error instead of estimation error. Although intuitive, the statistical benefit of the latter approach is not well understood yet is important to guide the prescriptive usage of machine learning. In this paper, we dissect the performance comparisons between these approaches in terms of the amount of model misspecification. In particular, we show how the integrated approach offers a “universal double benefit” on the top two dominating terms of regret when the underlying model is misspecified, while the traditional approach can be advantageous when the model is nearly well-specified. Our comparison is powered by finite-sample tail regret bounds that are derived via new higher-order expansions of regrets and the leveraging of a recent Berry-Esseen theorem.

In this study, we systematically evaluated the nanoLC-Zeno TOF platform using data-independent acquisition (DIA) and identified 3,142 to 7,059 protein groups from 1 to 200 ng of standard K562 digest. Additionally, a targeted acquisition strategy detected 184 peptides with enhanced signal using Zeno activation. To achieve spatial cell-type resolved proteomics, we developed a workflow integrating image-guided laser capture microdissection for analysis of gastric cancer tissue. A 30-minute nanoLC gradient facilitated the quantification of 2,514 proteins from 0.2 mm² regions representing cancer cells, cancer-associated fibroblasts (CAFs), and immune cells, as defined by multiplex immunohistochemistry. This spatial proteomics approach revealed cellular heterogeneity within the tumor microenvironment of a clinical gastric cancer sample, and targeted validation identified differential expression of 12 cancer-associated proteins comprising 103 peptides between cancer cell-enriched and CAFs-enriched regions. This workflow offers a practical strategy for spatial visual proteomic analysis in clinical tissue samples. ### Competing Interest Statement The authors have declared no competing interest. China State Key Basic Research Program Grants (2024YFA1307200, 2021YFA1301601, 2020YFE0202200, 2022YFC3401104, 2021YFA1301602 and 2021YFA1302603), the National Natural Science Foundation of China (92253304, 22125403, 32201218, and 22104047), the Shenzhen Innovation of Science and Technology Commission (JSGGZD20220822095200001, JCYJ20200109141212325, JCYJ20210324120210029 and JCYJ20200109140814408). Research Grants Council, Hong Kong SAR Government (Grant No. C5005-23W). China Postdoctoral Science Foundation (2021M701410), Shenzhen Medical Research Fund (2401008).

Ensemble learning is a popular technique to improve the accuracy of machine learning models. It traditionally hinges on the rationale that aggregating multiple weak models can lead to better models with lower variance and hence higher stability, especially for discontinuous base learners. In this paper, we provide a new perspective on ensembling. By selecting the most frequently generated model from the base learner when repeatedly applied to subsamples, we can attain exponentially decaying tails for the excess risk, even if the base learner suffers from slow (i.e., polynomial) decay rates. This tail enhancement power of ensembling applies to base learners that have reasonable predictive power to begin with and is stronger than variance reduction in the sense of exhibiting rate improvement. We demonstrate how our ensemble methods can substantially improve out-of-sample performances in a range of numerical examples involving heavy-tailed data or intrinsically slow rates.